Phylogenetic study documents different speciation mechanisms within the Russula globispora lineage in boreal and arctic environments of the Northern Hemisphere

The Russula globispora lineage is a morphologically and phylogenetically well-defined group of ectomycorrhizal fungi occurring in various climatic areas. In this study we performed a multi-locus phylogenetic study based on collections from boreal, alpine and arctic habitats of Europe and Western North America, subalpine collections from the southeast Himalayas and collections from subtropical coniferous forests of Pakistan. European and North American collections are nearly identical and probably represent a single species named R. dryadicola distributed from the Alps to the Rocky Mountains. Collections from the southeast Himalayas belong to two distinct species: R. abbottabadensis sp. nov. from subtropical monodominant forests of Pinus roxburghii and R. tengii sp. nov. from subalpine mixed forests of Abies and Betula. The results suggest that speciation in this group is driven by a climate disjunction and adaptation rather than a host switch and geographical distance. Electronic supplementary material The online version of this article (10.1186/s43008-019-0003-9) contains supplementary material, which is available to authorized users.


INTRODUCTION
Russula is a cosmopolitan genus of basidiomycetes comprising hundreds of species with a mainly agaric habit of the basidiomes (Looney et al. 2016). Members of the genus are important ectomycorrhizal (ECM) partners in forest ecosystems and are also used as commercially traded edible fungi (Looney et al. 2018). A study of global datasets of Russula shows only small overlaps in species diversity between geographically distant areas in the Northern Hemisphere (Looney et al. 2016). Current studies suggest that Russula species reported from multiple continents are mainly boreal-arctic taxa. Bazzicalupo et al. (2017) found that 17 of 72 Russula species from the Northwest of the USA (mostly from boreal-arctic areas) matched a sequence of European barcoded taxa. Accordingly, Geml et al. (2012) suggested that long distance dispersal and gene flow between boreal-arctic species of the Russulaceae (Russula and Lactarius) was common. A phylogenetic study of subsection Xerampelinae (Adamčík et al. 2016a) has demonstrated that R. subrubens and R. clavipes occur not only in distant areas of Europe, but that they also inhabit a wide range of ecosystems in very different climatic areas (from the temperate belt to boreal-arctic-alpine regions).
Little is known about the phylogeography of ECM fungi in Eurasia. For a long time, European Russula species were believed to be present in East Asia. Hongo (1960), for example, listed only four native Asian species among 40 Russulas reported from Japan, and 31 species reported from the country were originally described from Europe; the remaining species were from North America. Cao et al. (2013) demonstrated a 3-6% ITS nrDNA sequence divergence between Chinese collections and the most similar European species R. virescens, but despite the generally accepted 97% sequence similarity cut-off value for the recognition of species (Ryberg 2015), they still hesitated to describe them as a new species. More recently, there has been a considerable increase of new Russula species described from Southeast Asia. For example, 24 new Russula species have been described from China in the period 2005-16 (Li et al. 2016), and an additional eight new species have been recognized there in 2017, which indicates an increasing annual rate (Das et al. 2017a, b). However, only two new Russula taxa have so far been described from higher elevations of the montane belt (2500-4000 m a.s.l.) of the southeast Tibetan Plateau: R. wangii and R. amethystina subsp. tengii (Li et al. 2016). Accordingly, in the area of the southeast Himalayas, most taxa have been described from lower elevations of up to 2300 m (Joshi et al. 2012), and only R. shingbaensis and R. thindii have been found in subalpine coniferous forests (elevations above 3000 m) with Abies densa (Das et al. 2014).
The montane belt of the Himalayan Mountains has developed an autonomous flora differing from Eurasian boreal taiga forests because of isolation by distance and climatic disjunctions (Miehe et al. 2015). The montane coniferous forests that flank the Himalayas at their southeast to southwest edge are not connected to the Siberian boreal taiga forests. The closest areas with a connection to Eurasian taiga forests are the Tian Shan and Altai Mountains to the northwest and Inner Mongolia to the northeast. Singer (1938) reported 33 Russula taxa from the montane, subalpine to alpine belts of the Altai and four Russula species from similar areas of Tian Shan. To his Asiatic collections from the area, he assigned 21 European Russula names and described three new species, five new subspecies, four new varieties and one new form. His alpine species R. oreina is considered a synonym of the European R. pascua (Knudsen & Borgen 1992a) and his R. mesospora (associated with Betula) a synonym of the European R. intermedia (Ruotsalainen & Vauras 1994). Russula citrinochlora, described from the subalpine belt of the Altai Mountains, has also been reported from Greenland (Knudsen & Borgen 1992b) and Fennoscandia (Marstad 2004). Only a few Russula species have been reported from Inner Mongolia (e.g. Tian et al. 2014) and only one, R. jilinensis (Li et al. 2012), has been described as new from the region.
The assembly of phylogenetic community structure and species pool scaling of forest communities in East Asia revealed the important role of intercontinental migration during the Neogene and Quaternary in the formation of species diversity in the area (Feng et al. 2015). Study of the global biogeographic distribution of Alnus-associated ECM fungi (Põlme et al. 2013) indicated co-dispersal of hosts and their mycobionts, but at a regional scale, geographic distance or disjunctions may also largely account for the observed diversity (Bahram et al. 2012).
In this study we attempted to test the hypothesis that geographic distance, a host switch or a climatic disjunction contribute to the evolution of ECM fungi using the Russula globispora lineage. This lineage belongs to the subsection Maculatinae and is morphologically defined by a slightly acrid taste of the flesh, yellow spore print, yellow-brownish spots on the basidiomata and spores with isolated large warts or spines (Adamčík & Jančovičová 2013, Adamčík et al. 2016b). In the R. globispora lineage, R. globispora and R. dryadicola are the only species reported from boreal, arctic or alpine habitats of Europe (Sarnari 1998). The aim of this study was to analyse phylogenetic relationships among collections and available sequence data of the R. globispora lineage originating from coniferous subtropical, mixed boreal and alpine forests of Europe, the southeast Himalayas, and North America.

Sampling
We analysed 10 collections of Russula globispora and R. dryadicola from various habitats across Europe, including temperate deciduous forests, boreal mixed forests and alpine habitats. Asian material from the south Himalayan Mountains is represented by four morphologically similar collections from subtropical coniferous rainforests of Pinus roxburghii in Pakistan, and four other similar collections from montane mixed forests of China associated with Abies spectabilis and Betula delavayi. Three representatives of subsection Cupreinae, identified and sequenced by us, were used as the outgroup. This study is supplemented by sequences with high (99%) similarity retrieved from the Genbank (https://www.ncbi.nlm.nih.gov) and UNITE (https:// unite.ut.ee) databases. All samples and sequences used are listed in Additional file 1: Table S1.

Molecular analysis
Total genomic DNA was extracted from dried material using a variety of protocols (Eberhardt 2012, Adamčík et al. 2016a). In addition to these protocols, the EZNA Fungal DNA Mini Kit (Omega) was used following the manufacturer's recommendations, but with prolonged incubation time of up to 1 h after addition of the RNA-lytic enzyme. Three molecular markers were amplified: (1) the internal transcribed spacer regions of ribosomal DNA (ITS); (2) the partial mitochondrial small subunit of ribosomal DNA (mtSSU); and (3) the region between domains six and seven of the nuclear gene encoding the second largest subunit of RNA polymerase II (rpb2). The ITS region was amplified using the primers ITS1F-ITS4 (White et al. 1990, Gardes & Bruns 1993. The mtSSU region was amplified using the primer pair MS1 and MS2 (White et al. 1990). The rpb2 was amplified using the primers A-Russ-F-frpb2-7CR (Matheny 2005, Caboň et al. 2017). All three molecular markers were amplified with Hot Start Firepol Polymerase (Solis Biodyne, Tartu, Estonia) using cycling protocols according to Caboň et al. (2017). The PCR products were purified using Exo-Sap enzymes (Thermo Fisher Scientific, Wilmington, DE) or the Qiaquick PCR Purification Kit (Qiagen, Hilden, Germany). Samples were sequenced by the Seqme company (Dobříš, Czech Republic).

Phylogenetic analysis
Raw sequences were edited in Geneious version R10 (Kearse et al. 2012). Intra-individual polymorphic sites having more than one signal were marked with NC-IUPAC ambiguity codes. All three single-locus datasets were aligned by MAFFT version 7 using the strategy E-INS-i (Katoh & Standley 2013), manually improved in Geneious version R10 (Kearse et al. 2012) and concatenated into one multi-locus dataset using Sea-View version 4.5.1 (Gouy et al. 2010). The multi-locus dataset was analysed using two different methods: Bayesian inference (BI) and the maximum likelihood (ML) approach. For the ML analysis, the concatenated alignment was loaded as a PHYLIP file into RAxMLGUI v. 1.2 (Silvestro & Michalak 2012) and analysed as a partitioned dataset under the GTR + I model with 1000 bootstrap iterations. In addition, the GTR + G model was used as recommended by Stamatakis (2008). Phylogenetic analyses generated using both models resulted in trees with identical topology and very similar bootstrap support. Thus, only the GTR + G + I analysis is presented. For the BI analysis, the dataset was divided into six partitions: ITS, mtSSU, intronic region 7 of rpb2, and the 1st, 2nd and 3rd codon positions of rpb2. The best substitution model for each partition was computed jointly in PartitionFinder v. 1.1.1 (Lanfear et al. 2012). The BI was computed independently twice in MrBayes version 3.2.6 (Ronquist et al. 2012) with four MCMC chains for 10,000,000 iterations until the standard deviation of split frequencies fell below the 0.01 threshold. The convergence of runs was visually assessed using the Trace function in Tracer version 1.6 (Rambaut et al. 2013).

Morphological analysis
Micromorphological characteristics were observed using an Olympus CX-43 microscope with an oil-immersion lens at a magnification of 1000×. All drawings of microscopic structures, with the exception of spores, were made with a 'camera lucida' using an Olympus U-DA drawing attachment at a projection scale of 2000×. The contents of hymenial cystidia and pileocystidia were illustrated as observed in Congo red preparations from dried material, with the exception of some pileocystidia for which the contents are indicated schematically (by plus signs). Spores were observed on the lamellae stained with Melzer's reagent. All other microscopic observations were made in ammoniacal Congo red after a short treatment in a warm aqueous KOH solution to dissolve the gelatinous matrix and improve tissue dissociation. All tissues were also examined in Cresyl blue to verify the presence of ortho-or metachromatic reactions as explained in Buyck (1989). The trama and cystidia were examined in a sulfovanillin solution. Acidoresistant incrustations of primordial hyphae were stained with carbolfuchsin and observed in distilled water after incubation for a few seconds in a 10% solution of HCl (cf. Romagnesi 1967). Spores were scanned with an Artray Artcam 300MI camera and measured by Quick Micro Photo version 2.1 software with an accuracy of 0.1 μm. Spore measurements excluded ornamentation and their line drawings were made using enlarged, scanned pictures. The Q value indicates the length/width ratio of the spores. Spore ornamentation density was estimated following Adamčík & Marhold (2000). Estimates of the density of cystidia estimates follow Buyck (1991). Statistics for the measurements of microscopic characteristics were based on 30 measurements per specimen and are based on all examined material of the described species. The range of measured values is expressed as the mean ± standard deviation; presented in parentheses are the 5th and 95th percentiles.

Phylogenetic analysis
In total, 54 sequences were newly generated and the final dataset was supplemented by 28 published sequences, including 20 ITS sequences with high sequence similarity to Russula dryadicola retrieved from public databases (Additional file 1: Table S1). They altogether correspond to 42 collections in the tree (Fig. 1), and 38 of them are placed in the R. globispora clade with full support. Three collections of R. globispora from temperate deciduous forests are placed at the basal position of the clade. Collections from all other habitats (coniferous and boreal forests or arctic-alpine habitats) are clustered in a clade with moderate ML and full BI support (70/1). This clade is composed of three well supported clades with various branch length, each corresponding to samples from different habitats and geographic areas. Our collections from Pinus roxburghii mono-dominant forests of the southwest Himalayan Mountains (Pakistan) are grouped in a fully supported clade and below we describe them as a new species bearing the name R. abbottabadensis. Another well supported species clade is formed by collections from mixed montane forests of the southeast Himalayas (China) and corresponds to the second new species, R. tengii, described below. The third clade groups samples that originated from a large geographic area of Europe and Alaska, collected in boreal, arctic and alpine habitats. Alaskan samples retrieved from the Genbank database received strong support (81/0.96), but they are represented only by ITS sequences and differ from other samples within this clade of European origin only in 1-2 parsimony-informative positions. All European collections correspond morphologically and ecologically to the concept of R. dryadicola.

Morphological analysis
All four species in the Russula globispora clade are morphologically defined by yellow-brownish spots on the surface of the basidiomata, yellow spore print and relatively large spores, which are characters that correspond to the circumscription of R. subsect. Maculatinae. In addition, all studied members of the R. globispora clade have large isolated spines on the spores and a weakly acrid to mild taste that distinguishes them from other members of the subsection. In the field, R. abbottabadensis is easily recognized by the bright red colours on the pileus showing only a little discolouration in spots near the centre, whereas the other two species usually have dull pink to violet colours and are strongly discoloured near the centre. The basidiomata of R. dryadicola are more robust and larger.
Under the microscope, the most important differences were observed in the pileipellis. Russula abbottabadensis has long (on average longer than 29.5 μm), narrow (on average to 3 μm), cylindrical terminal cells of hyphae in the pileipellis near the pileus margin with usually branched subterminal cells and narrower (on average to 5 μm), prevailingly 1-2-celled pileocystidia. Russula tengii has shorter (on average to 29.5 μm) terminal cells of hyphae in the pileipellis near the pileus margin, less branched subterminal cells, and wider (on average more than 5 μm) and more septate (mainly 2-or more-celled) pileocystidia. Russula dryadicola is well characterized by wider (on average more than 3.5 μm), frequently fusiform or lanceolate terminal cells of hyphae in the pileipellis near the pileus margin that are often moniliform.
There are no other species known from southeast Himalayas that remind morphologically or in sequence R. globispora lineage. However, according to our knowledge (K. Wisitrassameewong, personal communication), there are probably more undescribed members of the  Spores (7.9-)8.7-9.5-10.2(− 11) × (6.9-)7.6-8.  Diagnosis: Basidiomata small to medium sized, thickfleshed, similar to R. dryadicola in pileus colour and presence of yellow-brownish spots on surface of pileus and stipe, from which differs by mainly cylindrical, apically obtuse and not constricted, in average up to 3.5 μm wide hyphal terminations in pileipellis near the pileus margin and pileocystidia that are mainly two and more celled. Description: Basidiomata small to medium-sized, robust.
Pileus 28-78 mm diam, first hemispherical, then convex to applanate, sometimes slightly depressed; margin first slightly deflexed, when mature straight to slightly inflexed, sometimes undulate, indistinctly striate; pileus cuticle slightly viscous when wet; colour near the pileus margin light brownish vinaceous to light russet vinaceous, when old brick-red, and madder-brown when dry, near the pileus centre mustard-yellow to aniline-yellow when old and dry, variegated with dark blackish red and dark yellow spots. Lamellae relatively dense, 10-16 per cm near the pileus margin, adnate-emarginated, 4-6 mm broad, occasionally forked near the pileus margin, rarely near the stipe, slightly anastomosed, without lamellulae, first creamy to yellowish, when mature baryta yellow to buff yellow, pale yellow-orange when dry; edge even, concolourous. Stipe 35-47 × 12-15 mm, central or slightly eccentric, cylindrical, often swollen toward the base, solid to hollow when mature, smooth; white, soon with a buff yellow tint, especially near the base. Context 3-7 mm thick at the attachment of lamellae to the stipe, compact but fragile in gills, white, slightly turning pale creamy or yellowish when wounded or old. Smell indistinct. Taste mild. Spore print yellow (IVc-IVd according to Romagnesi 1967).
Diagnosis: Basidiomata small to medium sized, pileus red, stipe with yellowish brown rusty spots, taste mild to slightly acrid, spore print yellow, spores large, similar as in R. globispora but with occasionally fused spines, Hyphal terminations in pileipellis near the pileus margin with narrow (in average up to 3 μm), often attenuated terminal cells and mainly branched subterminal cells, pileocystidia narrow (in average up to 5 μm) and usually 1-2 celled. Description: Basidiomata small to medium-sized, relatively thin-fleshed. Pileus 13-73 mm diam, first convex, then plano-convex to applanate, sometimes slightly centrally depressed; margin first deflexed, when mature usually inflexed, entire, hardly striated; pileus cuticle matt when dry, sometimes viscid and shiny near the centre; colour moderate red, strong red, deep red, vivid red or vivid reddish orange, rarely with light orange-yellow spots near the centre. Lamellae adnexed or adnate-emarginate, brittle, without lamellulae; when young white to pale yellow, when mature yellow; edge even, concolourus. Stipe 16-35 × 7-12 mm, central, clavate or cylindrical, longitudinally rugulose, not pruinose, white, occasionally with pink flush, with yellowish brown or rusty spots, solid. Context compact, white, slowly turning yellow-brownish when wounded or old. Smell indistinct. Taste mild to slightly acrid. Spore print yellow.

DISCUSSION
This study confirms that Russula dryadicola is a well-defined species in accordance with Fellner & Landa (1993), contrary to some previous studies that treat it as an infraspecific taxon of R. maculata (Singer 1926, Knudsen & Borgen 1992b. The species was first reported as R. maculata subsp. alpina from the Alps (Singer 1925(Singer , 1926 and considered an alpine-arctic element distributed in Greenland, Scandinavia, mountains of Central Europe and the Ural Mts (Kühner 1975, Knudsen & Borgen 1992b (Sarnari 1998). We have studied seven collections sent by J. Vauras from similar boreal habitats in Finland, and they all morphologically correspond to R. dryadicola. In our phylogenetic tree, boreal collections are represented by TURA 152390 and TU 101835, and they are clustered with other alpine and arctic collections of R. dryadicola. This species is an example of alpine-arctic habitats being incorrectly applied for species delimitation and identification. Similarly, R. nuoljae, previously described and recorded from alpine and arctic habitats, was recently reported also from boreal forests, and R. subrubens, previously reported only from lowlands of the temperate belt vegetated by willows, was synonymized with the alpine-arctic R. chamitae (Adamčík et al. 2016a). As new morphological characters of R. dryadicola, differentiating it from R. globispora (the latter species is re-described in Adamčík & Jančovičová 2013), we introduce here shorter (on average to 33 μm long) and often fusiform terminal cells of hyphae in the pileipellis near the pileus margin, narrower pileocystidia (on average to 6.5 wide; Fig. 3) and narrower pleurocystidia (on average to 12.5 μm wide; Fig. 4). Alaskan samples from boreal forests have nearly identical ITS sequences and any decision about their taxonomic status requires support of more DNA loci or other molecular support. There are also two ITS sequences from environmental samples associated with the alpine sedge Kobresia myosuroides; these sequences are of unspecified origin and have apparently been sampled for the Niwot Ridge project (http://niwot.colorado.edu) dealing with research of permanent plots in the Colorado alpine ecosystem in the Rocky Mountains (USA). Nevertheless, this affinity of North American and European collections suggests a wide distribution of R. dryadicola in Eurasian taiga and tundra habitats, as demonstrated by Geml et al. (2012). North American ITS sequences from environmental DNA samples, placed in this study in the R. globispora lineage, are the first reports of the group from America, but not a single basidiome has been collected in that region to date. The BLAST result of the ex-holotype ITS sequence of R. tengii shows 99% identity (under 100% query coverage) with the sequence of the R. dryadicola collection TURA 151632. This is an example of high similarity in the ITS region within the R. globispora clade. However, all three species recognized in our multi-locus phylogenetic analysis (Fig. 1) received good support (BS ≥ 65, PP ≥ 0.98). This is analogous to the high ITS sequence similarity within the R. clavipes species complex (Adamčík et al. 2016a).
Russula tengii is represented only by collections from mixed subalpine forests of the southwest Himalayan Mountains in China. This type of habitat falls within the known ecological amplitude of R. dryadicola. The existence of allopatric boreal-arctic species within the R. globispora clade suggests that glaciation events and geographic and climate disjunctions were the evolutionary drivers of speciation in the group on the Eurasian scale. It is possible that the common ancestor of R. dryadicola and R. tengii migrated via available corridors for movement to newly emerging environments during climatic changes in the Quaternary (Donoghue 2008), but we do not know the original area of this ancestor species.
The collections of R. abbottabadensis under study originate from subtropical coniferous forests of Pinus roxburghii in northeast Pakistan at the foothills of the Himalayas. Coniferous trees seem to be obvious hosts of this species, and they possibly also host R. dryadicola and R. tengii, which grow in mixed boreal forests. The long branch of this species and lack of support for R. globispora suggest missing data on possibly unrepresented species of the R. globispora lineage from other environments and regions of Eurasia. It is also possible that gene flow and species migrations between Southeast Asia and Europe were multidirectional and that R. abbottabadensis may represent another lineage of so far unknown origin. Both Asian species described here represent, to our knowledge, the first reports of the R. globispora lineage from Southeast Asia and we hope that this study will boost more studies from different parts of the continent to answer questions on origin and speciation pattern of this Russula group. Looney et al. (2016) suggested that host switching plays major role for species diversification of the genus Russula especially at the rank of clades. However, the speciation mechanism of Russulaceae and even ECM fungi in general is not well understood especially within closely related species lineages. This study suggests that host switch did not contribute to speciation of studied species of R. globispora lineage and more probable drivers in the evolution of R. tengii and R. abbottabadensis were climate disjunctions and adaptations.

CONCLUSIONS
We recognised three species of R. globispora lineage that either grow in mixed forests in cold climate (boreal, montane, arctic or alpine) or are associated with conifers in mountains. Russula dryadicola is widely distributed in alpine, arctic and boreal areas of Eurasia and Northern America, although North American collections show little genetic divergence.