Beringian origins and cryptic speciation events in the fly agaric (Amanita muscaria)

Amanita muscaria, the fly agaric, is a widely distributed ectomycorrhizal fungus known for marked morphological variation (notably pileus and wart colour) and occurrence with a broad range of host trees. Earlier molecular work identified three major geographic groups within A. muscaria that were interpreted as Eurasian, Eurasian subalpine (alpine) and North American clades, but the historical dispersal and population-structure processes that produced these groups remained unclear. Beringia (Alaska and northeastern Siberia) is biogeographically important because much of it remained ice-free during glacial maxima and because it served as a major terrestrial connection between Asia and North America; however, no A. muscaria specimens from Alaska had been included in previous phylogeographic studies. G E M L and colleagues set out to clarify the phylogenetic and phylogeographic structure of the A. muscaria species complex by sampling across Arctic, boreal and humid temperate regions of Alaska and combining newly generated sequences with homologous sequences from GenBank. They generated sequence data from three loci (protein-coding beta-tubulin, and the ribosomal ITS and LSU regions), applied multiple phylogenetic approaches (maximum parsimony, Bayesian inference, and a matrix-representation parsimony supertree), and used nested clade analysis (NCA), coalescent simulations and molecular-clock methods to infer historical demography, divergence times and the likely centre of origin for the complex.

Twenty A. muscaria specimens were collected from various regions of Alaska and deposited in the University of Alaska Fairbanks herbarium. DNA was extracted from dried sporocarps; ITS and beta-tubulin sequences from additional isolates were obtained from GenBank. Amanita pantherina sequences were used to root phylogenetic trees. PCR amplification and Sanger sequencing were performed for a portion of the beta-tubulin gene and for the ITS plus partial LSU regions using established primers and protocols. Sequence reads were edited and assembled; alignments were initiated with Clustal W and manually corrected. A hypervariable region in the beta-tubulin alignment (positions 60–86) contained a 21-bp deletion and other indels; this region was recoded to retain phylogenetic information without over-weighting deletions. Phylogenetic inference used both maximum-parsimony (MP, paup*) and Bayesian methods (MrBayes). The partition homogeneity test examined combinability of loci. Model selection for Bayesian analyses used the Akaike information criterion. MP searches used tree-bisection-reconnection (TBR) with bootstrap assessment (1000 replicates for single-locus, 100 for combined); Bayesian runs sampled trees across many generations with burn-in discarded and convergence assessed by stabilisation of likelihood scores. Kishino-Hasegawa tests compared alternative topologies. A supertree was constructed using matrix representation with parsimony (MRP). The researchers sampled 100 random Bayesian trees per locus after convergence, encoded source-tree nodes as binary characters, combined them into an MRP matrix and carried out MP analyses; node stability was assessed by bootstrap. Phylogeographic inference used nested clade analysis (NCA). Maximum-parsimony haplotype networks (tcs) defined nested clades that were analysed in geodis to compute clade distances and test genotype–geography associations via permutations (10 000 replicates). Nucleotide diversity (π) was estimated with Arlequin. Coalescent-based methods collapsed identical sequences to haplotypes and tested neutrality (Tajima’s D, Fu and Li’s D* and F*), examined recombination blocks (site compatibility matrices), and used mdiv (an MCMC coalescent method) to estimate mutation rate, divergence time, migration rate and time to most recent common ancestor (TMRCA). Genealogies, mutation ages and TMRCA were further explored with genetree. Molecular-clock analyses used ML on LSU sequences with and without clock enforcement; the Ustilaginomycetes/Hymenomycetes split was fixed at 430 Ma to estimate absolute node ages and paml was used to estimate branch lengths and standard errors.

Sequence data sets comprised ITS (717 characters; 36 parsimony-informative), beta-tubulin (468; 14 informative), LSU (625; 12 informative) and the combined alignment (1810; 62 informative). The Tamura–Nei + I model was selected as best fit for individual loci. Bayesian and MP analyses of ITS, LSU and the combined data sets consistently recovered three major clades, labelled Clades I–III; relationships among those clades were unresolved. The beta-tubulin locus provided good support for Clades I and III but did not robustly resolve Clade II’s monophyly; nevertheless, no strong topological conflicts were detected between beta-tubulin and the ribosomal loci. Each clade carried a distinctive signature in the hypervariable beta-tubulin region: Clade I showed a 21-bp deletion, Clade III had several small indels, and Clade II shared a GT motif at positions 82–83 that was absent in the other clades. In the combined data set the three clades received strong support (MP bootstrap 96%, 99% and 100% and Bayesian posterior probabilities of 1.0 for Clades I, II and III, respectively). A southeast Alaskan subclade within Clade II (II/A) was also strongly supported. Supertree (MRP) analyses—using 100 Bayesian trees per locus recoded as node characters—produced a most parsimonious tree that recovered the major clades and subclade II/A with good resolution. Tests of morphological monophyly showed that commonly recognised colour varieties (var. alba, var. formosa, var. regalis) are not monophyletic: Kishino–Hasegawa constrained-tree tests indicated that forcing each variety to be monophyletic produced significantly worse trees (all P < 0.01), consistent with ancestral polymorphism of pileus and wart colour. Phylogeographic analyses identified 25 haplotypes in Northern Hemisphere isolates; the networks for the three clades formed separate 95% connection networks but could be connected at a 92% limit. NCA results varied by clade: the total network for Clade I was consistent with contiguous range expansion (CRE). Clade II showed significant genotype–geography associations, with evidence of CRE in one nested clade (2-2) and either allopatric fragmentation (AF) or isolation by distance (IBD) in another (2-1); overall Clade II’s total cladogram supported CRE. Clade III was represented by few haplotypes and NCA could not distinguish between AF and IBD. Across the combined cladogram, NCA detected patterns consistent with ancient allopatric fragmentation explaining deep divergence among the major clades. Coalescent analyses (after removing indels) produced 18 distinct ITS haplotypes for genealogy reconstruction. The coalescent-based genealogy suggested Clade III is sister to Clade I and provided within-clade radiation age estimates expressed in coalescent units (2N): Clade I mean age 0.128, Clade II 0.276, Clade III 0.507, implying the oldest within-clade radiation in Clade III and the youngest in Clade I. The authors noted lower apparent mutation rate in Clade III. Molecular-clock ML comparisons on LSU with and without enforcing a clock did not differ significantly (χ2 test). Using the calibration constraint, the first separation within A. muscaria (between Clades I and II) was estimated at 7.48 ± 4.53 million years ago (Ma). Nucleotide diversity (π) was highest in Alaskan samples (π = 0.013094 ± 0.00702, n = 20), followed by Eurasia (π = 0.011446 ± 0.006216, n = 18) and other North American samples (π = 0.009614 ± 0.005676, n = 9). The data provided no evidence of postglacial migration of southern North American haplotypes back into Alaska; instead, Alaskan populations appear to have persisted in refugia and given rise to eastern and western North American lineages prior to the Quaternary.

Phylogenetic, supertree and coalescent results concordantly identify three distinct phylogenetic species within the Amanita muscaria complex. Rather than being strictly allopatric, representatives of all three clades occur sympatrically in Alaska; the congruent gene genealogies indicate little or no gene flow among these clades, supporting recognition as separate phylogenetic species. Multiple morphological varieties are shared across clades, and Kishino–Hasegawa tests rejected monophyly of key colour varieties; the authors interpret this pattern as ancestral polymorphism for pileus and wart colour that pre-dated the clade divergences, although environmental factors may also influence expression of colour traits. G E M L and colleagues argue that the combined phylogeographic, coalescent and nucleotide-diversity evidence points to a Beringian (Siberian–Beringian) centre of origin for the A. muscaria complex. They hypothesise that an ancestral A. muscaria population occupied humid temperate forests across Beringia in the late Tertiary, underwent fragmentation (consistent with the opening of the Bering Strait and subsequent cooling) and later diversified into the three clades. NCA and coalescent inferences indicate southward range expansions from Beringia: in North America Clade I expanded along western and eastern routes of the Rocky Mountains to give rise to western US and eastern US populations, respectively. The authors find no evidence in their sample for migration from Eurasia into Alaska; instead, Alaskan populations apparently persisted through glacial maxima in local refugia and contributed lineages to postglacial North American populations. Key uncertainties acknowledged by the investigators include incomplete geographic sampling—particularly in Asia and Siberia—which limits resolution of some phylogeographic patterns, and wide confidence around molecular-clock estimates (the LSU-based estimate for the initial split is 7.48 ± 4.53 Ma). They note limited haplotype sampling within Clade III reduced the power of NCA to distinguish among historical processes for that clade. The authors also highlight the need to verify coalescent assumptions (neutrality, absence of recombination), steps they took to test neutrality and recombination but which nevertheless constrain interpretation. Implications discussed include that Beringia may have been not only a centre of origin but a long-term refugium and biodiversity hotspot for high-latitude ectomycorrhizal fungi. The phylogeographic patterns documented for A. muscaria could be paralleled in other boreal ECM fungi. The authors recommend further sampling, especially from Siberia and other unsampled Asian regions, to refine the historical biogeography and to test the hypotheses of Beringian origin and subsequent southward expansions.