Genetic diversity and population structure of endemic mushroom Leucocalocybe mongolica in Mongolian Plateau uncovered by EST-SSR markers

Abstract Leucocalocybe mongolica (S. Imai) X.D. Yu &Y.J. Yao is an endemic mushroom with a distribution in Mongolian Plateau. It is an important edible and medicinal mushroom which is often collected unscrupulously. Wild L. mongolica is facing an unprecedented threat and is in a critically endangered state. Exploration of the genetic diversity and population structure of L. mongolica is important for conservation biology. The main distribution areas were investigated. In total, 402 alleles of 223 individuals from 17 wild geographical populations were obtained based on the eight simple sequence repeat loci. The total number of alleles from each population ranged from 11 to 31 with an average of 24 per population. The average observed heterozygosity (HO) and expected heterozygosity (HE) ranged from 0.300 to 0.875 and from 0.188 to 0.621, respectively. Analysis of molecular variance (AMOVA) indicated that 14.19% of the total molecular variance could be attributed to differences between populations (p <.001), whereas 85.81% to differences within populations. The mean genetic distance (FST) value of 0.171 suggested low genetic differentiation among populations. The UPGMA dendrogram generated from the SSR data suggested that the 17 populations separated into three main sister clusters on the Nei's genetic identity coefficient 0.8. The 223 L. mongolica individuals were divided into three groups by STRUCTURE analysis. Keeping and enriching the population diversity is crucial for the conservation of L. mongolica. The data about the population structure and genetic diversity obtained in this study will facilitate the establishment of conservation measures. Abbreviations AMOVA: analysis of molecular variance; Dest: estimator of actual differentiation; EST-SSR: Expressed sequence tag- simple sequence repeats; Fis: Inbreeding coefficient; Fit: total inbreeding coefficient; Fst: fixation index; GST coefficient of gene differentiation; G’ST: standardized measure of genetic differentiation; HE: expected Heterozygosity; H0: observed heterozygosity; HS: mean genetic diversity within populations; HT: total genetic diversity for species; I: Shannon’s diversity index; km: kilometer; L. mongolica: Leucocalocybe mongolica; N: sample size; NA: number of average alleles; NE: number of effective alleles; Nm: gene flow; NP: number of private alleles; NR: allelic richness; NT: number of total alleles; PCoA: principal coordinate analysis; PIC: polymorphism information content; PPB: Percentage of polymorphic loci; TP-M13-SSR: tailed primer M13; UHE: unbiased expected Heterozygosity


Introduction
Leucocalocybe mongolica (S. Imai) X.D. Yu &Y.J. Yao is an endemic species in Mongolian Plateau. Photographs of its ecological environment are shown in Supplementary Figure S1. Leucocalocybe mongolica is distributed in the steppe of Inner Mongolia of China, Mongolia. It can be used to make delicacies and has high popularity throughout East Asia. Leucocalocybe mongolica was once called Tricholoma mongolicum S. Imai [1]. The genus of Leucocalocybe was established because of some peculiarities in the morphological and molecular characteristics of L. mongolica [2,3]. As the only member of a mono-species genus, L. mongolica is of great value. At the same time, the basidiocarp is used as a folk medicine by Mongolians and Evenks believe that the basidiocarps can be used in the treatment of botulism and wounds, detoxification and anesis of blood-heat syndrome [4]. Natural products of potential pharmaceutical value have been extracted from basidiocarps or fermentation mycelia of L. mongolica, such as ergosterol, ergosterol peroxide, polysaccharide and lectins, which have been tested for biological activity, with some encouraging results [5][6][7][8][9][10][11][12][13][14].
Owing to its special habitat, L. mongolica has a restricted area of distribution [15]. Lacking a substantial breakthrough in its artificial domesticated cultivation, the sustainable use of resources has been seriously restricted. Leucocalocybe mongolica forms fruitbody only within a short timeframe, from June to September, accompanied by copious rainfall. Because the habitat is grassland and belongs to the arid or semiarid areas, the precipitation here does not meet the demand for fruiting. The habitat of L. mongolica is facing severe challenges, such as grassland degradation for overgrazing and grassland desertification [16]. In addition, L. mongolica is constantly being gathered by the local people and petty dealers, as the L. mongolica is in short supply and cannot meet the demand in a commercial trade. Leucocalocybe mongolica faces an unprecedented threat with scarcer and scarcer wild resources. The mushroom is now extremely endangered.
Exploration of the genetic diversity and population structure is a key step in conservation biology towards accelerating the development of a reasonable conservation strategy. However, to our knowledge, few efforts in this direction have been made in L. mongolica. In this study, we used microsatellites, or simple sequence repeats (SSRs), to study the population structure and genetic diversity of L. mongolica. As molecular markers, SSRs possesses the advantages of high polymorphism, abundant informativeness, simplicity and presence throughout the genome. Here, we used a set of SSRs identified and characterized in our previous work [16].

Ecological investigation and sample collection
This study included 223 L. mongolica samples from 17 wild geographical populations. The geographical populations were separated from one to another by at least 40 km, except for population 5 (ZK). The geographic distance between the central points of the 17 populations is presented in Supplementary Table S1). The geographic distance between the central points of pop5 (ZK) and pop4 (XL) is 28.41 km. Pop5 is a special sampling site, because it is located in the Inner Mongolia Grassland Ecosystem Research Station of the Chinese Academy of Sciences and as an isolated sample plot for continuous scientific research. All samples were collected in China and Mongolia from July 2015 to August 2017 (Table 1) and preserved in the Herbarium of Mycology of Jilin Agriculture University (HMJAU), under the Engineering Research Center of Chinese Ministry of Education for Edible and Medicinal Fungi in Jilin Agricultural University. The closest distance between individuals was supposed to be greater than 30 m.

DNA extraction and PCR amplification
Genomic DNA of basidiocarps was extracted using the Rapid Fungi Genomic DNA Isolation Kit (Sangon Biotech Shanghai Co. Ltd). A Quawell micro volume spectrophotometer Q5000 was used for detecting the DNA concentration and purity. DNA integrity was analysed through electrophoresis on 1.2% agarose gels. The DNA quality requirements included optical density of OD 260/280 !1.8, an electrophoresis result with a clear band and DNA concentration of 4 ng/mL or more.
Previously, 1860 SSRs were identified as molecular markers for L. mongolica by de novo assembly of the transcriptome [16]. One hundred and eight SSRs were  Table S2). SSRs with a tailed primer M13 (TP-M13-SSR) [17][18][19] were used in detecting the size of the PCR products. The PCR amplification system and program of fluorescence primer are described in Tables 3 and 4. The PCR products were analyzed by ABI 3730xl DNA analyzer (capillary type: ABI 4331246,50 cm,96 channels).

Data analysis
The parameters related to genetic diversity, including N A (number of average alleles), N E (number of effective alleles), I (Shannon's diversity index), H O (observed heterozygosity), H E (expected heterozygosity), U HE (unbiased expected heterozygosity), F is (inbreeding coefficient), F it (total inbreeding coefficient), F st (fixation index), PPB (percentage of polymorphic loci), H S (mean genetic diversity within populations), H T (total genetic diversity for species), G ST (coefficient of gene differentiation), G' ST (standardised measure of genetic differentiation), D est (estimator of actual differentiation) were calculated using GenAlEx v6.5 [20]. N T (number of total alleles) and N P (number of private alleles) were calculated by Arlequin V3.5.2.2 [21]. N R (allelic richness) was calculated by ADZE [22]. Pic (polymorphism information content) was determined using Cervus3.0.7 [23]. Gene flow (Nm) was computed by using the for- Eight loci were used to test whether they were suitable for the linkage disequilibrium and Hardy-Weinberg equilibrium in GENEPOP Web Version 4.0.10 [24,25]. Microchecker software was used to detect null alleles and typographic errors.
The genetic structure of 17 populations was analysed by STRCTURE 2.3.4 [26] with the admixture mode. K-values ranging from 2 to 11 were tested by setting the length of burnin period at 100,000 and the number of MCMC reps after burnin at 10000. The deltaK method [27] was implemented in Structure Harvester and the results were processed by CLUMPP [28]. The optimal clusters of populations were reflected on the map by the Pie function of R language.
Furthermore, Nei's genetic identity coefficient between all pairs of populations was determined by POPULATIONS Version 1.2.31 (http://bioinformatics. org/$tryphon/populations/). The dendrogram of the genetic relationship among the 17 populations was acquired based on Nei's genetic identity by applying the unweighted pair group method with arithmetic means (UPGMA) using NTSYSpc v2.2 [29].

Index of genetic diversity
The genetic parameters for Leucocalocybe mongolica wild populations are shown in Table 5. In total, 402     All quantitative values of the unbiased expected heterozygosity (U HE ) showed an overall trend similarly to the expected heterozygosity (H E ). The coefficient of inbreeding (F IS ) ranged from -1 to 0.187. The percentage of polymorphic loci (PPB) can be used as a marker to reflect the genetic diversity richness: DW, AB, XL and ZK were with lower richness than the average value 92.6%.
The characteristics on the eight microsatellite loci are given in Table 6. The numbers of total alleles per locus ranged from four (S11 and S87) to eight (S08) with an average of five per locus. The locus S08 had two private alleles; S03, S05, S09, S32 and S87 each had only one, while S11 and S59 had none. S08 had the highest average number of alleles (N A ), number of effective alleles (N E ) and Shannon's information index (I). Whereas H O ranged from 0.343 (S87) to 0.733 (S11), H E ranged from 0.359 (S32) to 0.577 (S08) and U HE ranged from 0.373 (S32) to 0.603 (S08). The inbreeding coefficient (F IS ) ranged from -0.276 to 0.198; the total inbreeding coefficient (F IT ) ranged from -0.135 to 0.299; the fixation index (F ST ), which reflects the degree of genetic differentiation, ranged from 0.109 to 0.259 with an average of 0.171 per locus, indicating a lower level of genetic diversity.
The genetic diversity reflects the adaptiveness of each species to environmental changes. High level of genetic diversity implies higher survivability to harsh natural circumstances. Individuals with more abundant variation will have greater variability of alleles associated with adaptation to the environmental changes. With changing environmental conditions, the population will spread on account of the individuals with more variation [32][33][34]. The mean Fst value of 0.171based on 8 loci in L. mongolica suggests low genetic difference among the 17 populations. It is lower than that of Armillaria luteovirens (Fst ¼ 0.176) in Qinghai-Tibet plateau [35] and much lower than that of Laccaria amethystina (Fst ¼ 0.416) between Japanese populations and European populations [36].
Nm at at each locus ranged from 0.716 (S32) to 2.04 (S59), with an average of 1.372 per locus. The highest polymorphism information content (PIC) was 0.693 (S08), and the lowest value was 0.385 (S05). The mean genetic diversity within populations (H S ) ranged from

Genetic structure of populations
The genetic structure analysis based on the clustering algorithm by STRUCTURE suggested a possible sectionalisation was K ¼ 3 (Figure 1), which justified by the fact that the 223 individuals were divided into three groups (Figure 2). The integral distributional situation of the three groups was mapped to the 17 geographic populations on the map (Figure 3). The genetic identity coefficient between populations varied from 0.5249 (AB,XW) to 0.9641 (EQ,WS).
The UPGMA dendrogram (Figure 4) indicated that the 17 populations were divided into three main sister clusters based on Nei's genetic identity coefficient 0.8. The first cluster contained the populations AB and DW. The lowest value of Shannon's information index (I) was attributed to population DW (0.260) and AB (0.442). In addition, the unbiased expected heterozygosity (U HE ) also had the lowest value in these two populations (DW 0.205 and AB 0.284). They are marginal populations on the boundary of the distribution area in Xilin Gol grassland. Although the sample collection plan was designed so as the closest individual distance to be greater than 30 m, the sampling regions of the geographic populations AB and DW are quite restricted. They are isolated and small. The phenomenon of habitat fragmentation emerged. The population AB situated on the edge of Otindag sandy land, is also frequently plagued by ecological stress. The population DW located on the edge of Ujimqin flat grassland, suffers from stresses such as drought and overgrazing. Habitat fragmentation is mostly depicted as a landscape-scale process concerning both habitat loss and the division of habitat. Empirical studies so far reveal that habitat loss has significant, continuingly negative impacts on biodiversity. On the contrary, habitat fragmentation has relatively lesser effects on biodiversity [37].
The second cluster contained populations XL and XW. They are geographically neighboring populations with a similar population structure; their main distribution area is Xilin Gol grassland. This cluster occupied a larger area in the pie graphs ( Figure 3). Wild individuals from the same geographic origins possess similar genetic components [38].
The remaining 13 populations could be roughly assigned into a cluster by Nei's genetic identity coefficient 0.8 and the 13 populations can also be divided into 5 subclusters by Nei's genetic identity coefficient 0.89. The 5 subclusters were respectively   As another method to explore the relationships among the 17 populations, we conducted principal coordinate analysis (PCoA). The two principal axes explained 49.05% and 82.18% of the total variance, respectively ( Figure 5). The Mantel test ( Figure 6) revealed no statistical significance between the genetic distance (F ST ) and geographical distance (r ¼ 0.049, p ¼.28); in other words, the genetic distances (F ST ) were not correlated with the geographical distances between the populations of L. mongolica. According to the UPGMA dendrogram (Figure 4), the clusters of populations poorly corresponded to the geographical distance. However, the genetic diversity, population structure and habitat conditions are considered more relevant. When the habitats of L. mongolica are in a state of destruction, ecological stress can cause changes of population structure. Some allele loss may occur  under the influence of unfavourable environmental conditions. Ecological stress plays a vital role in founder effect when habitat degradation causes ecological fragmentation. This will accelerate the change of population structure. We speculate that such environmental change can drive natural evolution in fungi as previously reported in annual plants [39].
Likewise, AMOVA analysis grouped by lineages indicated that 11.57% of the molecular variance was between groups and only 7.9% was between populations. For total populations, 14.19% of the total molecular variance was attributable to variance between populations (p < 0.001), whereas the rest (85.81%) was due to differences within populations (Supplementary Table S3). A similar proportion was observed in Tricholoma matsutake: 80% and 20% of the overall genetic variation were respectively partitioned within and among populations [40]. The vast majority of total molecular variance existed within populations. It could be because of a relatively greater gene flow (Nm ¼ 1.375) between L. mongolica populations (geographic distance between populations of 28.41-771.29 km). The distribution range of L. mongolica is limited. In this region, the basidiospores of L. mongolica are competent for transport by wind and living beings, so as to decrease the genetic differences among populations. That could further explain why the genetic distances (F ST ) were not correlated with the geographical distance of the populations in our study. Figure  S2 shows the average pairwise differences in the populations.
Macrofungi, due to their distinctive life cycle, grow for the larger part of their life cycle in the form of mycelium in deep, humus-rich soil form. It is difficult to make a clear-cut distinction between the population and an individual. The macroscopical basidiocarp occurs only for a brief period, for the purpose of sexual reproduction. Dicaryophytic mycelium growth and the dissemination of basidiospores, all contribute to the growth of mushroom populations [41]. The definition of a population and an individual mushroom is a pendent challenge in mycology. In our survey, the sampling scheme was developed in such a way that the closest basidiocarp distance was supposed to be greater than 30 m. However, since most of the fungal thallus consists of a network of anastomosing hyphae embedded in the substratum, the clonal "individual" is difficult to define [42]. After studying the natural population of the Oyster Mushroom Pleurotus ostreatus, the basidiospore dispersal genetic mechanism needs further studies to determine the extent of different mating patterns between individuals in the population [43].

Conclusions
In this study, out of 108 SSRs, 8 polymorphic SSRs were detected by acrylamide gel electrophoresis. The products amplified by 8 polymorphic SSR primers were separated by capillary electrophoresis as a more accurate test method suitable for the size of PCR products. The EST-SSRs were successfully applied to reveal genetic diversity and population structure of Leucocalocybe mongolica. Keeping and enriching the population diversity is crucial for the conservation of L. mongolica. The obtained data will accelerate the establishment of appropriate conservation measures.