Phylogeny, biogeography and taxonomic re-assessment of Multifurca (Russulaceae, Russulales) using three-locus data

Multifurca is a small genus newly established to accommodate lactarioid and russuloid species with some characters reminiscent of corticoid members of Russulaceae. It shows an amphi-pacific distribution with strong preference for the tropical zone of the Northern Hemisphere and thus has particular significance for biogeographical study. Using worldwide samples and three loci (ITS, 28S rDNA and rpb2), we demonstrated that Multifurca is split into two highly supported major clades that are here recognized at the subgeneric level: subg. Furcata subg. nov. exclusively includes lactarioid species, while subg. Multifurca includes species with a russuloid habit. Using phylogenetic species recognition and comparison of genetic distances we recognize five new and six previously described species, almost double the known number of species before this study. Molecular dating using a Bayesian method suggested that Multifurca originated in early Paleocene and diversified in the Eocene. The most recent interspecific divergences occurred both in Asia and America, roughly at the same time around the Pliocene. Ancestral area reconstruction and comparisons of genetic distances and morphology suggested an early divergence within Australasia or tropical Asia. From the early Miocene to Pliocene, multiple dispersals/migrations to Australasia and North America by island hopping or land bridge likely happened. Vicariance at the late Tertiary might be the most likely mechanism accounting for the eastern Asia—southeastern North America and Australasia—tropical Asia disjunct distributions. The shared polymorphisms in the ITS alignment, numerous degenerated base pairs in the rpb2 sequences and weak conflict between the ITS and LSU genealogies of M. subg. Furcata suggest recent speciation. Host specificity of Multifurca species or species pairs is relatively low. Host shifts are believed to have aided establishment in new territories during the dispersals and migrations.

eleven M. furcata-like specimens, two of M. zonaria [19] and two of M. roxburghiae from central and southern China, three of M. aurantiophylla from New Caledonia, six of M. furcata and four of M. ochricompacta from central-southern and southeastern USA (Arkansas, Mississippi, Texas), and one of an undescribed species from eastern Australia were sampled to amplify the ITS region and 28S and rpb2 genes. In the DNA dataset of Lebel et al. [16], only ITS and 28S were sequenced for M. stenophylla and the ITS sequences were incomplete. We amplified the complete ITS region and rpb2 gene for four of the six specimens and one additional specimen of M. stenophylla collected from Tasmania, Australia (TL2462, MEL) to cover a broader geographical distribution.
The samples/data used in this study covered all previously described species of Multifurca and included five holotypes (M. orientalis, M. pseudofurcata, M. zonaria and two new species described below) and two epitypes (M. furcata and M. stenophylla). In addition, we retrieved three ITS sequences from GenBank by BLASTn: KP071201 labeled as R. japonica (from a Malaysian sample), FJ656009 from a New Caledonian sporocarp sample by Prin et al. (unpublished) and KR364125 from a Chinese sample labelled as Multifurca sp. [24]. These three sequences were included in the ITS dataset.
To test if the monophyly of Multifurca could still be retrieved after adding newly introduced taxa, seven species of Lactarius, 11 of Lactifluus and eight of Russula were added to construct a phylogeny. These taxa represented the various subgenera of these three genera based on the phylogeny of Russulaceae and DNA data published in former literatures [11,24,25] (S1 Fig).
After testing monophyly of the genus, we reduced the outgroup sampling to infer internal relationships within Multifurca.

DNA extraction, PCR amplifications and sequencing
For nine newly collected Chinese samples (less than ten years old), DNA extraction and PCR protocols followed Wang et al. [26]. For the other samples, DNA was extracted using a CTAB protocol [27] and purified with GeneClean1 II Kit (MP Biomedicals), following the manufacturer's instruction. The complete ITS region of the sample XHW2844 could not be successfully obtained by direct sequencing. For this sample, PCR product was cloned using the Takara1 pMD TM 18T cloning kit (Dalian, China) following the manufacturer's instruction. DNA insert of positive bacterial colonies was amplified and sequenced using the primer pair of the vector (M13F+M13R). Four clones with the desired length of PCR product were sequenced. Sequences of clones differing only in base substitutions were merged into one sequence by replacing the substitutions with degenerate bases; for clones with different INDELs, we present different sequences with these INDELs under the same sample (Table 1). The same cloning method was used to phase the heterozygosity of the ITS2 region of DPL4293 and to obtain sequences with weak PCR products.

Alignments and molecular phylogenetic analyses
Alignments were made using the online version of the multiple sequence alignment program MAFFT v7 [32], applying the L-INS-I strategy and manually adjusted in BioEdit [33]. To test if potential ambiguous sites in the ITS alignment would add noise to the phylogenetic analyses, we used different sets of parameters, from conservative to tolerant to select sites in Gblocks v0.91b [34] and compared the generated Maximum Likelihood (ML) trees and the ML Bootstrap Proportions (ML-BP) from each parameter set with the 28S and rpb2 trees and corresponding ML-BP. When testing the reciprocal monophyletic relationships of Multifurca, Lactarius, Lactifluus and Russula, only 28S and exons of rpb2 were used. For the subsequent analyses, we kept the alignable intron of rpb2 and used the three loci. In the concatenated ITS-28S-rpb2 dataset, we used duplicate sequences of 28S and rpb2 for the ITS copies of XHW2844 and DPL4293 respectively and removed the GenBank-retrieved sample xp2-20120922-01 of M. pseudofurcata because it has only ITS sequence (KR364125) that is identical with some of our sequences. The ITS, 28S and rpb2 alignments are deposited in TreeBASE as S21024.
To test for combinability between the ITS, 28S and rpb2 genealogies, we conducted ML bootstrap and Bayesian analyses on individual datasets sampling only the taxa for which the three loci were available. Phylogenetic conflict was assumed to be significant when two different relationships (one monophyletic and the other non-monophyletic) for the same set of taxa were both supported with ML-BP � 70% and posterior probability of Bayesian Inference (BI-PP) � 95%.
Phylogenetic analyses were conducted using ML in RAxML v7.2x [35] and BI in MrBayes v3.2.6 [36]. The dataset rpb2 was partitioned as follows: first and second codon positions, third codon position and intron. The ITS-28S-rpb2 combined dataset was partitioned as follows: ITS, 28S and the above three partitions for rpb2. ML analyses were conducted applying the Rapid Bootstrapping algorithm with 1000 replicates, followed by an ML tree search. For BI analyses, the best-fit models were selected using MrModeltest [37]. The BI analyses were conducted using four runs with four chains each for one million generations sampling one tree every 100 th tree. Runs were inspected to make sure the average standard deviation of split frequencies went below 0.01 and effective sampling sizes were > 200 in Tracer v1.7.0 [38]. A 50% majority rule consensus tree was built after discarding trees from the burn-in (first 25% of the trees). Trees generated by the two analyses were viewed and exported in FigTree v1.3.1. Trees were rooted with Lactarius pubescens.

Species delimitation using DNA data
Genealogical Concordance Phylogenetic Species Recognition (GCPSR) [39,40] was used to recognize the terminal independent evolutionary lineages, using only samples for which all three loci were obtained. A cluster of individuals was taken as an independent evolutionary lineage if it met with one of the two following criteria: i) the clade was present in at least two of the three ML genealogies; ii) its monophyly was highly supported by both ML-BP �70% and BI-PP �0.95 in at least one of the three genealogies and was not contradicted in any of the other genealogies at the same level of support. Exhaustive subdivision was then used for deciding which independent evolutionary lineages represented phylogenetic species [40]. That is to say, in the combined tree of the three loci, if any individual was not included in one of the terminal evolutionary lineages, we traced down the nodes of the tree from that individual until all individuals were included in evolutionary lineages.
To determine if the singletons "2676" and "RH10009" (that GCPSR does not fit) represent real species, we calculated two sets of genetic distances: within-group mean distances of all the species recognized by GCPSR and between-group mean distances between those two singletons and the other species. If the lowest between-group mean distance was bigger than the highest within-species distance, we separated them into two species, otherwise we merged them into their sister species. For "2676", we used ITS data and for RH10009 the combined ITS-28S-rpb2 data to calculate mean distances in MEGA5 [41], using the following settings: Bootstrap method (1000 replicates) as the variance estimation method, Maximum Composite Likelihood model as the substitution model, Gamma distribution as the rates among sites, different patterns among lineages, and pairwise deletion for gaps treatment.

Biogeographical analysis
We followed with minor modifications the methods of Dentinger et al. [42] and Feng et al. [43] to estimate the divergence times and infer the center of origin of Multifurca. Specifically, we used 28S and rpb2 datasets composed of 26 Multifurca and 36 outgroup taxa (S2 Table). Rhizopus oryzae Went & Prins. Geerl. (Mucoromycotina) was used to set the time of the most recent common ancestor (tMRCA) of Ascomycota and Basidiomycota, four species to represent the major clades of Ascomycota and 31 species representative of the major clades of Basidiomycota [44,45] [49,50].
For the rpb2 alignment, homogeneity of translated proteins was referred to improve the nucleotide alignment. The only intron of rpb2 was excluded from analyses. BEAST 1.8.0 [51] was used to estimate divergence times. To produce the executable xml file by BEAUTi (implemented in BEAST), we set parameters as follows: 28S and rpb2 datasets were set as two partitions, with substitution and molecular clock models unlinked and trees linked; taxa were constrained to monophyletic sets by referring to the SSU+28S+rpb2 or nrDNA+rpb2+ef1a phylogenies of fungi [44,52], the multi-gene phylogenies of Lactarius [46], Lactifluus [24] and Russula [11] and our three-locus phylogeny of Multifurca; Following the outcomes of model selection, GTR+I+G model was used as the substitution and site heterogeneity model and best frequencies were set as estimated; a lognormal relaxed clock was used for molecular clock analysis; Yule process was used as speciation for tree priors; ucld.mean of the 28S and rpb2 datasets were set to 0.033, with upper value = 1.0 and lower value = 0; for fossil calibration, we set mean = 575 Ma with a standard deviation of 38.26 Ma for the tMRCA of Ascomycota and Basidiomycota, inferred from the estimated time range of 500-650 Ma of the stem base of Ascomycota [53] based on the well-preserved ascomycetes fossil Paleopyrenomycites [54,55]. The reason why we selected this more ancient time origin instead of the more recent one of 452 Ma inferred by previous studies [56,57] is that the distribution of Multifurca in the South Pacific seems to suggest a Gondwana affinity and an earlier origin will be in favor of testing such an assumption. Two independent runs were conducted in BEAST for 20 million generations. Logfiles of the runs were inspected in Tracer v1.7.0 [38], to monitor if the ESS values of all parameters exceeded 200. Logfiles of the two runs were combined in Tracer v.1.7.0 by setting 10% logs as the burn-in. Tree files of the two runs were combined by in LogCombiner v1.8.0 [51] by setting 10% as the burn-in and processed in TreeAnnotator v1.8.0 [51] and then viewed and exported in FigTree 1.3.1.
To infer the center of origin for Multifurca and its species, we selected Multifurca from the chronogram generated by BEAST above for further analysis. We used ML analysis and Dispersal-Extinction-Cladogenesis (DEC) model implemented in LAGRANGE [58] provided in RASP v3.2 [59] to do the ancestral area reconstruction (AAR). The geographical areas of the extant species were defined as three regions: Australasia (A), Mainland and southeast Asia (B), and North and Central America (C). In LAGRANGE, for range constraints, we set the maximum areas as two and excluded the constrained area AC considering the two individual regions have never been adjacent. Dispersal probabilities between different areas were constrained following Clayton et al. [60], except that for the period 0-5 Ma we set the dispersal probabilities between Asia and America as zero, considering the cold climate at that time and the thermophilic habit of Multifurca species.

Morphology
Descriptions of sporocarps were from fresh material, and micro-morphological study was performed on dried material. Basidiospores were observed and measured in Melzer's reagent in side view, excluding ornamentation and apiculus. All other micro-morphological structures were revived in 5% KOH then mounted with Congo-red (aqueous reagent). All drawings except those of the basidiospores were made with the aid of a drawing tube installed on a Nikon E400 microscope. Colour codes refer to Kornerup & Wanscher [61] and abbreviations of herbaria follow Index Herbariorum. Measurements of spores were given as follows: (MIN) mean-SD-mean-mean+SD (MAX), and this for length, width and Q ratio; SD is standard deviation, while MIN or MAX is the absolute minimum, resp. maximum value measured. These values were only added (in between brackets) if they were lower, resp. higher, than the mean values minus, resp. plus, standard deviation.

Nomenclature
The electronic version of this article in Portable Document Format (PDF) in a work with an ISSN or ISBN will represent a published work according to the International Code of Nomenclature for algae, fungi, and plants, and hence the new names contained in the electronic publication of a PLOS article are effectively published under that Code from the electronic edition alone, so there is no longer any need to provide printed copies.
In addition, new names contained in this work have been submitted to MycoBank from where they will be made available to the Global Names Index. The unique MycoBank number can be resolved and the associated information viewed through any standard web browser by searching the MycoBank numbers contained in this publication to the website http://www. mycobank.org. The online version of this work is archived and available from the following digital repositories: PubMed Central and LOCKSS.  (Fig 1). In addition, all taxa of the lactarioid clade more or less share INDELs with at least one of the other species (Fig 1).

DNA alignments and characters
In comparison with the ITS dataset, the 28S and rpb2 alignments lack three samples, Multifurca sp 2676, K18C of M. aurantiophylla and xp2-20120922-01 of M. pseudofurcata due to unavailability of the other two loci (Table 1). All the samples of Multifurca have the same length of intron in rpb2, except that LPT1093 (M. pseudofurcata) is three base pairs shorter. Heterozygous loci in the rpb2 dataset are much more common than in the ITS and 28S datasets: among the 202 polymorphic sites of Multifurca, 69 sites are degenerated.

Phylogenetic relationships
The three ITS datasets with sites selected by conservative to tolerant parameters in Gblocks produced similar topologies (trees not shown). The dataset produced by the most conservative setting gave less resolved topology among the samples of M. furcata, M. mesoamericana, M. orientalis and M. pseudofurcata and among M. aurantiophylla, M. australis and M. zonaria and produced the lowest support values on nearly all the branches. In contrast, the support values produced by the dataset selected by the tolerant setting and un-Gblocked dataset are more comparable with that by the 28S and rpb2 datasets. It seems that the ambiguous sites in the ITS alignment did not cause much noise in the phylogeny and based on this we decided to use the original ITS dataset for further analyses.
MrModeltest selected HKY+G, GTR+I+G, SYM+G and SYM+G as the best-fit models for the ITS, 28S, rpb2 and ITS-28S-rpb2 datasets for the BI analyses. ML and BI analyses of the respective ITS, 28S and rpb2 datasets produced almost identical topologies with comparable support values. The ITS, 28S and rpb2 phylogenetic analyses inferred all identical deep and middle rank relationships but their respective topology solved different terminal relationships for M. furcata, M. mesoamericana, M. orientalis and M. pseudofurcata (refer to S2 Fig). However no significant conflict was detected between the topologies inferred by individual locus. The ML and BI analyses of the three-locus combined dataset produced identical topologies with comparable support values (Fig 2).
Analyses of the concatenated dataset produced a highly resolved phylogeny, with two major clades supported with ML-BP 100% and BI-PP 1.00 (A and B, Fig 2). The two major clades are represented by latex producing species (clade A) and species lacking latex (clade B) respectively. In clade A, samples of the Australian species M. stenophylla formed a basal clade, sister to the clade formed by all the Asian and American samples with an "M. furcata" morphology.  Bootstrap proportions higher than 70% in the ML analysis (ML-BP) and posterior probabilities of the Bayesian Inference (BI-PP) higher than 95% are indicated above and below the branches respectively or as ML-BP/BI-PP by the node. New species are in bold. Initials proceeding sample numbers refer to the collectors (see Table 1). as in the three-locus phylogeny, where they formed a clade sister to that formed by the North American M. ochricompacta and Asian M. roxburghiae. Compared to the species of clade A, relationships among all the species of clade B are well resolved. Distance analysis on the ITS-28S-rpb2 combined dataset showed that the within-group mean distance of clade B (0.035 ± 0.003) is almost double that of clade A (0.018 ± 0.002).

Species delimitation
Using the two criteria for determining evolutionary lineages in GCPSR, i.e. genealogical concordance and/or genealogical non-disconcordance [40], ten terminal independent evolutionary lineages were recognized. The tenth lineage, only supported by ITS analyses (ML-BP/ BI-PP 98/1.00), comprised four out of the eight M. pseudofurcata sampled in this study (S2 Fig: group 2  For the ITS dataset (S1 Text), the highest within-group mean distance (0.005) was for M. aurantiophylla. The lowest between-group value for the Malaysian sample 2676 (ITS: KP071201) is 0.022 (between M. mesoamericana) and based on this, we took it as a distinct species (Multifurca sp). For the ITS-28S-rpb2 dataset (S1 Text), the highest within-group mean distance (0.003) is for M. pseudofurcata group 2. The lowest between-group value for the Australian sample RH10009 is 0.037 (between M. zonaria) and based on this, we recognized it as another species (M. australis). In total, we delimitated twelve putative species in Multifurca using DNA data, five in M. subg. Multifurca and seven in M. subg. Furcata.

Divergence time estimation
Analyses calibrated with the stem base age of Ascomycota at 500-650 Ma estimated the stem base age of Multifurca at 60.32 ± 9.11 Ma (mean ± stdev) and diversification at 41

Ancestral area reconstruction (AAR)
The DEC model of LAGRANGE did not give a clear inference of a (single) center of origin for Multifurca with significant ML support (Fig 3). Instead, it gave three alternative ancestral

Synapomorphy, phylogeny and worldwide species diversity of Multifurca
Including worldwide samples in our morphological and phylogenetic analyses, we confirmed the synapomorphy of Multifurca to be a unique combination of three features as concluded by Buyck et al. [3]: i. a pileus context reproducing a much more distinct zonation than the one eventually observed on the pileus surface. Such a zonate context is very rarely known from other Russulaceae.
ii. regularly forking gills producing an orange spore print. Such regularly forking and yellow gills have also been reported from certain tropical Lactifluus [e.g. the Malagasy L. phlebophyllus (R. Heim) Buyck], while in Russula regularly forking gills are also found in some species of R. subg. Malodora Buyck & V. Hofst. However, none of these species produce an orange spore print.
iii. remarkably small spores with an inamyloid suprahilar spot. These spores are distinctly smaller than those formed by nearly all other agaricoid Russulaceae, yet comparable in size or even slightly larger than those formed by species of Russula subg. Archaea Buyck & V.
Hofst., which is still assumed to represent the oldest lineage in that genus [11].
In our three-locus phylogeny of Multifurca, the genus is composed of two fully supported major clades, in line with previous topologies [3,16,63]. These two clades separate all lactarioid, i.e. latex-producing, species from those that do not produce latex and could thus be described as russuloid. To accentuate the high similarity of all Multifurca species, including features of below-ground structures, we prefer to maintain a single genus and recognize each major clade as a distinct subgenus. Our choice of maintaining a single genus is also supported by the placement of M. zonaria in subg. Multifurca, notwithstanding the fact that its morphology would indicate a position amongst the lactarioid species.
GCPSR delimited ten phylogenetic species and analyses on genetic distances added two species in Multifurca. In total, twelve species were recognized with DNA data. In M. subg. Multifurca, morphological differentiation supported all five DNA-based species. In M. subg. Furcata, we kept five of the seven DNA-based species and merged the two groups of M. pseudofurcata into one species. The reasons for this "two-in-one" decision are: 1. Although M. pseudofurcata group 1 met the criterion of an independent lineage (genealogical non-disconcordance, Table 2), they did not form a monophyletic clade in the ITS-28S-rpb2 tree (S2 Fig); 2. Two different ITS sequences were obtained for sample M. pseudofurcata XHW2844, one occupying the most basal position in M. pseudofurcata clade while the other clustering in group 2 (Fig 2). This suggests that recombination may still possibly occur between group 1 and group 2; 3. We could find neither morphological differences between the two groups nor geographical structure within M. pseudofurcata. Inversely, we think there is sufficient genetic divergence and geographical segregation between the two lactarioid DNA-based species in America, M. furcata and M. mesoamericana, to recognize them as two species considering the monophyletic relationship of the two groups lacked significant support in all ML analyses and was only supported in the BI analysis (S2 Fig). The worldwide number of Multifurca species recognized in this study is eleven. Each continent has its endemic species.

Tropical Asia or Australasia as early divergence center of Multifurca
AAR did not infer a single center of divergence for Multifurca (Fig 3). This is most likely due to the unavailability of the origin center of its sibling genus Lactarius and the balanced distribution of its component major clades. However, North and Central America was excluded from being the divergence center and Australasia and Asia received medium (68.6%) to high (89%) marginal likelihood respectively. Analysis of genetic distance in each area showed that Asia and Australasia harbor comparable genetic diversity, higher than North America. In both major clades of Multifurca, there is a gradient of genetic distances across the known localities, i.e., from Australasia across Asia to North America, genetic distances decrease. Morphologically, species from Australasia and Asia show higher morphological diversification than those in North America. For instance, the Asian M. zonaria is the only species in the genus that has both hymenial pseudocystidia and gloeocystidia [3,14,19]. New Caledonian M. aurantiophylla is the only species in the genus that has mucronate pileocystidia [3,15]. In the lactarioid clade, the Australian M. stenophylla is the only one that can be morphologically recognized by the larger spores with a well-developed ornamentation [16]. Moreover, in this clade an ITS sequence of a Malaysian sample (GenBank accession KP071201) shows relatively bigger genetic distance than those between the other Asian species and even Asian and American species. All these facts make us believe that tropical Asia, Australasia or somewhere between the two regions might be the key region in the early divergence of Multifurca and its two major clades. The final answer to this question will probably depend on an ancestral area reconstruction of its sister genus Lactarius or even on a family wide approach.

Post Cretaceous migrations, dispersals and recent vicariance as possible mechanisms shaping the distribution of Multifurca
Multifurca is found under Nothofagaceae in New Caledonia and can also occasionally be found with this same host family in Australia [3,16]. However, it was never found in Patagonia, even though these Argentinian Nothofagaceae forests have been explored by quite a number of mycologists. There is also no record of Multifurca in Africa or Madagascar. A scenario that Multifurca has been extinct in this part of the world seems unlikely given the high local diversity of other ectomycorrhizal genera with a similar tropical distribution such as Lactifluus [24] and Cantharellus [64]. We can therefore safely conclude that species of Multifurca display a typical amphi-Pacific distribution with each continent having its endemic species (Fig 4).
Halling et al. [23] suggested three possible hypotheses that could account for such a distribution: i) long distance dispersal by basidiospores, ii) post-Cretaceous migration of co-symbionts over land bridges with a change or shift in symbiotic partners, and iii) an ancient, Pangaean distribution with little migration since the breakup of Pangaea in the Cretaceous. Previous molecular studies have suggested that these three scenarios are all possible, e.g. Pisolithus and some species pairs of Hysterangiales for the first hypothesis, Solioccasus and Hydnum p.p. for the second, and Cyttaria and Tuberaceae for both the first and third) [65][66][67][68][69]. Our analysis using a more ancient calibration point (500-560 Ma for the stem base age of Ascomycota) estimated the origin of Multifurca to be centered at the early Paleocene, well after the breakup of Gondwana. Our estimated ages of Multifurca and its relatives in Russulaceae and Russulales are slightly older than those given by Looney et al. [63] and Wisitrassameewong et al. [46], although those of the mycenoid and suilloid clades, Agaricales and Boletales are largely in line with the mininum age of the corresponding fossils and estimated ages of these orders using genome data [47,48,70]. Even with these older ages, the third hypothesis, i.e, an ancient, Pangaean distribution with little migration still does not currently fit Multifurca.
The decrease in genetic distance from Australasia across Asia to North America argues against a panmictic distribution, and thus against random long distance spore dispersal as the sole mechanism (hypothesis one). Following our AAR results, the corresponding times of divergence and the symbiotic relationships with hosts, we have to accept a complex mechanism that incorporates both medium distance dispersal and post-Cretaceous migration to explain the current distribution pattern.
Multifurca comprises two noteworthy biogeographic disjuncts in both subgenera: one between mainland and southeast Asia and Australasia, and the other one between East Asia and eastern North America (Figs 3 and 4). Multifurca zonaria, M. aurantiophylla and M. australis formed a New Caledonia-Australia-mainland and southeast Asia disjunct pattern, which is an intriguing biogeographic topic [71]. Our AAR inferred (with probability 100%) that the distribution of M. aurantiophylla in New Caledonia is due to a dispersal event from the joint area of Australasia and mainland and southeast Asia to the current territory (Fig 3). This inference is in line with the consensus that New Caledonia was separated from other land masses early in the upper Cretaceous and has never been connected with them again. Similarly some published data may suggest the same mechanism in ectomycorrhizal fungi between the islands in the South Pacific and other landmasses. For instance, in Lactifluus, some samples from American Samoa are found as conspecific with a Papua New Guinean environmental sample [72]. Wang et al. [25] showed that in the Lactifluus leoninus lineage, southern China-Thai L. aff. leoninus and a Papua New Guinean environmental sample shared 97% similarity in ITS sequences. In the genus Hydnum, it was found that one sample from Papua New Guinea is nearly identical to some Chinese samples [70]. More recently Han et al. [73] showed at least three independent dispersal events for ancestors of Strobilomyces from Southeast Asia to Australasia. Considering the relatively low host specificity of Multifurca species or species pairs [3,16,19] we postulate that the dispersal of M. aurantiophylla might be achieved by island-hopping with the aid of host shifts (Fig 4). Similarly the disjunct distribution between the Australian M. stenophylla and the rest species of M. subg. Furcata might be also by means of islandhopping. The northward movement of the Australian landmass during the Miocene is believed to have facilitated these dispersals.
Another interesting biogeographical topic involves the disjunct distribution pattern of the genus between East Asia and eastern North America [20,21,74]. In Multifurca, there are two species pairs/complexes that have such a distribution pattern. For the russuloid species pair M. roxburghiae-M. ochricompacta, AAR inferred the disjunct distribution to be due to vicariance from an ancestral population widely distributed in Asia and America. For the lactarioid species complex M. furcata-M. mesoamericana-M. pseudofurcata-M. orientalis, AAR suggested a first dispersal event from a widely distributed population in Asia-North America to Asia (for M. pseudofurcata) and then a vicariance from an ancestral population in Asia and America. This inference must be strongly influenced by the sibling relationship between American M. furcata and M. mesoamericana and Asian M. orientalis. Actually this topology was only produced by BEAST analysis and never received support in any other phylogenetic analysis. Therefore we would suggest this topology is an artifact and the relationship among these four species is actually equivocal. Taking account of such relationship and the shared INDELs and degenerated base pairs among these species (Fig 1), we would also like to suggest a scenario that a first vicariance between Asia and North America and a subsequent but very close migrations between eastern and western parts of Asia and southeastern North America and Central America respectively. In both species pairs/complexes, the widely distributed ancestral Asia-North American populations seem to have formed at the time between Miocene to Pliocene. The relatively warm climate in Miocene [75] and the presence of Beringia at that time [21] are believed to keep the Tertiary populations possible. The time of vicariance between Asia and North America was estimated to start around at the end of Tertiary, when the warm and humid climate was cooling down. The vicariance and migration must have taken place very fast, so that we see quite many shared polymorphisms in the ITS data (Fig 1), the weak conflicts between the ITS and 28S genealogies (S2 Fig) and the numerous degenerated sites and poor differentiation in the rpb2 data. This could also explain why in the American-Asian M. furcata complex there is no absolute split between Asia and Eastern North America, quite different from the findings in plant examples [76][77][78]. In Russulaceae there are at least two more groups that show such distribution pattern: Lactifluus sect. Gerardii (A.H. Sm. & Hesler) Stubbe and Russula compacta species complex [7,79]. Moreover in L. sect. Gerardii, similar equivocal relationship was found in some terminal clades [7]. The biogeographical analyses of the species pairs of Multifurca here will be of special reference for these biogeographically interesting groups.

Ectomycorrhizal associations of Multifurca species
The presumed ectomycorrhizal nature of the species composing Multifurca is still circumstantial and, notwithstanding several efforts by one of us (BB), dual symbiotic root structures have not yet been demonstrated, but the morphology of rhizomorphs emanating from the stipe base confirm the typical features of the family [11]. All the Multifurca specimens studied in this work were collected near putative ectomycorrhizal trees: in the northern hemisphere mainly Pinaceae and Fagaceae and to a lesser extent Dipterocarpaceae, while Nothofagaceae and Myrtaceae were likely tree hosts in the southern hemisphere. The Chinese specimens were all found under mixed forests with trees of Pinus subg. Pinus and Fagaceae or in pure pine forest, such as M. roxburghiae, which was originally described from Pinus roxburghii forests in Himalayan India [3] and re-collected in this same habitat at the moment of submission of this paper (K. Das pers. comm.). Multifurca zonaria was associated with fagaceous trees (Castanopsis spp.) in China [19], but with dipterocarp hosts in Thailand [14]. The American M. furcata and M. ochricompacta were all collected in fagaceous forests dominated by Quercus spp., although for both North American taxa, beech and in some cases also pine were close by [18] (Buyck pers. obs.). Nothofagaceae were the likely hosts for the New Caledonian M. aurantiophylla and for some collections of Australian M. stenophylla, although Myrtaceae (including Eucalyptus spp.) were clearly associated with some other collections of the latter, and apparently the only possible host for the single specimen of M. australis described here. It is most parsimonious to infer that Multifurca is a putative ectomycorrhizal group, using these two evidences: 1). Both our phylogenetic analysis and previous phylogeny of Russulaceae [3,7,9,24] show that Multifurca is nested in the core clade of Russulaceae that exclusively includes ectomycorrhizal species; 2). All the tree families/genus mentioned above have ages older than Multifurca [71,80,81].
Although based on the available data we are not able to infer the ancestral host of the genus and its major clades, the remote relationships of the putative hosts of the species pairs and even the same species (i.e. relatively low host specificity) strongly suggest host shifts during the diversification of the genus, which have been reported in many ectomycorrhizal groups [64,82,83]. If the symbiotic relationships of Multifurca species and these trees are eventually clearly demonstrated, the medium distance dispersals and migrations among America, Asia and Australasia and the present amphi-pacific distribution pattern of Multifurca will be better understood. Fig 5A-5E Although the species of this subgenus can be considered 'russuloid', as they do not exude latex upon injury, they do share several features that make them more similar to Lactarius, rather than to Russula. Indeed, the often zonate cap and scrobiculate stipe are characteristic of Lactarius, and so is the lamellar trama that is mainly composed of hyphae whereas sphaerocytes are dominant in the lamellar trama of Russula and most Lactifluus. Multifurca zonaria even has unequivocal hymenial pseudocystidia. Multifurca aurantiophylla is the only species of the genus that has well-differentiated pileocystidia, similar to many species of Russula.

Multifurca subgenus Multifurca,
The hymenial macrocystidia in Multifurca subg. Multifurca (except M. zonaria) are of two types. The first type concerns very voluminous cystidia that are deeply imbedded in the lamellar trama where these may sometimes inflate remarkably and emerge at the surface through a long neck. This type is unique among agaricoid Russulaceae as they are identical to those of some corticioid Gloeocystidiellum species. The second type is 'normal' minutely mucronate hymenial gloeocystidia that are so typical for most Russula. This second type can sometimes be extremely rare and is therefore easily overlooked, e.g. in some specimens of M. ochricompacta and M. roxburghiae. It is interesting to note that the minutely mucronate gloeocystidium type is absent from other parts of the fruiting body, but becomes abundant in the below-ground structures of all species studied in this respect: M. aurantiophylla [11], M. australis (this paper), M. zonaria (this paper) and even in the lactarioid M. furcata (this paper). Commentary: This American species might be the most (locally) common species in the genus. One of us (BB) has found it frequently in bottomland hardwoods and Beech-Magnolia-Loblolly Pine slope forests. It differs from all other species in this subgenus by its lower spore ornamentation.
Additional species: 1. Commentary: This species is by far the smallest species of all Multifurca (Fig 5A), in our opinion a consequence of the toxic metallic soils in which the host trees are living, because the exceptional small size is an almost general phenomenon for most of the local ectomycorrhizal fungi [84]. This species often has a strongly eccentric to laterally implanted stipe, which again, might be related to the often steep to vertical slopes on which this species often sporulates. This is the only species in the genus that has well-developed pileocystidia. Diagnosis: Differs from M. aurantiophylla in the field because of the larger basidioma size, the more central stipe and its association with Myrtaceae; differs from all other whitish species of the genus because of the one-type gleoecystida that are never ventricose.
Spores shortly ellipsoid, small, (5.0) 5.1-5.36-5.6 (5.8) × (3.7) 3.9-4.17-4.4 (4.6) μm, Q = (1.14) 1.21-1.29-1.37 (1.44); ornamentation subreticulate, low (ca 0.2 μm high), predominantly composed of ramifying crests, but often also with occasional laterally prolonged, convex warts and subtle line connections, the whole forming a subreticulate to almost reticulate pattern of variable density; some clearly ill-developed spores lacking the subreticulate ornamentation but beset with widely dispersed and strongly amyloid, blunt, isolated warts or even droplet-like pustules; suprahilar spot relatively small, inamyloid. Basidia (35) 40-50 × 7-9 μm, Phylogeny, biogeography and taxonomy of Multifurca clavulate, four-spored, with rather slender, long sterigmata, 5-6 × 1 μm; basidiola slender, subcylindrical to clavulate; intermixed with many very slender, irregular, sometimes terminally branching cells that are clearly not basidiola but clearly visible at the hymenial surface. Hymenial gloeocystidia very abundant on sides and edges of the gills, varying from hardly emergent to projecting up to 20-30 μm beyond the basidia, robust and of very variable size, originating sometimes from deep within the trama, (40) 80-150 (250) × (6) 8-12 μm, fusiform, subcylindrical, often slightly constricted subapically and subcapitate to rarely appendiculate, with a usually short neck, shortly pedicellate, filled with coarsely refringent-amorphous contents organized around circular cavities, that become then slowly coarsely crystalline, partially blackening in SV. Marginal cells slender, in length similar to basidiola, but not wider than 6 μm, somewhat undulate or repeatedly constricted in outline, usually tapering toward the tip, sometimes diverticulate. Lamellar trama a dense tissue of slender and narrow, strongly branching, filamentous hyphae of similar diameter (ca 2-3 μm), intermixed with the protruding bases of the gloeocystidia and some dispersed oleiferous fragments that are strongly and irregularly inflated-nodulose or variously constricted, also with some rare cystidioid hyphae. Pileipellis one-layered, poorly differentiated, hardly gelatinized, entirely orthochromatic in cresyl blue, near the very pileus surface forming a relatively thick tissue of very densely interwoven and very narrow, hardly 3 μm wide, and regularly cylindrical hyphae with simple, obtuse-rounded extremities; there also thin-walled or with slightly refringent walls near the very tips, but towards the lower pileipellis and underlying context often distinctly coated with refringent, mucus-like substance; inside the hyphae often with sparse refringent granules or inclusions. Well-differentiated pileocystidia and lactifers not observed. Stipitipellis equally undifferentiated, but in its lower portion developing large trichoids composed of very long and slender, cylindrical hyphae of identical diameter, some long caulocystidia of ca. 4 μm diam. present in these trichoid fascicules, but particularly abundant below the surface in the lower stipe and there often very large (similar in diam. to hymenial cystidia) and with first granular-amorphous, refringent, then coarsely crystalline contents; thick-walled vessel-like hyphae of variable diam. and length present inside lower stipe, also continuing in onset of rhizomorphal tissue below, some empty, others filled with homogeneous dense substance. Oleiferous hyphae abundant, similar as in pileus trama. Ladder-like hyphae not observed. Clamp connections absent. Below the stipe developing mycelial tissue with long and slender, aculeate and minutely mucronate gloeocystidia, 4-5 μm wide.
Commentary: The hymenial gloeocystidia of M. australis are never ventricose, i.e. they do not inflate in the part that is embedded in the subhymenium or lamellar trama, while they inflate moderately in M. aurantiophylla, and strongly so in M. ochricompacta and M. roxburghiae. The Australian species differs further from M. aurantiophylla in the field because of the (much) larger basidioma size and more central stipe (Fig 5B), as well as in its association with Myrtaceae; microscopically also in the smaller size and lower ornamentation of spores, the more regular and more slender hyphal endings in the pileipellis and absence of minutely mucronate, small hymenial gloeocystidia. Although M. australis does not possess minutely mucronate gloeocystidia in its fruiting body, it produces these on the below-ground parts (mycelium, rhizomorphs). They are nearly identical in form to those of M. aurantiophylla but lack the secondary septation as observed in the latter (Fig 7C, compare with ref. 11). The gloeocystidia in the lower stipe of M. australis do possess distinct contents (Fig 7A and 7B) and are therefore clearly differentiated contrary to the pileus where these are not apparent.  Cap up to 9 cm broad, when young slightly depressed in the center and already irregular in outline with a lobed or wavy and strongly inrolled margin, becoming infundibuliform at maturity with the margin plane or narrowly revolute, not becoming striate; surface covered by a subtle, white, felted tomentum that can locally take web-like to almost membranous aspects over the disc but is distinctly hirsute-fibrillose at the margin, only moderately viscid when wet, yellowish ochraceous with locally salmon orange tints, with faint concentrical zonation. Gills decurrent, quite crowded, very narrow, < 2 mm high, repeatedly forking, off-white when young, then fleshy-ochraceous with a tinge of salmon. Stem 2-3 cm. long and about 1.3 cm. thick, firm, solid and tough but the central part often narrowly and irregularly hollowing out over most of its length, surface covered by the same white tomentum as the cap but densely pitted-scrobiculate with light yellowish to ochraceous spots of variable size, stipe base not rounded, but clearly irregular. Flesh white and firm, in the cap with distinct concentrical watery-gray zonation throughout, up to 7 mm thick near stem; Odor weak. Milk moderately abundant, white, then slowly a distinct glaucous gray and remaining this color indefinitely; moderately acrid. Spore print pale salmon.
Additional material examined: UNITED STATES OF AMERICA. Arkansas: Perry Co., near the Lake Sylvia Spur, Ouachita National Recreation Trail, at the crossing with Highway 324, across the highway from the paved parking lot adjacent to the Nature trail, GIS 34  Commentary: Although M. furcata shares the same habitat and area of distribution with M. ochricompacta, it is much rarer than the latter species. Buyck et al. [3] discussed M. furcata on the basis of sequence data obtained for M. mesoamericana, which is described as a new species below. Given the prevalence of cryptic species in the M. furcata species complex and even several ITS types of M. furcata in North America, to fix the precise application of the name M. furcata, epitypification based on DNA data is needed.
When studying the mycelium-stipe transition, we found long and slender pseudocystidia with apices reminiscent of the typical aculeate-mucronate gloeocystidia in below-ground parts of species of subg. Multifurca. However, we were unable to find well-developed rhizomorphs. So far, there are no other observations made for below-ground structures of species in subg. Furcata. Commentary: Some of the morphological differences (more strongly zonate and viscose) mentioned in the diagnosis of this Costan Rican species may in part be due to the different climatic conditions inside the mountainous oak forests where it is found. However, there appears to exist a clear color difference between the two species and the very bright yellow tinges that characterize most M. furcata collections have not yet been observed in M. mesoamericana.  (DOCX) S1 Text. Matrices of within-group and between-group genetic distances estimated using ITS and ITS-LSU-rpb2 sequence data.