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.

Bondarzewia Singer (Bondarzewiaceae, Russulales) is a globally distributed genus of mushroom forming fungi. Some species are edible and have medicinal potential12, whereas some are considered to be forest pathogens3. Bondarzewia can be mistaken for the mycorrhizal genus Lactarius4. Phylogenetically, Bondarzewia forms sister relationship with the genus Heterobasidion in Bondarzewiaceae, but Lactarius is closed to Russula in Russulaceae5. Species of Bondarzewia are not mycorrhizal5, and eleven species are currently accepted in the genus: B. dickinsii (Berk.) Jia J. Chen, B.K. Cui & Y.C. Dai, B. podocarpi Y.C. Dai & B.K. Cui, B. submesenterica Jia J. Chen, B.K. Cui & Y.C. Dai and B. tibetica B.K. Cui, J. Song & Jia J. Chen (reported from Asia)678, B. mesenterica (Schaeff.) Kreisel (from Europe)9, B. berkeleyi (Fr.) Bondartsev & Singer and B. occidentalis Jia J. Chen, B.K. Cui & Y.C. Dai (from North America)710, B. kirkii J.A. Cooper, Jia J. Chen & B.K. Cui, B. propria (Lloyd) J.A. Cooper and B. retipora (Cooke) M.D. Barrett (from Oceania), and B. guaitecasensis (Henn.) J.E. Wright from (South America)711.

Previous studies of Bondarzewia were mainly based on morphological characters, and three species from North and South America were confirmed101112. Recent, phylogenetic analyses of Bondarzewia have been based on the sequences of the internal transcribed spacer (ITS) and the large nuclear ribosomal RNA subunit (nLSU)678. Five new species and three new combinations were established with clear interspecific affinities: B. dickinsii and B. occidentalis group together and show a close relationship; B. podocarpi, B. propria, B. kirkii and B. retipora cluste together and from two sister groups with different hosts; B. submesenterica and B. tibetica are closely related with B. mesenterica.

Bondarzewia is a marcofungal genus with a few species and comparatively complete records worldwide, and the interspecific affinities are clear. However, a comprehensive estimation of divergence time is lacking, and the biogeography of the Bondarzewia mushrooms is not well understood. Molecular phylogeny has been widely used to delineate the lineages and their biogeographic distribution, and major ecological and geological events can be dated more accurately by applying the molecular clock, with the proper calibrations determined by fossils to gene phylogenies13. Fungal fossil records are rare and not frequent in the evolutionary history of fungi. Recent molecular studies on basidiomycetous fungi determined the divergence between Basidiomycota and Ascomycota (i.e., 582 Mya) based on a 400-million-year-old fossil of Paleopyrenomycites devonicus13141516.

The aims of the present study were to assess the divergence and biogeography of Bondarzewia. Therefore, we carried out multi-locus phylogenetic analyses with five ribosomal DNA genes and one mitochondrial gene, i.e., ITS, nLSU, elongation factor 1-α (tef1), the largest and second largest subunits of RNA polymerase II (RPB1 and RPB2, respectively), and mitochondrial small subunit rDNA (mtSSU). These genes are common in fungi identification which could delimitate the species very well.

Results

Phylogenetic analysis with the combined dataset

The ITS dataset included sequences from 27 fungal samples representing 13 taxa. The dataset had an aligned length of 576 characters, of which 418 are constant, 17 are variable and parsimony uninformative, and 141 are parsimony informative. The nLSU dataset included sequences from 26 fungal samples representing 12 taxa. The dataset had an aligned length of 881 characters, of which 819 are constant, 2 are variable and parsimony uninformative, and 60 are parsimony informative. The mtSSU dataset included sequences from 22 fungal samples representing 11 taxa. The dataset had an aligned length of 485 characters, of which 389 are constant, 45 are variable and parsimony uninformative, and 53 are parsimony informative. The tef1 dataset included sequences from 23 fungal samples representing 11 taxa. The dataset had an aligned length of 540 characters, of which 389 are constant, 8 are variable and parsimony uninformative, and 143 are parsimony informative. The best fitting model identified for the single gene datasets (ITS, nLSU, mtSSU and tef1) were the same: General Time Reversible + Proportion Invariant + Gamma (GTR + I + G). Maximum parsimony (MP), Maximum likehood (ML), and Bayesian inference (BI) analyses yielded similar tree topologies for each gene, and only only the support values in the nodes are different. Only the ML tree with maximum parsimony (MP), maximum likelihood (BS) and Bayesian posterior probabilities (BPP) value were provided (see Supplementary Fig. S1). The phylogenetic results based on mtSSU did not clearly differentiate taxa of Bondarzewia, and revealed weakly supported clades or subclades (Fig. S1). ITS, nLSU or tef1 dataset could clearly differentiate the taxa of Bondarzewia (Fig. S1).

The ITS + nLSU + mtSSU + tef1 sequence matrix contained 28 taxa and 2483 aligned base pairs (bp), of which 576 bp, 881 bp, 485 bp and 540 bp from ITS, nLSU, mtSSU and tef1, respectively. The best fitting model identified for the combined dataset were (GTR + I + G). The MP analysis yielded 4 equally parsimonious trees (TL = 661, CI = 0.834, RI = 0.917, RC = 0.764, HI = 0.166), and only the position of species in each clades behave little different. ML analysis and BI yielded similar tree topologies to the one inferred by MP analysis and only the support values in the nodes are different (Fig. 1). The taxa of Bondarzewia could differentiate clearly and the backbones of the phylogenetic tree were highly supported.

Figure 1

Phylogenetic tree within Bondarzewia inferred from the Maximum Parsimony (MP) analysis based on the ITS + nLSU + mtSSU + tef1 dataset.

Branches are labeled for MP/BS and BPP values greater than 50% and 0.95, respectively.

Based on the combined dataset analyses, our results showed that the genus Bondarzewia forms a group with strongly support (100% MP, 100% BS, 1.00 BPP) and can be divided into three distinctive clades along with B. berkeleyi (Fig. 1). Clade I is moderately supported by MP analyses (55% MP). This clade comprises five species from East Asia, Oceania and South America. Clade II is composed of two species and is moderately supported by MP analyses (50% MP). Within this clade, B. dickinsii is distributed in East Asia and B. occidentale covers the East Coast of America. Clade III is composed of three species from Eurasia and is well-supported (99% MP, 100% BS, 1.00 BPP).

Divergence of Bondarzewia lineages

The alignment of the two datasets (RPB1 and RPB2), which are 1554 and 1462 bp in length respectively, consisted of 39 taxa. Analyses calibrated by P. devonicus (Fig. 2), 582 Mya between Ascomycota and Basidiomycota, estimate the divergence time of Russulales at 173.78 ± 0.47 Mya (127.42–220.43 Mya, 95% HPD) that meet the constraint of Russulales. The initial diversification of Bondarzewia is at Late-Oligocene, 25.54 ± 0.17 Mya (15.65–37.18 Mya, 95% HPD). The estimated divergence time of two fossil records point for comparison, Archaeomarasmius leggetti Hibbett, D. Grimaldi & Donoghue and Quatsinoporites cranhamii S.Y. Sm., Currah & Stockey, are 159.45 ± 0.45 Mya (118.75–208.64 Mya, 95% HPD) and 143.34 ± 0.42 Mya (93.13–195.41 Mya, 95% HPD). The estimated divergence times for other nodes are summarized in Table 1.

Figure 2

Chronogram and estimated divergence times of Bondarzewia generated by molecular clock analysis using the RPB1 and RPB2 dataset.

The chronogram was obtained using the Ascomycota–Basidiomycota divergence time of 582 Mya as the calibration point. The calibration point and objects of this study are marked in the chronogram. The geological time scale is in millions of years ago (Mya).

Table 1

Estimated divergence times of each time.

Node	Mean ± standard error	95% HPD	Node	Mean ± standard error	95% HPD
A: Ascomycota/Basidiomycota	564.76 ± 0.34	467.45–666.26	H: Boletales/Agaricomycetes	207.67 ± 0.60	154.78–260.51
B: Ustilaginomycotina/Basidiomycota	429.61 ± 0.53	335.21–529.36	I: Agaricales/Russulales	205.64 ± 0.59	153.73–258.69
C: Pucciniomycotina/Agaricomycotina	404.45 ± 0.53	312.45–497.78	J: Russulales	173.78 ± 0.47	127.42–220.43
D: Tremellomycetes/Agaricomycotina	382.31 ± 0.53	293.36–473.40	K: Bondarzewia	25.54 ± 0.17	15.65–37.18
E: Dacrymycetes/Agaricomycetes	346.66 ± 0.53	261.10–430.92	*C1: Hymenochaetaceae	143.34 ± 0.42	93.13–195.41
F: Gomphales/Agaricomycetes	261.17 ± 0.58	196.72–329.66	*C2: mushroom-forming fungi	159.45 ± 0.45	118.75–208.64
G: Hymenochaetales/Agaricomycetes	247.18 ± 0.57	184.07–310.49

Fossil record (for comparison): *C1: Quatsinoporites cranhamii, 125–130 Mya ; *C2: Archaeomarasmius leggetti, 90 Mya (minimum age).

Historical biogeography of Bondarzewia

The inferred historical biogeographic scenarios from the analyses conducted using LAGRANGE17 and RASP18 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).

Figure 3

Divergence time estimation and ancestral area reconstruction of Bondarzewia using the ITS dataset.

The chronogram was obtained by molecular clock analysis using BEAST. The pie chart in each node indicates the possible ancestral distributions inferred from Bayesian Binary MCMC analysis implemented in RASP. The characters above and below each branch identify the possible ancestral distribution estimated by maximum likelihood-based program LAGRANGE. The color key lists the possible ancestral ranges at different nodes; black with an asterisk represents other ancestral ranges.

Figure 4

Map of the geographic distribution of Bondarzewia and its possible dispersal routes generated by ArcGIS v10.1 (http://esri.com/arcgis).

A hypothetical schematic depicting of the place of origin, their migration routes, and the rapid radiation and speciation of Bondarzewia.

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 Hemisphere67. B. propria was first found in New Zealand and was associated with the endemic plants Dacrydium cupressinum and Agathis australis7. 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 structure711. 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 basidiospores7. 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 far7.

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 organisms19. The host plants such as Abies, Picea and Pinus also have this kind of distribution pattern2021. 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 Acer710. 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 ongoing2223. 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. Meanwhile, belts of relatively arid climates in the centers of Eurasia and North America formed24. Previous study on organisms, including the host plants of Bondarzewia such as Abies, Salix and Populus, proved that dispersal, extinction and speciation occurred during this time because of the dramatic climate changes202526272829.

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 Gondwana3031323334. 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 Castanea78. 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 boundary3536. The evolutionary studies of organism on either side of Wallace’s Line could clearly confirm the migration between Oceania and Asian plates37. 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 spores38. 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 Mya3940, 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 organisms1415414243. 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 speciation222344.

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 Larix202125282945. We speculate that the severe climate changes that occurred since 15 Mya and the sebsequent aridification in Central Eurasia22 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 Bondarzewia8.

Table 2

Information of sequences used in this study.

Species	Sample no.	GenBank accessions
		ITS	nLSU	RPB1	RPB2	tef1	mtSSU
Albatrellus higanensis	AFTOL-ID 774	—	—	AY788846	AY780935	—	—
Aleurodiscus wakefieldiae	He 2580	—	—	KX577720a	KX577723a	—	—
Agaricostilbum hyphaenes	AFTOL-ID 675	—	—	AY788845	AY780933	—	—
Boletus edulis	HMJAU 4637	—	—	KF112586	KF112704	—	—
Bondarzewia berkeleyi	Dai 12759	KJ583202	KJ583216	KX066152a	KX066162a	KX066138a	KX066169a
B. berkeleyi	Dai 16052	KX263720a	KX263722a	—	—	—	—
B. dickinsii	Cui 8682	KJ583209	KJ583223	KX066150a	KX066160a	KX066136a	KX066167a
B. dickinsii	Dai 13413	KJ583210	KJ583224	KX066151a	KX066161a	KX066137a	KX066168a
B. dickinsii	Li 150909/19	KX263721a	KX263723a	—	—	—	—
B. dickinsii	Li 1097	FJ644288	—	—	—	—	—
B. guaitecasensis	Rajchenberg 11898	FJ644287	—	—	—	—	KX066175a
B. kirkii	PDD 94520	KJ583215	KJ583229	—	—	KX252748	KX066180a
B. kirkii	JAC 10839	KJ734674	KM067469	—	—	KX252747	KX066179a
B. kirkii	PL 450211	KJ583214	KJ583228a	—	—	KX252746a	KX066178a
B. mesenterica	DD 348/06	KM243328	KM243331	—	—	KX066147a	KX066182a
B. mesenterica	Niemelä 5374	KM067468	KM067470	—	—	KX066146a	KX066181a
B. occidentalis	HHB 14803	KM243329	KM243332	KX066156a	KX066163a	KX066142a	KX066176a
B. occidentalis	Lowe 7887	KM243330	KM243333	KX066157a	KX066164a	KX066143a	KX066177a
B. occidentalis	AFTOL-ID 452	DQ200923	DQ234539	DQ256049	AY218474	DQ059044
B. podocarpi	Cui 6380	KJ583206	KJ583220	KX577718a	KX577721a	KX252745aa	KX066174a
B. podocarpi	Dai 9261	KJ583207	KJ583221			KX252743a	KX066172a
B. podocarpi	Dai 12038	KJ583208	KJ583222			KX252744a	KX066173a
B. propria	PDD 60293	KJ583213	KJ583227
B. retipora	LNP 68	KJ747633	KJ747630			KX066144a
B. retipora	LNP 75	KJ747632	KJ747629			KX066145a
B. submesenterica	Cui 9809	KJ583203	KJ583217			KX066141a	KX066171a
B. submesenterica	Cui 10345	KJ583204	KJ583218			KX066140a	KX066170a
B. submesenterica	Cui 10724	KJ583205	KJ583219	KJ651627	KJ651720	KX066139a
B. tibetica	Cui 12078	KT693202	KT693204	KX066159a	KX066166a	KX066149aa	KX066184a
B. tibetica	Yu 56	KT693203	KT693205	KX066158a	KX066165a	KX066148a	KX066183a
Calocera cornea	AFTOL-ID 438	—	—	AY857980	AY536286	—	—
Clavaria zollingeri	AFTOL-ID 563	—	—	AY857987	AY780940	—	—
Coltricia perennis	AFTOL-ID 447	—	—	AY864867	AY218526	—	—
Coniophora puteana	MUCL 1000	—	—	GU187451	GU187778	—	—
Coprinopsis cinerea	okayama 7#130	—	—	XM_001828525	XM_001829088	—	—
Cryptococcus humicola	AFTOL-ID 1552	—	—	DQ645518	DQ645517	—	—
Dacryopinax spathularia	AFTOL-ID 454	—	—	AY857981	AY786054	—	—
Echinodontium tinctorium	AFTOL-ID 455	—	—	AY864882	AY218482	—	—
Fomitiporia hartigii	Dai 11766	—	—	KJ651628	KJ651721	—	—
F. mediterranea	AFTOL-ID 688	—	—	AY864869	AY803748	—	—
Gautieria otthii	AFTOL-ID 466	—	—	AY864865	AY218486	—	—
Heterobasidion annosum	06129/6	KJ583211	KJ583225	KF006499	KF033133	KX252741a	KJ651577
H. australe	04164/3	—	—	KF033134	KF006500	—	—
H. parviporum	H 091605	KJ651503	KJ651561	KJ651657	KJ651750	KX252742a	KJ651622
Lactarius deceptivus	AFTOL-ID 682	—	—	AY864883	AY803749	—	—
Malassezia globosa	CBS 7966			KF706493	KF706518
Marasmius rotula	AFTOL-ID 1505	—	—	DQ447922	DQ474118	—	—
Mycena amabilissima	AFTOL-ID 1686	—	—	DQ447926	DQ474121	—	—
M. aurantiidisca	AFTOL-ID 1685	—	—	DQ447927	DQ474122	—	—
Neurospora crassa	OR74A	—	—	XM959004	AF107789	—	—
Pseudowrightoporia japonica	Dai 12086			KX577719a	KX577722a
Ramaria rubella	AFTOL-ID 724	—	—	AFTOL database	AY786064	—	—
Rhizopus stolonifer	AFTOL-ID 632	—	—	AFTOL database	AFTOL database	—	—
Russula sp.	H-2009BT109A	—	—	JN389202	JN389208	—	—
Schizophyllum commune	TUB 012156			DQ068011	KC904262	—	—
Schizosaccharomyces pombe	972h-	—	—	NM001021568	NM001018498
Serpula lacrymans	REG 383	—	—	GU187485	GU187809	—	—
Stereum hirsutum	AFTOL-ID 492	—	—	AY864885	AY218520	—	—
Suillus pictus	AFTOL-ID 717			AY858965	AY786066
Ustilago maydis	AFTOL-ID 505	—	—	AFTOL database	AY485636	—	—

aNewly generated sequences.

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 group46. 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/).

Table 3

PCR primers used in this study.

Gene*	Primer	Primer sequences (5′-3′)a	Reference
ITS	ITS5	GGA AGT AAA AGT CGT AAC AAG G	White et al. (1990)
	ITS4	TCC TCC GCT TAT TGA TAT GC	White et al. (1990)
nLSU	LR0R	ACC CGC TGA ACT TAA GC	http://www.biology.duke.edu/fungi/mycolab/primers.htm
	LR7	TAC TAC CAC CAA GAT CT	http://www.biology.duke.edu/fungi/mycolab/primers.htm
RPB1	RPB1-Af	GAR TGY CCD GGD CAY TTY GG	Matheny et al. (2002)
	RPB1-Cf	CCN GCD ATN TCR TTR TCC ATR TA	Matheny et al. (2002)
RPB2	fRPB2-5F	GAY GAY MGW GAT CAY TTY GG	Liu et al. (1999); Matheny (2005)
	fRPB2-7cR	CCC ATR GCT TGY TTR CCC AT	Liu et al. (1999); Matheny (2005)
mtSSU	MS1	CAG CAG TCA AGA ATA TTA GTC AAT G	White et al. (1990)
	MS2	GCG GAT TAT CGA ATT AAA TAA C	White et al. (1990)
tef1	983F	GCY CCY GGH CAY CGT GAY TTY AT	http://ocid.NACSE.ORG/research/deephyphae/EF1primer.pdf
	1567R	ACH GTR CCR ATA CCA CCR ATC TT	http://ocid.NACSE.ORG/research/deephyphae/EF1primer.pdf

aDegeneracr 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.

Sequence alignments and phylogenetic analyses

The sequences of Heterobasidion annosum (Fr.) Bref. and Heterobasidion parviporum Niemelä & Korhonen were used as outgroups8. 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 BioEdit48. Finally, the four gene fragments were concatenated with SEAVIEW 449 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.350. PAUP* 4.0b1051 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-HPC252 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.253 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 strategies141541. 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 & Hanlin16. A normal distribution was applied by setting the mean and the standard deviation to 582.5 and 50.15, respectively16. 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 Pinaceae54. 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 Pinaceae2155.

The BEAST 1.8.0 software package was used to estimate divergence times56. The two gene fragment, RPB1 and RPB2, were concatenated for molecular dating. We retrieved the sequences of six additional species—Marasmius rotula (Scop.) Fr., Mycena amabilissima Peck, M. aurantiidisca (Murrill) Murrill, Fomitiporia hartigii (Allesch. & Schnabl) Fiasson & Niemelä, F. mediterranea M. Fisch., and Coltricia perennis (L.) Murrill—as representative taxa of the initial diversification of mushroom-forming fungi (based on the 90-million-year-old fossil, A. leggetti Hibbett, D. Grimaldi & Donoghue)57 and the divergence of the Hymenochaetaceae (based on the 125-million-year-old fossil, Q. cranhamii S.Y. Sm., Currah & Stockey58).

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 analysis59. 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 summarized60. 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.

Additional Information

How to cite this article: Song, J. et al. Phylogeny and biogeography of the remarkable genus Bondarzewia (Basidiomycota, Russulales). Sci. Rep. 6, 34568; doi: 10.1038/srep34568 (2016).