Genomic DNA was isolated from CTAB-preserved tissues or dry specimens using a CTAB isolation procedure adapted from Doyle and Doyle (1990). PCR amplification of the ITS region (nuclear ribosomal internal transcribed spacer) and LSU (large subunit ribosomal DNA) was performed using the primer pairs ITS4/ITS5 or ITS1-F/ITS4 and LR0R/LR5, respectively (http://biology.duke.edu/fungi/mycolab/primers.htm). Parts of the protein-coding genes β-tubulin, rpb2 (second largest subunit of RNA polymerase II) and tef-1 (translation elongation factor 1 alpha) were amplified using the primer pairs Am-β-tub-F/Am-β-tub-R, Am-6F/Am-7R and EF1-983F/EF1-1567R, respectively (Zhang et al. 2010). PCR products were purified by adding 1 U of Exonuclease I and 0.5 U FastAP Alkaline Phosphatase (Thermo Scientific, St. Leon-Rot, Germany) and incubating at 37 °C for 1 h, followed by inactivation at 80 °C for 15 min.

Sequencing was performed by Macrogen Inc. (Korea and The Netherlands) using the same primer combinations as for PCR, except for Am-β-tub-F, which was replaced by the shorter primer Am-β-tub-F-Seq (5’-CGGAGCRGGTAACAAYTG-3’) following Thongbai et al. (2017). The sequences were assembled in Geneious Pro v. 6.0.6 (Biomatters).

Thirty-five sequences of Amanita specimens were newly generated for this study and deposited in GenBank (http://www.ncbi.nlm.nih.gov/; Table 1). Initial BLAST searches (http://blast.ncbi.nlm.nih.gov) of both LSU and ITS-5.8S sequences were performed to estimate similarity with Amanita sequences already present in Genbank database (Table 1). Additional sequences were selected from previously published phylogenies and from GenBank (Table 1). The quality of the sequences was taken into account in selecting the sequences for the phylogenetic analyses. Materials and sequences used in this study are listed in Table 1.

A combined dataset (including nuclear ribosomal partial LSU and ITS-5.8S, partial tef1-α, rpb2 and β-tubulin genes), comprising sequences from 94 collections including the outgroup and an ITS-5.8S / LSU dataset of 69 sequences, including several clones derived from the same collections and the outgroup, were constructed and used for further phylogenetic analyses.

Amanita cf. spissacea voucher OR1214 and Amanita subjunquillea voucher HKAS63418 were used as outgroups for the combined and ITS-LSU datasets, respectively (Thongbai et al. 2017, Cui et al. 2018).

Nucleotide sequences were automatically aligned using the MUSCLE algorithm (Edgar 2004) with default settings. The alignment was further optimised and manually adjusted as necessary by direct examination with the software Se-Al v. 2.0a11 (University of Oxford).

The assignment of codon positions in the protein-coding sequences was confirmed by translating nucleotide sequences into predicted amino acid sequences using MacClade 4.0 (Maddison and Maddison 2000) and then compared with the annotated Amanita brunnescens sequences AFTOL-ID 673.

Potential ambiguously aligned segments, especially in the three introns present in tef-1 and β-tubulin gene sequences and in the ITS-5.8S alignment, were detected by Gblocks v0.91b (Castresana 2000; http://molevol.cmima.csic.es/castresana/Gblocks.html) with the following parameter settings: minimum number of sequences for a conserved position = 24 (minimum possible); minimum number of sequences for a flank position = 24 (minimum possible); maximum number of contiguous non-conserved positions = 4 bp, minimum block size = 4 bp and gaps allowed within selected blocks in half of the sequences.

To detect the possible bias from substitution saturation and to evaluate the phylogenetic signal, we tested each partition of the combined dataset and the ITS-LSU dataset by using Xia’s test (Xia et al. 2003, Xia and Lemey 2009), as implemented in DAMBE (Xia and Xie 2001). As the Iss.c is based on simulation results, there is a problem with more than 32 species. To circumvent this problem, DAMBE was used to randomly sample subsets of 4, 8, 16 and 32 OTUs multiple times and to perform the test for each subset to see if substitution saturation exists for these subsets of sequences. In order to confirm the results of the Xia’s method, we also plotted the raw number of transversions and transitions against Tamura-Nei genetic distances with the aid of the DAMBE package, with an asymptotic relationship indicating the presence of saturation.

Models of evolution for BI were estimated using the Akaike Information Criterion (AIC) as implemented in Modeltest 3.7 (Posada and Crandall 1998).

The dataset was subdivided into 10 data partitions: tef-1 1^st and -2^nd codon positions, tef-1 -3^rd codon positions, tef-1 introns and rpb2 1^st and -2^nd codon positions, rpb2 -3^rd codon positions, β-tubulin 1^st and -2^nd codon positions, β-tubulin -3^rd codon positions, β-tubulin intron, ITS, LSU. Phylogenetic analyses were performed separately for each individual and concatenated loci using Bayesian Inference (BI) as implemented in MrBayes v3. 2 (Ronquist et al. 2012) and Maximum Likelihood (ML) as implemented in RAxML 7.2.7 (Stamatakis et al. 2008).

The best-fit models for each partition were implemented as partition specific models within partitioned mixed-model analyses of the combined dataset (Table 2). All parameters were unlinked across partitions. Bayesian analyses were implemented with two independent runs, each with four simultaneous independent chains for ten million generations, starting from random trees and keeping one tree every 1000th generation. All trees sampled after convergence (average standard deviation of split frequencies < 0.01 and confirmed using Tracer v1.4 [Rambaut and Drummond 2007]) were used to reconstruct a 50% majority-rule consensus tree (BC) and to calculate Bayesian Posterior Probabilities (BPP). BPP of each node was estimated based on the frequency at which the node was resolved amongst the sampled trees with the consensus option of 50% majority-rule (Simmons et al. 2004). A probability of 0.95 was considered significant. Maximum Likelihood (ML) searches conducted with RAxML involved 1000 replicates under the GTRGAMMAI model, with all model parameters estimated by the programme. In addition, 1000 bootstrap (ML BS) replicates were run with the same GTRGAMMAI model. We provided an additional alignment partition file to force RAxML software to search for a separate evolution model for each dataset. Clades with Maximum Likelihood bootstrap values of 75% or greater were considered supported by the data.

To detect topological conflicts amongst data partitions, the nodes between the majority-rule consensus trees obtained in the ML analysis from the individual datasets were compared with the software compat.py (available at www.lutzonilab.net/downloads). Paired trees were examined for conflicts only involving nodes with ML BS > 75% (Mason-Gamer and Kellogg 1996, Lutzoni et al. 2004, Reeb et al. 2004). A conflict was assumed to be significant if two different relationships for the same set of taxa (one being monophyletic and the other not) were observed in rival trees. Sequence data and statistical analysis for each individual dataset and combined analysis are provided in Table 2.

Two major toxin-encoding genes, AMA1 and PHA1, directly encode for α-amanitin and the related bicyclic heptapeptide phallacidin, the lethal peptide toxins of poisonous mushrooms in the genus Amanita. α-Amanitin and phallacidin are synthesised as pro-proteins of 35 and 34 amino acids, respectively, in the ribosomes and are later cleaved by a prolyl oligopeptidase (Hallen et al. 2007, Luo et al. 2009, Li et al. 2014). In these pro-proteins, the amino acid sequences found in the mature toxins are flanked by conserved amino acid sequences, with an invariant Pro residue immediately upstream of the toxin regions and as the last amino acid in the toxin regions.

The toxins genes and MSDIN (cyclic peptide precursor) family members and related sequences were amplified from total genomic DNA with two consecutive PCR reactions, using the products of the first PCR as templates for the second one. For the first PCR, we used degenerated primers forward (5’ATGTCNGAYATYAAYGCNACNCG3’) and the reverse primer (5’CCAAGCCTRAYAWRGTCMACAAC3’), following the cycling condition detailed in Li et al. (2014).

For the nested PCR amplification (using the PCR products above as the amplification template of AMA1 and PHA1 genes), primers targeting conserved regions of MSDINs family were obtained from previous studies (Hallen et al. 2007, Luo et al. 2009, Li et al. 2014, Wołoszyn and Kotłowski 2017) or designed ad hoc against the conserved upstream and downstream sequences of AMA1 and PHA1 available on Genbank and tested in different combinations. For α-amanitin, we used 5’CCATCTGGGGCATCGGTTGCAACC3’ as forward primer (Li et al. 2014) in combination with the reverse primers 5’CTACGTYYGAGTCAGGACAACTGCC3’ (Li et al. 2014) and the newly generated AMA-α-R2 (5’GTCAAAGTCAGTGCGACTGCCTTGT3’) and AMA-α-R3 (5’CTGCATTTGAGTTAGGATAACGACA3’). We also tested primer pairs AMAF and AMAR 5 (Wołoszyn and Kotłowski 2017). For β-amanitin, we used forward primer AMA-β-F (5’CCATMTGGGGMATMGGTTGYRACC3’) in combination with reverse primers AMA-β-R (5’GTCMACAACTYGTATYGKCCACTACT3’), AMA-β-R2 (5’GTCMACAACTYRTATYGKCCACMGCT3’) and AMA-β-R3 (5’CCTRAYAWRGTCMACAACT3’). For PHA genes, we used forward primer 5’CCTGCYTGGCTYGTAGAYTGCCCA3’ (Li et al. 2014) in combination with the reverse primers 5’CGTCCACTACTAYDTCMARGTCAGTAC3’ (Li et al. 2014) and AMA-PHA-R2 (5’AGTCACGACTACATCGAGGTCAGTACA3’). Primer pairs FALF and FALR (Wołoszyn and Kotłowski 2017) were also tested for amplification of the phallotoxins genes.

Thermal cycling conditions were: initial denaturation at 94°C for 4 min, followed by 33 cycles of denaturation at 94°C for 30 s, annealing at 59°C for 30 s, extension at 72°C for 30 s and a final extension at 72°C for 7 min, for all reactions except for the ones involving primers from Wołoszyn and Kotłowski (2017), for which an annealing temperature of 68°C for 1 min was used.

Two groups of dried mushrooms, i.e. with cuticle (n=3) and without cuticle (n=3), were analysed. For each specimen, 100 mg of dry tissues were ground and homogenised in 3 ml extraction medium (methanol:water:0.01 M HCl [5:4:1, v/v/v]) using a tissue homogeniser. After 1 hour of incubation, all extracts were centrifuged at 5000 rpm for 5 min, the supernatant was filtered using a 0.45 mm syringe filter and 20 µl of this supernatant was injected in the RP-HPLC device for toxin detection.

The α-amanitin and phalloidin standards were obtained from Sigma-Aldrich (USA). The β-amanitin, γ-amanitin and phallacidin standards were obtained from Enzo Life Sciences (Farmingdale, NY, USA). The solvents used in this study were all HPLC grade. Stock solutions of all toxins (100 µg/ml) were prepared in methanol. The calibration standards of all toxins were diluted in the extraction fluid in concentrations of 1, 5, 20, 100, 200, 500 ng/ml. Calibration curves were produced for each toxin; they were linear over the range of interest (R² > 0.99).

Chromatography conditions for the procedure followed in this study were reported by Kaya et al. (2013, 2015). In short, the authors reported excellent separation of amatoxins and phallotoxins with the RP-HPLC and UV detection. In the laboratory, an RP-HPLC analysis of mushroom extracts was performed on a Shimadzu (Japan) HPLC system. The RP-HPLC analysis of standard solutions of α-amanitin, β-amanitin, γ-amanitin, phalloidin, phallacidin and subsequent quantification of mushroom extracts were performed on 150 × 4.6 mm, 5 mm particle, C18 column (Agilent Technologies, Palo Alto, CA) with 302 nm (for amatoxins) and 290 nm (for phallotoxins) at the UV detector. The mobile phase was used in isocratic profile with a flow rate of 1 ml/min. The content of the mobile phase was 0.05 M ammonium acetate (pH 5.5 with acetic acid)/acetonitrile (90:10 v/v). The detection limits were set at 0.6 ng/g for all toxins.

By comparing the tree topologies obtained for the individual datasets, no significant conflict, involving significantly supported nodes, was found using the 75% ML BP criterion; the datasets were therefore combined.

The test of substitution saturation (Table 2) showed that the observed index of substitution saturation (Iss) for the ITS-LSU dataset (ITS and LSU partition considered individually) the tef-1, rpb2, LSU and β-tubulin alignments of the combined dataset was significantly lower than the corresponding critical index substitution saturation (Iss.c), indicating that there was little saturation in our sequences (P < 0.001). On the other hand, the ITS partition of the combined dataset, the tef-1 introns and the β-tubulin intron showed sign of substitution saturation, indicating the unsuitability of these data for phylogenetic analysis. Nevertheless, re-analysing the ITS-LSU partition with DAMBE, after the exclusion of the 378 sites (40% of a total of 935 sites) retained by Gblocks, the substitution saturation test revealed an Iss value that was significantly (P < 0.001) lower than the Iss.c (Table 2), indicating the suitability of this data for further phylogenetic analysis. We therefore included an ITS partition, excluding the poorly aligned positions identified by Gblocks, in the combined dataset. Regarding the introns partitions, according to Thongbai et al. (2017), A. zangii and the A. ballerina clade (Fig. 1), should be considered to belong in a different section (Amanita incertae sedis), sister to the Phalloideae sensu Bas (1969). We therefore tested the combined dataset for substitution saturation by using A. zangii as the outgroup and excluding from the analysis the A. ballerina clade and the outgroup. In this case, no sign of saturation was evidenced, which supports the consistency of the phylogenetic signal in the main Phalloideae clade. We therefore decided to include the introns partitions in the phylogenetic analyses in order to increase the resolution at species level.

Figure 1. 

The 50% majority-rule consensus tree from Bayesian inference of the combined dataset. Thickened branches in bold represent ML BS support greater than 75% and BPP greater than 0.95; thickened branches in grey denote branches supported by either ML BS or BPP. For selected nodes ML BS support value and BPP are, respectively, indicated to the left and right of slashes. The new taxa are highlighted in the shaded box. AMA and PHA indicate the presence of amatoxins and phallotoxins, respectively, detected by HPLC. NT indicates not tested.

The ITS-LSU dataset and the final combined DNA sequence alignments of all loci (β-tubulin, rpb2, ITS, LSU, tef-1) alignments contained 15 and 35 OTUs and were 1575 and 3133 sites long including gaps, respectively. Sequence data and statistical analysis for each dataset are provided in Table 2.

The topologies obtained by analysing the combined dataset and the ITS-LSU dataset were highly congruent with published trees (Zhang et al. 2010, Cai et al. 2016, Thongbai et al. 2017), at least for what concerns significantly supported branches, and the Bayesian consensus trees (Figs 1 and 2) were almost identical to the optimal trees inferred under the Maximum Likelihood criterion. Several collections from tropical Africa clustered together in a well-supported clade. So far, this clade remains isolated but is notably distantly related to all other Amanita species, as yet reported from Africa (Zhang et al. 2010, Cai et al. 2016, Thongbai et al. 2017) or elsewhere and for which sequences are known (Figs 1 and 2), suggesting a common phylogenetic background. Amanita alliiodora clustered together with the two unnamed species from tropical Africa in all phylogenetic inferences considered individually or concatenated (i.e. phylogenetic species, Figs 1 and 2, shaded box).

Figure 2. 

The 50% majority-rule consensus tree from Bayesian inference of the combined nuclear ITS-5.8S and LSU sequences. Thickened branches in bold indicate bootstrap support greater than 70% and Bayesian posterior probability greater than 0.95. For selected nodes, parsimony bootstrap support value and Bayesian posterior probabilities are, respectively, indicated to the left and right of slashes. The new taxa are highlighted in the shaded box.

Morphological examination showed combinations of morphological features unique to and characteristic of each, thereby defining two morphotypes. The critical morphological features that differentiate them are the following. The first species grows under Eucalyptus. Its bulb at stipe base is (sub-)globose, neither pointed nor rooting. The ring is striated and the smell sweetish and conspicuous. The second species is not bound with Eucalyptus and has been collected in Miombo woodland and in a garden. The bulb at the stipe base is turnip-shaped to rooting. The ring is smooth or vaguely plicate and the smell weak, resembling raw potato. We therefore concluded that these two morphotypes / clades represent two distinct new species, which we describe below resp. as A. bweyeyensis sp. nov. and A. harkoneniana sp. nov.

By using a combination of the degenerated primers cited above, we obtained a complete 17-mer sequence of phallacidin precursor for the three specimens of A. bweyeyensis and the two specimens of A. harkoneniana studied (Table 3), comprising the mature toxin region sequence of phallacidin (AWLVDCP) and both the invariant Pro residues immediately preceding the mature peptide sequence and the last amino acid of the toxin. Surprisingly, this is the first time that a complete PHA sequence has been found in a species of Amanita sect. Phalloideae that does not produce this toxin. This finding is in contrast with the study of Hallen et al. (2007), concluding that all of the species synthesising amatoxins and phallotoxins, but none of the other species, hybridised to AMA and PHA genes probes (based on the same primers used in this study). However, while successful PCR amplification proves the presence of a gene (PHA gene in this case), an unsuccessful PCR, possibly due to primer mismatches, cannot be used to prove the absence of the genes encoding α- and β-amanitin, whose exact DNA sequence for these specimens is not known.