Genome of lethal Lepiota venenata and insights into the evolution of toxin-biosynthetic genes

Background Genomes of lethal Amanita and Galerina mushrooms have gradually become available in the past ten years; in contrast the other known amanitin-producing genus, Lepiota, is still vacant in this aspect. A fatal mushroom poisoning case in China has led to acquisition of fresh L. venenata fruiting bodies, based on which a draft genome was obtained through PacBio and Illumina sequencing platforms. Toxin-biosynthetic MSDIN family and Porlyl oligopeptidase B (POPB) genes were mined from the genome and used for phylogenetic and statistical studies to gain insights into the evolution of the biosynthetic pathway. Results The analysis of the genome data illustrated that only one MSDIN, named LvAMA1, exits in the genome, along with a POPB gene. No POPA homolog was identified by direct homology searching, however, one additional POP gene, named LvPOPC, was cloned and the gene structure determined. Similar to ApAMA1 in A. phalloides and GmAMA1 in G. marginata, LvAMA1 directly encodes α-amanitin. The two toxin genes were mapped to the draft genome, and the structures analyzed. Furthermore, phylogenetic and statistical analyses were conducted to study the evolution history of the POPB genes. Compared to our previous report, the phylogenetic trees unambiguously showed that a monophyletic POPB lineage clearly conflicted with the species phylogeny. In contrast, phylogeny of POPA genes resembled the species phylogeny. Topology and divergence tests showed that the POPB lineage was robust and these genes exhibited significantly shorter genetic distances than those of the house-keeping rbp2, a characteristic feature of genes with horizontal gene transfer (HGT) background. Consistently, same scenario applied to the only MSDIN, LvAMA1, in the genome. Conclusions To the best of our knowledge, this is the first reported genome of Lepiota. The analyses of the toxin genes indicate that the cyclic peptides are synthesized through a ribosomal mechanism. The toxin genes, LvAMA1 and LvPOPB, are not in the vicinity of each other. Phylogenetic and evolutionary studies suggest that HGT is the underlining cause for the occurrence of POPB and MSDIN in Amanita, Galerina and Lepiota, which are allocated in three distantly-related families. Electronic supplementary material The online version of this article (10.1186/s12864-019-5575-7) contains supplementary material, which is available to authorized users.

Whether Conocybe can produce similar cyclic peptide toxins is still in question [14,15]. The bicyclic toxin is highly stable and resistant to high temperature, acids, alkalis and salts; general cooking methods do not destroy its toxicity [16]. Amanitins, including α-Amanitin, are synthesized through a ribosomal mechanism in Amanita and Galerina [17,18]. However, the pathway of amanitin biosynthesis in Lepiota is unknown. With the advent of the genome era, a few genomes of lethal mushrooms containing α-amanitin have been published [17][18][19][20], but no genome of Lepiota species has been reported. In contrast to amanitin-producing agarics, genomes of other Amanita species including A. muscaria [21], A. jacksonii [22] and A. thiersii [23] have become available in recent years. These non-amanitin-producing Amanita species do not produce cyclic peptides and our extensive BLAST search in those genomes did not return any known toxin-biosynthetic genes, namely the MSDIN family and prolyl oligopeptidase B (POPB) genes.
In recent years, there have been many poisoning tragedies all over the world caused by accidental consumptions of poisonous mushrooms containing α-amanitin. It was reported in May 2018 that serious mushroom poisoning incidents occurred in Iran in a short period, which resulted in 18 deaths out of 1151 patients (http:// www.tehrantimes.com/news/423947/Mushroom-poisoning-kills-18-in-Iran).
Amanita virosa, a α-amanitin-containing species commonly known in Europe as the destroying angel, is likely responsible for this poisoning case. In lethal Amanita species, a mature individual may contain α-amanitin that exceeds the adult lethal dose [15,24], and 144 deaths were caused by Amanita species in China during 1994-2016 [11]. In September 2017, two individuals died of L. venenata poisoning in Jingzhou City, Hubei Province of Central China [25]. After the incidence, we collected the fresh fruiting bodies of L. venenata from the locality, and these mushrooms were sent for genome sequencing with combined platform of PacBio Sequel and Illumina HiSeq X10. The draft genome was in turn used to study MSDIN and POPB genes that associate with amanitin biosynthesis. The biosynthesis of this toxin and related cyclic peptides begins with activation of MSDIN genes that encode a precursor peptide of 34-37 amino acids, with the first 5 highly conserved residues mostly being MSDIN [17]. The precursor peptides are cleaved and macrocyclized into 7-10 amino acid cyclic peptides by a specialized prolyl oligopeptidase, POPB [19,26].
Our previous report found that the three disjunct genera, Amanita, Galerina and Lepiota, all possess a similar biosynthetic pathway for amanitin production, which is probably due to horizontal gene transfer (HGT) [14]. At the time, the phylogenetic approach did not fully rule out massive gene loss (over 2000 genes counted) as one of the possibilities. In this report, MSDIN and POP (POPA and POPB) genes were investigated based on a draft genome of the lethal L. venenata. The genes were mined and the structures determined. Using combined approaches, we found a new homolog of POP genes. With this new POP gene, we once again constructed phylogenetic trees of POP genes for more insights into the evolutionary history. Subsequently, the topology was examined, and the POPB branch tested for reliability. In this report, we mainly searched for conflict between POPB and species phylogenies, as this approach is considered as the most reliable way of assessing HGT. Other supporting methods, including genetic distance (applied to both POPB and MSDIN), topology test, and comparison of gene and species trees, were included to further support our conclusion.

Results
Overview of L. venenata genome and toxin genes Genome sequencing of L. venenata was performed using single molecule real-time (SMRT) DNA sequencing. A draft genome sequence was generated on the PacBio Sequel platforms and Illumina HiSeq X10 (GenBank: RCFS00000000 Ver. RCFS01000000). The total amount of clean data was 4.23 Gb, and the size of assembled genome was 49.25 Mb. The main genomic features were shown in the table (Table 1). BUSCO prediction with Basidiomycota settings implied that 89.5% genes were found, and missing BUSCOs were at low level of 5.0%. Only one MSDIN gene was found after both tBLASTn and BLASTn searches using multiple queries from both Amanita and Galerina. This gene was located on contig 27, and in total we obtained 571 bp in length with 154 bp from ATG to TAA. The deduced amino acid sequence had 32 amino acids, and the 10th to the 17th amino acids, IWGIGCNP, coded for α-amanitin. The full amino acid sequence was MDANATRLPIW-GIGCNPWTPESVNDTLTKDLS with the core peptide in bold and underlined. Due to the similarity to other AMA1 genes, it was named LvAMA1 (Additional file 1). The POPB gene (denoted as LvPOPB later) was found on contig 4 with 3118 bp in length (ATG to TGA) (Additional file 2). The genome was visualized as a circular diagram using Circos software (Fig. 1). The tracks (a to f ) represented the contigs (a), GC content (b), GC skew (c), the location of predicted genes (d), the location of exons (e), and the distribution of the two toxin genes in the genome (f ). As can be seen from the figure, the predicted distribution of genes and exons in the whole genome were relatively uniform. LvPOPB and LvAMA1 were located on contigs 4 and 27 respectively, and were far apart (at least 525 Kb) with over 100 predicted genes in between. Clearly, these two toxin genes did not form a cluster.

Analysis of toxin(s) by LC-MS
Only one cyclic peptide, α-amanitin, was found in the extract (Fig. 2 ] ion is 920.3620, which was not detected in the extract. The LC analyses detected a peak of 4.76 min, presenting the expected mass of protonized α-amanitin, and no beta-amanitin was found throughout the elution range (0-18 min). The masses of two other major toxins, phalloidin and phallacidin, were also applied in the scan, and no corresponding masses were identified. The result indicated that only one toxin, This result showed a different toxin profile compared to previous studies in other Lepiota species. Usually lethal Lepiota contains more than one toxin, for example, L. brunneoincarnata has a certain amount of both α-amanitin and β-amanitin [27], and αand γ-amanitins are present in L. josserandii [6].

Predicted structures of LvAMA1 and LvPOPB genes
The gene structures of LvAMA1 and LvPOPB were predicted on Splign website (Fig. 3). The LvAMA1 gene contained four exons and three introns. The coding sequence spanned the first intron, consisting of 99 bp. The LvPOPB gene was composed of 18 exons and 17 introns. The LvAMA1 protein had 32 amino acids, and the LvPOPB had 731 amino acids. Over all, the LvPOPB is very similar to other know POPB genes, in structure and in sequence.

Phylogenetic analysis of POP
POP phylogenetic tree was constructed with three outgroups (Ophiocordyceps sinensis, Marssonina brunnea and Colletotrichum nymphaeae). As shown in Figs. 4 and 5, whether from the nucleotide or the amino acid phylogenetic tree, POPB genes were clustered together forming a highly monophyletic clade with strong support. However, POPA distribution was more complicated: POPAs of Amanita clustered together as expected, while that of G. marginata was separated from them nesting in a group containing Gymnopilus chrysopellus and Cortinarius glaucopus. "POPA", referred as "generic POP" from now on, is considered as a house-keeping gene present in most basidiomycetes, which reflects the widely accepted phylogeny of the species included in the analyses [28][29][30]. Based on the result, "POPA" was used later in this report to refer to only Amanita POPAs.
A new POP gene was found in the genome of L. venenata, but surprisingly, it was not closely related to either POPB or POPA genes. As shown in Figs. 4 and 5, it was clustered with the POP analogs of Laccaria bicolor and Macrolepiota fuliginosa. We named it LvPOPC for reference purposes. The DNA and CDS sequences of LvPOPC were obtained through cloning and the detailed structure was shown in Fig. 6, which was compared with of LvPOPC sequence was the shortest while AbPOPA was the longest. One of the most significant differences was that the last exon length of AbPOPA was roughly twice the size of those in AbPOPB and LvPOPC. On the other hand, the pattern of exon-intron structure of LvPOPC was similar to that of AbPOPA and other POPAs, with 19 exons and 18 introns in similar arrangement, but significantly shorter. On the level of sequence similarity, LvPOPC was more divergent from both POPA and POPB, as it showed significantly lower identity during BLAST search. This explains why LvPOPC did not cluster with either POPA or POPB.
According to the nucleotide phylogenetic tree of POPs (Fig. 4), the taxa shaded in blue indicated the positions of generic POPs (POPAs, GmPOP and LvPOPC) in the three genera that produce α-amanitin. These phylogenetic positions were consistent with the typical species tree for Agaricales constructed by combined rpb1, rpb1-intron2, rpb2, 18S, 25S and 5.8S nucleotide sequences [30]. The results showed that generic POPs conformed to the agaric phylogeny. The taxa marked in red were the monophyletic branch of POPB with strong supportive statistics (Maximum likelihood bootstraps: 100% and Bayesian posterior probabilities: 1). The blue bars on the right of the figure indicated that all the taxa included in the analysis were consistent with the species phylogeny, except the POPB lineage (red bar). The POPBs of three different families violated the species phylogenetic relationships and were clustered into one monophyletic branch. Previous studies have shown that evolutionary tree analysis can be used as a powerful method to assess whether genes occur through HGT [31][32][33]. If the evolution of one gene with high sequence similarity in different species does not conform to the phylogenetic relationship, the gene may have HGT [34][35][36]. Our results showed that POPBs from three agaric families had high similarities in sequence and structure, but their phylogeny did not conform to the species phylogeny.
In order to further analyze the congruency between POPB lineage and the species phylogeny, we compared the lineage with a species tree based on rpb2 (Fig. 7). In POPB lineage, the branches of Galerina and Amanita formed a sister group and the Lepiota species produced a monophyletic branch, clearly contradicting the species phylogeny indicated in the species tree, where Galerina and Lepiota were clustered together. Unlike our previous report, in this study the topology incongruency of the gene tree and the species tree is clear as described in the following, largely due to the newly found gene LvPOPC. POPA (from Amanita species), GmPOP and LvPOPC precisely reflected their phylogenetic positions in the widely accepted agaric phylogeny [30], and when massive gene loss is considered (i.e., species in between the genera removed from analysis), the collapsed

Topology test
To verify the robustness of the POPB branch, we used PAUP program to establish three hypothetical trees (see Methods), calculated their site-wise log-likelihoods, and then conducted approximately unbiased test with Consel. The results (Table 2) showed that, the approximately unbiased p-values (AU) and bootstrap probability (NP) of the three hypothetical trees were close to 0, while the value of best tree (the CDS POP phylogenetic tree) was 1.
From the statistical point of view, the probability of occurrence of the best tree was 100%, while the three hypothetical trees were rejected. The evidence suggested that the POPB branch was robust and reliable.

Gene tree and species tree
In order to evaluate evolutionary history of POPB, we established a species tree and a gene tree similar to those in our previous report [14]. The species tree faithfully reflected the phylogenetic relationship of Amanita, Galerina and Lepiota. For the gene tree, the branches of POPA and POPB were clustered with their own homologs, respectively. With Notung, the Divergence-Loss models (DL model) returned the following general statistics: Event Score = 38.5, Dups = 5, Losses = 31, and Numbers of optimal solutions = 1.  The Divergence-Transfer-Loss (DTL model) produced: Event Score = 24.0, Dups = 0, Transfers = 5, and Losses = 9. The DTL Event score was significantly smaller than that of DL model. This indicated that the probability of DTL model with HGT in POPB was significantly higher than that in gene duplication and loss only assumption. In addition, a POP and a MSDIN gene trees similar to those in our recent report [14] were constructed. These trees were analyzed by PAML software. The species tree represented the phylogeny of the housekeeping gene (rpb2), while the gene trees represented the evolutionary relationships of the POP and MSDIN genes. The results showed that there were significant differences in distances, synonymous rates (dS) and nonsynonymous rates (dN) among the three amanitin-producing genera. The distances of rpb2 was significantly higher than those of POPB (2~6 times, Table 3) and MSDIN (up to 3 times, Table 3). The data in the table also showed that the dN/ dS values of POPB were small, and the dS ratio for A. phalloides and G. marginata was about 1:7, while the dN ratio was close to 2:1. These results demonstrated that the distances of POPB and MSDIN genes were significantly shorter than those of rpb2, and therefore the evolution of POPB and MSDIN was not in accordance with vertical inheritance. These data are consistent with HGT scenario as significantly lower rates are expected when compared to house-keeping genes.

Discussion
Lepiota genome α-Amanitin, the major cyclic peptide toxin, is responsible for the vast majority (> 90%) of deadly mushroom poisonings worldwide [37]. The organisms containing this toxin are known to be distributed in the three disjunct agaric genera, Amanita, Galerina and Lepiota [4,[38][39][40][41]. Draft genomes are available for both Amanita and Galerina, but not for Lepiota due to the difficulty of obtaining high quality material. In China, there has been more than a few deadly poisoning cases caused by Lepiota in the past three years [25], allowing acquisition of samples suitable for genome sequencing. To the best of our knowledge, this is the first available genome of Lepiota, which now allows studies of the toxin biosynthesis on full spectrum, i.e., including all known α-amanitin-containing genera. Research on the mechanism of the toxin biosynthesis and toxin genes of lethal Lepiota species also has a practical significance for establishing prevention and control systems for the poisonings in the future. As an example, primers targeting the toxin genes in Lepiota have been developed for rapid PCR identification of these mushrooms by our group.
Most of the genomic statistics of L. venenata are within the ranges of those in related agaric genomes [14,17,[20][21][22][23]. One aspect does stand out: the genomic GC content of the L. venenata genome, 49%, is the highest of the three genera. In addition, there is a clear distinction between the saprotrophic and symbiotic agarics. The GC content of symbiotic Amanita is the lowest, such as A. phalloides and A. bisporigera, at 42.2 and 43.2%, respectively [20]. The genomic GC content of the saprotrophic mushrooms is higher: it is 48% for G. marginata [42], close to that of L. venenata at 49%. Our other ongoing sequencing projects with additional Amanita and Lepiota species conform to this trend. GC content often correlates in negative relationship with amount of repetitive sequences. Our ongoing sequencing of A. subjunquillea and A. pallidorosea yielded 53.26 and 55.07% of the respective genomes as repetitive sequences. For L. venenata the percentage dropped to 26.81%. The significant differences between L. venenata and the two Amanita species are consistent with the negative correlation.

Ribosomal nature of amanitin biosynthetic pathway
Lepiota venenata possesses both MSDIN and POPB genes, and the structures of these genes are consistent with those of their counterparts in both Amanita and Galerina [1,14,18]. Clearly, L. venenata uses the same ribosomal mechanism to produce α-amanitin as in Amanita [17,20] and Galerina [18]. The genes involved in the pathway, homology, and exon-intron structures, clearly indicate these pathways in the three genera are related. At this point, it may be safe to conclude that, for all the α-amanitin-producing mushrooms, the same ribosomal biosynthetic pathway based on MSDIN and POPB genes, is universally used for amanitin production.

LvPOPC and evolution of POPB genes
LvPOPC is critical for the phylogenetic reconstruction of POP history, and it was the reason that our previous report [14] was not able to fully rule out massive gene loss by comparing the gene tree with the species tree. Since phylogenetic conflict is still considered as the most powerful method to investigate HGT [43,44], we repeated the analysis with the addition of LvPOPC. As shown in Figs. 4 and 5, the POPB lineage was the only clade that conflicted with the species phylogeny (red bar in Fig. 4), and the topology conflict between the gene tree and the species tree became indisputable even when massive gene loss was considered. This result has lent strong support to the HGT hypothesis of POPB.
The reason we (and our colleagues) missed LvPOPC in the beginning was that it has low homology to the query POPA and POPB sequences. After the gene was cloned and structure determined, we realized that it is in a way similar to POPA, with similar exon-intron pattern, although shorter in length. The high divergence in the sequence may be due to its separation from other POPs long time ago.

MSDIN genes
Lepiota venenata only possesses one MSDIN sequence in the draft genome, which is consistent with the toxin profile of this mushroom (Fig. 2). The scenario is similar to that in G. marginata [18], and both of the mushrooms are saprotrophic agarics. Furthermore, two MSDIN sequences found in the genome of L. brunneoincarnata also encode α-amanitin [14]. In Amanita, subsets of MSDIN genes oftentimes reach the range of 20 to 30 [3,17,20,45]. MSDIN genes in Amanita are conserved in the leader peptide region, and they mostly start with the amino acid string of MSDIN. In contrast, the leader peptide in Galerina begins with MFDTN, and our ongoing sequencing of G. sulciceps confirms this (data not shown). In Lepiota, the second amino acid in the leader peptide region is missing, starting with M-DAN. These sequences are clearly related, and phylogenetic reconstruction of these genes showed that, on genus level they cluster by their taxonomic distribution, but within a genus they group by function, i.e., MSDIN sequences encoding same cyclic peptides group together (data not shown). This result is consistent with our recent report [14].

HGT scenarios: Transposon and gene cluster
Unlike vertical inheritance, HGT is transfer of genetic material between unicellular or multicellular organisms, rather than from parents to offspring [44,46]. HGT plays an important role in the evolution of many organisms [47,48]. Transposons are pieces of DNA that can be moved in a genome [49,50]. Horizontal transfer is an important way for transposons to avoid extinction in the host genome due to purifying selection, genetic drift or mutation inactivation [51][52][53], and it is considered to be an important driving force for genome mutation and biological transformation [54]. Transposons are divided into two categories, which can be described as copy-paste (class I) or cut-and-paste (class II) [55,56]. In L. venenata, transposons accounted for 26.81% of the whole genome and were distributed on every contig. In addition, we analyzed and counted the numbers of transposons within 50 Kb upstream and downstream of LvPOPB and LvAMA1. There are 5 transposons in the vicinity of POPB, which belong to Class I (retrotransposons). Near LvAMA1, there are 4 transposons, two of which belong to Class I, and the other two belong to Class II (DNA transposons). Although the mechanisms of HGT are not fully understood, transposon may play a role in the event, and further research is needed in the direction.
Supernumerary chromosome transfers can also be explained by interspecific mating rather than HGT [57,58], and similarly, horizontal chromosome transfer (HCT) can be the underlining cause as well [59,60]. These hypotheses can be readily tested via comparative genomic approaches, and our analyses using Symap were negative on both assumptions.
In the L. venenata genome, POPB gene is located far from the only MSDIN, LvAMA1, with at least 100 predicted genes in between. The result indicates it may be a case of single gene transfer instead of gene cluster transfer. Our other ongoing sequencing of L. brunneoincarnata returned the same result (at least 200 Kb apart, with hundreds of predicted genes in between).

Evolutionary history of the biosynthetic pathway
Collectively, the L. venenata genome offered new evidence that the key toxin-biosynthetic gene, POPB, is acquired through HGT. The genetic distances of MSDIN genes also support that MSDIN underwent the same process. POPB catalyzes the cyclization of the precursor peptide and is at the very center of the pathway. Therefore the HGT nature of this gene would likely apply to the entire pathway. Multiple lines of evidence support the HGT hypothesis. First, significant similarities are found in POPB and MSDIN gene sequences and structures across three distant families, indicating these genes are related and did not arise independently. From the phylogenetic and topology analyses, the POPB branch is robust and does not conform to the species phylogeny. The genetic distances of both POPB and MSDIN genes are considerably shorter than those of rbp2, consistent with the HGT scenario. The detailed comparison of POPB lineage to the species phylogeny illustrated that the conflict is clear, ruling out massive gene loss.

Conclusions
The genome of L. venenata illustrated that a ribosomal mechanism is the underlining mechanism of its amanitin biosynthesis. The toxin genes, LvAMA1 and LvPOPB, were far apart in the genome, not forming a gene cluster. The evidence suggested that, in the three disjunct genera, Amanita, Lepiota and Galerina, these toxin genes were acquired through HGT mechanism.

Biological materials
The mushroom fruiting bodies used in this study was L. venenata, collected in Jingzhou City, Hubei Province of central China. The fungal materials were identified by the second and the second-to-last authors who published the new species in Cai et al [25]. Specimens were deposited in the Herbarium of Cryptogams, Kunming Institute of Botany, Chinese Academy of Sciences (holotype, HKAS 101874), and Mycological Herbarium of Hunan Normal University (isotype, MHHNU 31031), respectively. The fruiting bodies were collected and wrapped in tin foil, immediately put on dry ice, and subsequently stored in − 80°C before use.

Genome sequencing and assembly
The sequencing platform for the genome of L. venenata used PacBio Sequel at NextOmics Biosciences, Wuhan, China. Sequencing and assembly were carried out using the company's standardized pipeline briefly described as the following. High quality DNA was extracted and examined for high molecular weight suitable for 20-Kb library construction. The genomic DNA was then randomly interrupted with Covaris g-TUBE. Large fragment of DNA was enriched and purified by magnetic beads. Then, the fragmented DNA was repaired. At the ends of DNA fragments, the stem circular sequencing joint was connected, and unconnected fragments removed by exonuclease. A 20-Kb library was constructed using a PacBio template prep kit and analyzed by Agilent 2100 Bioanalyzer for quality control. After the completion of the library, the DNA template and enzyme complex were transferred to the Sequel sequencer for real-time single molecule sequencing. Illumina HiSeq X10 platform was used for nucleotide level correction, based on a 350 bp library constructed, and the company's a standard method was applied.

Sequences of MSDIN and POP genes
Nucleotide sequences of MSDIN and POP genes from L. venenata genome were obtained by BLAST (NCBI BLAST + 2.4.0). Query MSDIN and POPB sequences came from A. bisporigera and G. marginata, which are well characterized by our previous molecular and biochemical approaches [1,17,18]. Comparison was also done with related sequences from A. phalloides [20]. After alignment of MSDIN and POP gene sequences, the introns and exon positions were predicted by MegAlign v7.1.0 following our previous method [14]. After the introns were removed, predicted MSDIN and POPB CDS sequences of L. venenata were obtained and named LvAMA1 and LvPOPB, respectively. One additional sequence with weak homology to both POPA and POPB was identified and denoted as LvPOPC. This gene was further cloned as described below.
Cloning and structure prediction of LvPOPC and LvAMA1 CDS sequences of LvPOPC and LvAMA1 were obtained by reverse transcription PCR (RT-PCR) using primers based on the genomic data. As BLAST result indicated, one stretch of genomic DNA sequence has homology (although weak) to both POPA and POPB. This particular DNA sequence was selected for primer design aiming to cover the full coding region, using those in POPA and POPB as the references. For the LvPOPC, the forward primer was 5′-CCCGGGTTGTAGTGGTGTAAGG-3′, and the reverse primer was 5′-ACATATTATCTCCC TGCTTTCACC-3′. For the LvAMA1, the forward primer was 5′-TCTCCAGGCCTCATTCACATTACC-3′ and the reverse primer was 5′-TGCCAGACACGGAA CAAATACATC-3′. The reactions were carried out under standard conditions. The structures of the genes were illustrated using the CDS and gDNA sequences on Splign website (https://www.ncbi.nlm.nih.gov/sutils/ splign/splign.cgi? The textpage = online&level = form).

Visualization of the genome and toxin genes
In this study, genome visualization software Circos [61] was chosen to map the genome and toxin genes. Python scripts for obtaining GC content and GC skew were generated. Genome annotation files were processed mainly through Excel, and the resulted track files were used for building information tracks in Biopython environment. The tracks were loaded into Circos to produce a genome overview in Perl environment. Coordinates of LvPOPB and LvAMA1 genes were loaded as a track to show their precise genomic locations.

Analysis of toxin by LC-MS
To verify if the genomic potential for cyclic peptide in L. venenata is consistent with the actual toxin phenotype, a liquid chromatography-mass spectrometry (LC-MS) method was applied on Aglient 6530 series system (Agilent Technologies, Palo Alto, CA, USA).
Toxins were extracted from fruiting bodies, using methanol: water: 0. 01 M hydrochloric acid (5: 4: 1) as the extraction buffer. Then, 0.06 g dried material was weighed and ground to fine powder in liquid nitrogen, 2 ml buffer were added, and the suspension transferred to 1.5 ml centrifuge tubes. The tubes sit at room temperature for 30 min, followed by centrifugation (12,000 rpm) for 3 min. Finally, the supernatant was boiled, centrifuged again, and supernatant transferred to fresh centrifuge tubes.

Phylogenetic analyses of POP genes
Coding sequences (CDS) and amino acid sequences of selected POP genes for the phylogenetic analyses were obtained from NCBI (https://www.ncbi.nlm.nih.gov/) and JGI (http://genome.jgi.doe.gov/programs/fungi/index.jsf ) (Additional file 3). Those sequences were aligned by MAFFT v7.304b [62] with default settings, and then manually adjusted with BioEdit [63] (Additional files 4 and 5). The CDS sequences and amino acid sequences of POP include 75 taxa, of which three are chosen as outgroups. For the amino acid alignment, the best model LG + G was detected by ProTest 3.4.2 [64] under Akaike Information Criterion (AIC). For the nucleotide alignment, GTR + I + G was selected as the best model for the CDSs of POP genes, using MrModeltest v2.3 [65] under AIC. Maximum Likelihood analyses and bootstrapping (1000 replicates) were performed using RAxML v7 [66]. For Bayesian inference analyses, MrBayes v3.2.6 [67] was used under the optimal substitution model calculated from MrModeltest. The posterior distributions of parameters were obtained using Markov chain Monte Carlo (MCMC) analysis for 20 million generations. Chain convergence was determined when the stopval equals to 0.01. The sampled trees were summarized after omitting the first 25% of trees as burn-in.

Topology test
In order to test whether the POPB branch in the POP phylogenetic tree is reliable, hypothetical trees were built to compare with the best tree (the CDS POP phylogenetic tree). Three different hypothetical trees were built by PAUP v4.0b10, and their site-wise log-likelihoods obtained. These hypotheses were (1) POPBs and Amanita POPA were monophyletic, (2) POPBs, Amanita POPA and Galerina POPs were monophyletic, and (3) POPBs were monophyletic with Galerina POPs. Site-wise log-likelihoods for the above hypothetical threes were generated and transferred into Consel v0.1i to perform approximately unbiased tests [68].

Gene tree vs. species tree
The POP gene and species trees were constructed from 20 taxa and 31 taxa, respectively. These taxa included representative species from Amanita, Galerina and Lepiota. For species tree, the 31 taxa (Additional file 6) contain all the 20 taxa in gene tree (Additional file 7). The rpb2 sequences [14] were obtained from our custom genomes, GenBank and JGI genomes (blastp) using A. subpallidorosea rpb2 (KP691703) as the query. Both gene tree and species tree were constructed by RAxML v7 and MrBayes v3.2.6, and the best model was GTR + I + G. Plicaturopsis crispa and Anomoporia bombycina were chosen as the outgroups. The resultant gene and species trees were analyzed by Notung 2.9 [69] with Divergence-Loss (DL) and Divergence-Transfer-Loss (DTL) models under default settings. In addition, in order to further study the relationship between the POP gene and species trees, we used codeml program in PMAL v4.9 [70] for post-phylogenetic analysis, including distance and divergence rate calculations.