Phylogeny and biogeography of the remarkable genus Bondarzewia (Basidiomycota, Russulales)

Bondarzewia is a conspicuous and widely distributed mushroom genus, but little is known about its origin and biogeography. Here, we investigated the systematics and biogeography of Bondarzewia species using multi-locus phylogenetic analysis. Four genetic markers, including the internal transcribed spacer (ITS), large nuclear ribosomal RNA subunit (nLSU), elongation factor 1-α (tef1) and mitochondrial small subunit rDNA (mtSSU), were used to infer the phylogenetic relationships of Bondarzewia. We performed Bayesian evolutionary analysis on the gene datasets of the largest and second largest subunits of RNA polymerase II (RPB1 and RPB2). From the results, we inferred that the maximum crown age of Bondarzewia is approximately 25.5 million-years-ago (Mya) and that tropical East Asia is likely to be its ancestral area, with three possible expansions leading to its distribution in North America, Europe and Oceania.

Historical biogeography of Bondarzewia. The inferred historical biogeographic scenarios from the analyses conducted using LAGRANGE 17 and RASP 18 are shown in Fig. 3. Our biogeographical analyses indicated that, of the 11 phylogenetic species that were identified in this lineage, 4 (37%), 3 (27%), 2 (18%), 1 (9%), and 1 (9%) of the species were reported to be from East Asia, Oceania, North America, Europe and South America, respectively. The Bayesian binary Markov chain Monte Carlo analysis showed the East Asia presented the highest probability (78%) of being the ancestral area of Bondarzewia. The maximum likelihood-based estimation also provided strong support for East Asia being the ancestral area. In addition, the basal species (B. guaitecasensis, B. kirkii, B. podocarpi, B. propria, and B. retipora) exhibited a pantropical distribution pattern (Figs 1 and 3). These findings support that Bondarzewia originated in tropical East Asia. Meanwhile, the East Asian and North American ancestral origins of the Holarctic Bondarzewia species are also supported by LAGRANGE. The three clades indicated three kinds of intercontinental distribution patterns respectively: East Asia-Oceania-South America, East Asia-North America, and East Asia-Europe (Fig. 4).

Discussion
In this study, we presented results from investigations on the phylogeny, origin and biogeography of the genus Bondarzewia. Three major clades and one isolate species comprising 11 species were identified within the genus (Fig. 1). In the following discussion, we focused on the features of the major clades and their distribution patterns.
Clade I is weakly supported and includes five taxa representing two sister groups (Fig. 1). B. podocarpi and B. propria formed a significantly supported gymnosperm-associated group (100% MP, 100% BS, 1.00 BPP); they share dimitic hyphal structure, have similarly sized pores (2 per mm) and cyanophilous basidiospores. To date, B. podocarpi has been found only in southern tropical China on Podocarpus and Dacrydium, which are native flora of the Southern Hemisphere 6,7 . B. propria was first found in New Zealand and was associated with the endemic plants Dacrydium cupressinum and Agathis australis 7 . B. retipora, B. kirkii and B. guaitecasensis formed an angiosperm-associated species group (91% MP). These three species grow on native evergreen angiosperm trees of the Nothofagaceae, Meliaceae or Lamiaceae families, and morphologically, they share orange pileal surfaces and a dimitic hyphal structure 7,11 . The data from the current work indicated close close relationships among these species. B. dickinsii and B. occidentale were grouped together with low support (50% MP) in our phylogeny. They both grow in temperate areas and have white pore surfaces and cyanophilous basidiospores 7 . B. dickinsii is associated with Fagaceae such as Quercus and Castanea, and has been found from East Asia; however, B. occidentale grows on gymnosperm trees such as Picea and Tsuga, and has been found only from West America so far 7 .
The East Asian species B. submesenterica and B. tibetica and the European species B. mesenterica were grouped together ( Fig. 1). They all produce brown pileal surfaces, cream pore surfaces, and grow mainly on conifers. B. submesenterica and B. tibetica were found in temperate areas of the Hengduan-Himalayan region, which is a global biodiversity hotspot located in China. B. mesenterica is common in Central-East and South Europe, which was also a refuge for organisms 19 . The host plants such as Abies, Picea and Pinus also have this kind of distribution pattern 20,21 . These results indicated that the appropriate growth environments of refuges may have guaranteed the survival of Bondarzewia during the Quaternary Ice Age.
Bondarzewia berkeleyi fromed a lineage in the multi-loci phylogeny (Fig. 1). This species has been delimited as an East American species and grows mainly on deciduous plants such as Quercus and Acer 7,10 . Phylogenetically, B. berkeleyi diverged the earliest of the Holarctic Bondarzewia species and occupies a basal and separate position.
The maximum crown age of Bondarzewia was estimated to be around the Late Oligocene (25.54 ± 0.17 Mya). During that period, the general landforms of the modern world had already formed, and dramatic climate changes was ongoing 22,23 . The resultant polar ice sheets advanced and retreated several times and covered most of North America and northern Europe until the Last Glacial Maximum (ca. 18,000 years ago) 23   Our results supported that Bondarzewia originated in tropical East Asia. However, the basal species were associated with Nothofagaceae, Meliaceae, Lamiaceae, Podocarpaceae and Araucariaceae; according to the co-evolution of fungi and host plants, the Gondwana origin cannot be rejected because Araucariaceae, Podocarpaceae, Nothofagaceae and Meliaceae are thought to have originated in Gondwana [30][31][32][33][34] . The Holarctic species including Clade II, Clade III and B. berkeleyi diverged later (8.5 Mya ) and were restricted to temperate zones and temperate plants such as Picea, Abies, Quercus and Castanea 7,8 . The temperate habit and late divergence time suggest that an adaptation to temperate climates occurred.
Species in the Clade I appeared the earliest in our phylogenies and were distributed in South China, New Zealand, Australia and South America (Figs 1 and 3). Species could have migrated within tropical South East Aisa, New Zealand and Australia (Fig. 4), when the Oceania and Asian plates were connected after their collision close to the Oligo-Miocene boundary 35,36 . The evolutionary studies of organism on either side of Wallace's Line could clearly confirm the migration between Oceania and Asian plates 37 . Interestingly, our data suggest that B. guaitecasensis, B. kirkii and B. retipora are closely related and diverged comparatively recently (3.6 Mya). As we all known, there are twenty thousand islands in South Pacific such as Solomon island archipelago, New Caledonia and Hawaii islands. The scaly tree ferns has been proved to dispersal between Oceania and South America by wind via the lightweight spores 38 . We speculate that the wind and ocean current drive the dispersal of Bondarzewia in South Pacific island by island.
Two sister species within Clade II, B. dickinsii and B. occidentale, exhibited an East Asia-Western North America temperate disjunct distribution (Figs 1 and 3), and their divergence time was estimated at approximately 5.9 Mya, implying that dispersals between East Asia and Western North American occurred, probably via the Bering Land Bridge (BLB) prior to their divergence (Fig. 4). The BLB separated from 5.5-5.4 Mya 39,40 , and connected the East Asia and North America before that time. The BLB has been confirmed as a route in the dispersal event of many other organisms 14,15,[41][42][43] . The divergence of B. dickinsii and B. occidentale may result from the separate of BLB and the extinction of species occurred. At the same time, the continuous significant climate changes between 15 Mya and the Last Glacial Maximum may have played important role in breaking the biotic connections between Tertiary floras of East Asia and North America and in facilitating allopatric speciation 22,23,44 .
An East Asian-European allopatric speciation was also inferred for B. submesenterica, B. tibetica and B. mesenterica (Figs 3 and 4). The estimated divergence occurred at approximately 3.1 Mya, and East Asia was inferred to be the most likely ancestral area. Long-distance dispersal between Europe and East Asia may have occurred, which was common in the immigration of plants such as Quercus, Salix, Populus, Picea, Abies and Larix 20,21,25,28,29,45 . We speculate that the severe climate changes that occurred since 15 Mya and the sebsequent aridification in Central Eurasia 22 may have been the causes of the divergence.
The East American species, B. berkeleyi, occupies the basal and separate position in the phylogenetic tree of the Holarctic Bondarzewia species (Figs 1 and 3). The dating analysis inferred an earlier divergence time (8.5 Mya). Its hosts such as Quercus and Acer have an extensive distribution in the Northern Hemisphere. We speculate its ambiguous affinities with other species may relate to an incomplete sampling and undescribed species in Central American. Alternatively, B. berkeleyi could be an earlier relic of Bondarzewia in the Western North America and the extinction of this species in the Holarctic region occurred after it originated from tropical East Asia. The separate position of East American B. berkeleyi deserves a detailed study.

Conclusion
The monophyletic genus Bondarzewia originated in the tropical zone of East Asia and has diverged since 25.54 ± 0.17 Mya (15.65-37.18 Mya, 95% HPD). The severe climate changes and the resulting reduction in sea-levels and aridification of the center of continents shaped the evolutionary history of Bondarzewia via the co-diversification of the fungi and their host plants. Three intercontinential distribution patterns were recognized: East Asia-Oceania-South America, East Asia-North America, and East Asia-Europe. More samples and sequences are needed to improve our understanding of the biodiversity and biogeography of Bondarzewia. The genus diversity still deserves a further study because the samples are incomplete in many places such as Central Asia, Siberia, Indonesia and South Africa, and the evolution history of some species are still ambiguous.

Materials and Methods
Taxa sampling. The sampled taxa, their genetic markers, and their GenBank accession numbers are provided ( Table 2). Specimens from East Asia (EA), Europe (EUR), North America (NA), Oceania (OC), and South America (SA) were studied for their morphological characters. The outgroup taxa were determined based on the previous phylogenetic study that focused on the diversity of Bondarzewia 8 . DNA extraction, PCR, and DNA sequencing. A rapid plant genome DNA extraction kit (Aidlab Biotechnologies Co., Ltd, Beijing, China) and the primers listed in Table 3 were used to obtain PCR products from dried specimens and cultures according to the manufacturer's instructions with modifications. The PCR protocols for ITS, nLSU, mtSSU, tef1, RPB1 and RPB2 have been described previously publication by the same research group 46 . The PCR products were purified and sequenced at the Beijing Genomics Institute (China) using the same primers. All newly generated sequences were deposited in GenBank (http://www.ncbi.nlm.nih.gov/).

Sequence alignments and phylogenetic analyses. The sequences of Heterobasidion annosum (Fr.)
Bref. and Heterobasidion parviporum Niemelä & Korhonen were used as outgroups 8 . Phylogenetic analyses was applied to single-locus genealogies for ITS, nLSU, mtSSU and tef1, and concatenated dataset that contained the ITS + nrLSU + mtSSU + tef1 sequences. Initially, the four genes were aligned using MAFFT 6 (http://mafft.cbrc. jp/alignment/server/) 47 with "G-INS-I" strategy and then the alignment was manually optimized in BioEdit 48 . Finally, the four gene fragments were concatenated with SEAVIEW 4 49 for further phylogenetic analysis. One thousand partition homogeneity test (PHT) replicates of ITS, nrLSU, mtSSU, and tef1 sequences were tested by PAUP* version 4.0b10 (Swofford, 2002) to determine whether the partitions were homogeneous. The PHT results indicated all the DNA sequences display a congruent phylogenetic signal (P value = 0.02).
ML, MP and BI methods were used to analyze the compiled datasets. A suitable substitution model for each partition of the dataset was determined using the Akaike Information Criterion implemented in MrMODELTEST2.3 50 . PAUP* 4.0b10 51 was used for MP analysis. All characters were equally weighted, and gaps were treated as missing data. Trees were inferred using the heuristic search option with TBR branch swapping and 1000 random sequence additions. Max-trees was set to 5000, branches of zero length were collapsed, and all parsimonious trees were saved. Clade robustness was assessed using bootstrap analysis with 1000 replicates. Descriptive tree statistics including the tree length (TL), consistency index (CI), retention index (RI), rescaled consistency index (RCI), and homoplasy index (HI), were calculated for each maximum parsimony tree generated.
ML searches conducted with RAxML-HPC2 52 on Abe through the Cipres Science Gateway (www.phylo.org) involved 100 ML searches under the GTR + GAMMA model; all model parameters were estimated by the program. Only the best maximum likelihood tree from all searches was kept. In addition 100 rapid bootstrap replicates were run with the GTR + CAT model to assess the reliability of the nodes.
BI was examined with MrBayes3.1.2 53 with a general time-reversible model of DNA substitution and an inverse-gamma distribution rate variation across sites. Four Markov chains were run from the random starting tree for 1 million generations for the combined datasets. Trees were sampled every 100 generations. The burn-in was set to discard the first 25% of the trees. A majority rule consensus tree of all the remaining trees was used to calculate BPP.
Branches that received bootstrap support for MP, BS and BPP greater than or equal to 75% (MP/BS) and 0.95 (BPP) were considered to be significantly supported.

Divergence time estimation.
Several studies have attempted to date the evolutionary splits of fungi using various calibration strategies 14,15,41 . Here, we used internal calibration to determine the divergence time between Ascomycota and Basidiomycota, 582 Mya, with the 400-million-years-old fossil Paleopyrenomycites devonicus Taylor, Hass, Kerp, M. Krings & Hanlin 16 . A normal distribution was applied by setting the mean and the standard deviation to 582.5 and 50.15, respectively 16 . One constraint was applied: the initial diversification of the Russulales was set at 189 Mya, consistent with the conservative estimated divergence time for the plant family Pinaceae 54 . Co-divergence between fungal lineages and their plant hosts suggest that Russulales and its allies occurred around, or slightly later than, the time of the diversification of Pinaceae 21,55 .
The BEAST 1.8.0 software package was used to estimate divergence times 56 . The two gene fragment, RPB1 and RPB2, were concatenated for molecular dating. We retrieved the sequences of six additional species-Mar-  First, we used BEAUti to generate xml files that were executable in BEAST. The RPB1 and RPB2 datasets were set as two partitions, the substitution and molecular clock models were set as unlinked, and the inferred trees were set as linked. For both partitions, the GTR model was chosen as the best substitution model by MrModelTest, and a relaxed lognormal model was employed for molecular clock analysis 59 . The tree prior was set to Yule speciation. For each analysis, two independent runs were conducted for 100 million generations. Log files of the two runs were combined using LogCombiner by setting the first 10% of the logs as burn-ins and then analyzed in Tracer 1.5 (http://tree.bio.ed.ac.uk/software/figtree/tracer). The resulting trees were also combined, interpreted in TreeAnnotator, and viewed in FigTree 1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/).
We also estimated the divergence time of the main nodes in Bondarzewia using a mini ITS dataset containing representatives of all 11 species. The estimated crown age of the genus Bondarzewia inferred by the RPB1 and RPB2 data was used to calibrate the ITS phylogeny by setting the prior to a normal distribution. The other procedures were the same as the ones applied in the estimation using the RPB1 and RPB2 datasets. The most recent common ancestors were only defined for the major clades that were well-supported in the ITS + nLSU + mtSSU + TEA phylogenies.
Inferring the geographic center of of Bondarzewia's origin. The phylogeny and divergence inferred from the ITS dataset were used to reconstruct the possible historical distributions of Bondarzewia lineages. Both the maximum likelihood-based estimations implemented in LAGRANGE and the Bayesian binary Markov chain Monte Carlo analysis provided by RASP v3.2 were used. The geographic distributions of Bondarzewia lineages were classified into five areas: East Asia, Europe, North America, South America and Oceania. The probabilities of dispersal were estimated according to the divergence times inferred earlier in this study between different areas as previously summarized 60 . Bayesian binary analysis was conducted in RASP by setting the generations to 10 million and by discarding the first 10% of samples as burn-ins; the other parameters used were the default settings. ArcGIS v10.1 (http://esri.com/arcgis) was used to visualize the geographic distribution and possible dispersal routes of Bondarzewia.   Table 3. PCR primers used in this study. a Degeneracr codes: S = G or C, W = A or T, R = A or G, Y = C or T, N = A or T or C or G, D = G or A or T, M = A or C. *ITS, internal transcribed spacer region; nLSU, the large nuclear ribosomal RNA subunit; RPB1, the largest subunit of RNA polymerase II; RPB2, the second subunit of RNA polymerase II.