Mengli Han1†
Jingxue Zhang1†Daxu Li2Shengnan Sun1Changbing Zhang2Chuanjie Zhang1Lijun Yan2Yongzhuo Guan3Lili Chen2
Yuxia Guo3Minghong You2Wenlong Gou2
Xinrui Li2Jiajun Yan2Shiqie Bai2
Xuebing Yan1*- 1College of Animal Science and Technology, Yangzhou University, Yangzhou, China
- 2Sichuan Academy of Grassland Sciences, Chengdu, China
- 3College of Animal Science and Technology, Henan Agricultural University, Zhengzhou, China
Elymus sibiricus L. is a perennial allotetraploid belonging to Triticeae of Poaceae, Elymus L., as the type species of genus Elymus L. The existing geographical distribution pattern and genetic spatial structure of E. sibiricus on Qinghai–Tibetan Plateau (QTP) are not yet clear. In this study, population genetic structure and demography history of 216 individuals from 44 E. sibiricus populations on QTP were studied used specific-locus amplified fragment sequencing (SLAF-seq). The result of genetic diversity showed that there was no single genetic diversity center was observed across all E. sibiricus populations. The results of genetic variation showed that 44 populations were clearly divided into the following three groups: Qinghai Plateau (Group I), South Tibet (Group II), and Hengduan Mountains (Group III). From the three analyses of AMOVA, Mantel test and Treemix, strong genetic differentiation across all populations and low genetic differentiation among populations within three groups. Molecular dating indicated that E. sibiricus diverged at 16.08 Ma (during the early Miocene) can be linked to the Himalayan Motion stage of QTP uplift. It is speculated that the reasons affecting the current phylogeographical pattern are as follows: (1) The environmental changes due to the uplift of the QTP; (2) The geographic distance between the populations (Groups I and III are close in geographic distance, and gene flow are frequent); (3) Geographical barriers (the Tanggula and Bayangela Mountains between Groups I and II). This study provides new evidence and historical perspective to the future exploration of the evolution and geographic distribution pattern of Elymus L.
Introduction
Understanding the underlying phylogeographical processes of the species distribution in some areas may help guide prioritization in conservation. It also may facilitate forecasts on the ecosystem services under the future climate conditions. It is believed that geological events and climate fluctuations have profoundly affected the distribution of species and the history of population dynamics in mountainous areas (Hewitt, 2004). The global climate change, especially the simultaneous climate that has fluctuated since the Quaternary Ice Age, has had a severe impact on the distribution of species, further forming plant differentiation and the current distribution pattern (Tian et al., 2010; Guo et al., 2012). In general, the response of species to Quaternary climate change includes changes in scope by forcing organisms into sanctuaries and re-colonizing under more favorable conditions (Taberlet et al., 1998; Hewitt, 2000). In the past 21–26 years, a large number of systematic geography studies in Europe and North America have summarized the evolution history of Quaternary flora in these regions (Abbott, 2000; Hewitt, 2004; Petit et al., 2005). Species conservation is an urgent global issue, especially in biodiversity hotspots, such as the QTP and the surrounding mountains, which are the highest priority areas for conservation research and implementation (Myers et al., 2000; Beger et al., 2015).
Qinghai–Tibetan Plateau is the highest plateau in the world covering an area of about 2.5 million km2 (Sun et al., 2014). It is one of the most important grazing land as well as natural reserve in China, which has harsh attributes and diverse ecological conditions, such as drought, low temperatures, strong winds, and high radiation (Yan et al., 2009). The QTP began to uplift 50 million years ago (Eocene), especially from Miocene to Pliocene (Miocene, 23–2.6 Ma), which greatly changed the surface pattern of the region. The rapid and extensive uplift of the QTP was considered a major driving force in shaping such high species diversity, and can generate geographic barriers between populations and restrict gene flow (Mao et al., 2010; Xu et al., 2010). Several cycles of climate change that produced mountain glaciers seem to have occurred in the Holocene (Zheng et al., 2002; Wang et al., 2012). This change of topography and geomorphology drives the genetic differentiation of plant related species and intra-specific populations. Climate change during the Quaternary Ice Age had a significant impact on the genetic distribution of plants on the QTP (Hewitt, 2000). The uplift, climate change of QTP, and restricted gene flow will lead to divergence between populations and new speciation over time due to the enhanced effects of genetic drift and/or natural selection (Hoskin et al., 2005). Phylogenetic, phylogeographic, and ecological studies support plant diversifications on the QTP through multiple mechanisms such as allopatric speciation via geographic isolation, climatic oscillations and divergences, pollinator-mediated isolation, diploid hybridization and introgression, and allopolyploidy. Nevertheless, a few recent studies examined the evolution of unique endemic plants and they have found that several monotypic plant genera in the alpine regions of the QTP are nested within larger genera (Friesen et al., 2000; Tian et al., 2011; Liu et al., 2012; Nie et al., 2013; Wen et al., 2014). Unfortunately, almost 90% of alpine grasslands on QTP region have been degraded due to overgrazing, rodent activities, and climate change in the recent years (Dong et al., 2013; Li et al., 2013). These factors can severely affect the feasibility of protecting Elymus L. germplasm resources. Much work especially on the widely distributed species with dominant composition and decisive function of grassland community is needed toward understanding the evolutionary mechanisms of plant diversifications on the QTP ecosystem.
With its diverse natural landscape and complex geological and climatic history, the QTP is a hotspot for biodiversity and a fascinating natural laboratory for studying how species have evolved and adapted (Li et al., 2012). Elymus sibiricus L. is the typical species of the genus Elymus L., which is the largest genus of Triticeae with about 150 species in temperate regions of the northern hemisphere (Barkworth and Dewey, 1984). Elymus sibiricus is a perennial, self-fertilizing grass, and an allotetraploid with the StStHH genome constitution, which is widely distributed on QTP (Zhou et al., 2016). It usually grows in wet meadows, sandy beaches by rivers, between open forests, on sunny or semi-shaded slopes of mountains, or valleys from 1,000 to 4,000 m (Ma et al., 2012). In addition to its native ecological function, E. sibiricus has become one of the most widely cultivated and used leading grass species on the QTP area in the recent years because of its high yield, good quality, and excellent adaptability to cold and arid climate. However, the related studies of E. sibiricus were mostly focused on ecological protection and genetic diversity by molecular markers, such as SSR (Simple Sequence Repeats) and AFLP(Amplified Fragment Length Polymorphism) (Xie et al., 2015; Zhou et al., 2016). Despite an increasing number of studies, the genetic divergence and phylogenomic analyses of the dominant community species, E. sibiricus, on the QTP mountain grassland remains poorly understood, which leads to the lack of scientific and targeted protection and utilization of the most important germplasm resources for grazing land.
In the recent years, many molecular markers have been used to infer the phylogeography and biogeography of species to understand their evolution patterns (Behura, 2010). Technological advances have made high-throughput sequencing a viable new strategy for studying population genetics (Hohenlohe et al., 2010). Specific-locus amplified fragment sequencing (SLAF-seq) is an efficient method of large-scale genotyping for genetic study, which is based on RRL (reduced representation library) and high-throughput sequencing (Sun et al., 2013). At present, SLAF-seq technology has been widely used in system evolution, population genetics, and adaptive evolution related to natural selection (Geng et al., 2016; Zeng et al., 2017). Consequently, in this study, the SNPs (single nucleotide polymorphism) in 216 individuals from 44 populations were explored using SLAF technology to (1) discover the distribution pattern of genetic diversity, spatial-temporal population genetic differentiation, and history of demography of E. sibiricus on the QTP; (2) speculate the population historical dynamics of E. sibiricus, and explore the historical reasons for the formation of the existing distribution pattern and possible sanctuaries, to clarify the effects of geological events such as uplift and glacial cover of the QTP and human activities on the distribution pattern of plant species diversity in this area; (3) provide a historical perspective for the phylogeographical formation of alpine meadow population structure and dynamics; (4) discuss the conservation strategy for improving collection, evaluation, and breeding efficiency and pertinence of E. sibiricus germplasm resources and also important for genetic evolution theory of polyploid populations and biodiversity conservation on the QTP. Consequently, the phylogeographical study of the dominant community species in mountain systems has motivated growing scientific interest in this research field.
Materials and Methods
Sampling Locality of Elymus sibiricus
A total of 216 individuals from 44 natural populations of E. sibiricus were collected from Gansu, Tibet, Sichuan, and Qinghai Provinces representing almost the entire distribution range of the species on QTP (Figure 1). The 44-specimen information about the sample id, sample number, latitude, longitude, and altitude was listed in Supplementary Table 1. The 44 populations have different geographic distribution. Five individuals were sampled for each population, and with individuals at least 30 m apart. Fresh leaves from each individual were collected and dried immediately in silica gel (Wang et al., 1996), and these individuals were stored in the Yangzhou University research lab for the next DNA attraction.
Figure 1. Geographic location of E. sibiricus populations in this study (the population geographic distribution of best results in clusters by Admixture).
The DNA Extraction and Sequencing
A total of 216 individuals genomic DNA was extracted from silica-dried leaf tissue using the DP350 Plant DNA Kit (Tiangen Biotech, Beijing, China). The concentration and purity of the DNA were assessed with ultra-micro photometer. Then the extracted DNA was stored at –20°C for further use.
According to the genomic size and GC content, the E. sibiricus genome (has not been published yet) was selected as the reference genome for predicting restriction enzyme sites. Finally, the optimal enzyme RsaI and EcoRV-HF was selected. Then, all DNA samples were digested separately, and a 474–504 bp restriction enzyme fragment length was defined as a SLAF tag. Rice has the characteristics of small genome; easy manipulation; and dense molecular markers, and has commonalities with other grass crops in gene sequence; structure; sequence and function, so it is hailed as a model species for genome analysis by the international scientific community (Upadhyaya, 2007). To evaluate the efficiency and quality of enzyme digestion, the Oryza sativa ssp. Japonica genome was selected as control group. After digestion, a single nucleotide (A) was added to the 30°C end using dATP at 37°C, and then dual-index adapters were ligated to the A-tailed DNA fragments. A PCR amplification was subsequently performed using diluted restriction-ligation DNA as the template. The E.Z.N.A.® Cycle Pure Kit (Omega, London, United Kingdom) was used to purify the PCR products. The purified PCR products were pooled and incubated at 37°C with RsaI, T4 DNA ligase, ATP, and Solexa adapter. Then, the products were purified with a quick spin column (Qiagen, Hilden, Germany) and separated on a 2% agarose gel. All amplified fragments were diluted for the paired-end sequencing on an Illumina-HiSeq™ 2500 sequencing platform (Illumina, Inc., San Diago, CA, United States) at Beijing Biomarker Technologies Corporation (Beijing, China).
SNP Calling
The raw sequence data were identified by dual-index sequencing, and the quality of the specimen sequence was evaluated after the adapter of the sequence data were filtered. The quality of sequencing (Q = –10 × loge10) and GC content were used to evaluate the single-base error rate of sequencing quality. Based on the sequence similarity, the high-quality paired-end reads were clustered using the BLAT software (Kent, 2002). More than 90% of the sequence similarity between individuals was identified as a SLAF locus (Sun et al., 2013). A total of 6,747,462 SLAF tags were predicted which were evenly distributed in the genome. All sequencing reads were aligned to the reference genome (E. sibiricus genome) by BWA software (Li and Durbin, 2010). The genome analysis GATK and SAMtools were used for SNP calling, and the resulting data were combined into a SNP data set (Mckenna et al., 2010). Eventually, the SNPs with a minor allele frequency (MAF) more than 5% and integrity more than 50% were selected. The raw data were submitted to the National Center for Biotechnology Information (NCBI) under the BioProject PRJNA811778.
Population Genetic Analyses
A total of 30,159,737 SNPs was developed to construct the population structure and calculate genetic diversity. The measures of genetic diversity included the expected heterozygous number (He), Nei’s diversity index (H), observed allele number (Na), observed heterozygous number (Ho), Shannon’s wiener index (I), and polymorphism information content. These indexes can be calculated by POPGENE (Yeh et al., 1999). The spatial distribution map of genetic diversity was drawn by Arcmap software based on inverse distance weighting (IDW) (Mattioni et al., 2017). The regression analyses were conducted with genetic diversity indexes against altitude, longitude, and latitude to identify geographic patterns in the geographic pattern of genetic diversity. Analysis of molecular variance (AMOVA) was used to assess the degree of population differentiation. Additionally, pairwise fixation index (FST) between populations was also computed the genetic differentiation. Analysis of molecular variance and FST were calculated by Arlequin (Excoffier et al., 2005). The SNPs were used calculate genetic distance among populations by the Kimura two-parameter model and the pairwise deletion option parameters.
The Phylogenetic tree was constructed by MEGA X program (Sudhir et al., 2018) with the following parameters: Neighbor-joining method and 10,000 bootstrap replicates. Based on the selected SNPs, the population structure was analyzed by Admixture software (Alexander et al., 2009). When the number of populations (K) was ranging from 1 to 10, clustering is carried out. The cross-validation error rate (ΔK) was used to determine K (minimum value). A principal component analysis (PCA) whose result can be used to verify with the results of other clustering methods and was calculated by EIGENSOFT (Price et al., 2006). To verify whether isolation-by-distance (IBD) exists among populations, Mantel test was used to examine for correlation between genetic distance and geographical distance in GenAlEx (Peakall and Smouse, 2010).
The linkage disequilibrium (LD) decay rate is crucial for maintenance of disequilibrium in a population (Pradhan et al., 2020). Because LD may affect both PCA and Admixture, the marker set was pruned by excluding SNPs in strong LD using PLINK software (Purcell et al., 2007). The linkage disequilibrium (LD) level was assessed using pair-wise LD as r2-value (Pearson’s squared correlation coefficient) by PLINK software (Purcell et al., 2007) based on the formula proposed by Hill and Robertson (1968). The r2-value was calculated between all pairs of SNPs with inter-SNP distances of less than 10 kb (r2 and LD-window parameters).
Demography History
To infer the split and Admixture history among groups of 44 populations, the TreeMix (Pickrell et al., 2012) was used based on the SNPs. TreeMix, which simultaneously infers a tree of relationships and “migration events,” used allele frequency data to approximate an unrooted maximum likelihood (Beger et al., 2015) population tree. The effect of migration on the residual covariance matrix was examined by the stepwise likelihood method, and the optimal location of migration events in the population tree was determined. There are few reports of fossil data of E. sibiricus; therefore, based on the data of Elymus species and Aegilops tauschii1, the MCMCtree program implemented in the PAML v4.8 package (Yang, 2007) was used to estimate the divergence times among the species using a maximum likelihood method (Beger et al., 2015). The independent rates model (clock = 2), which assumes uncorrelated relaxed molecular clock. It was used to specify the prior of rates among the internal nodes, which followed a log-normal distribution. The three parameters (birth rate, l; death rate, m; and sampling fraction, r) in the birth–death process with species sampling were specified as 1, 1, and 0, respectively. A loose maximum bound for the root was set between 3 Ma and 9.4 Ma. The first 50,000 cycles in MCMCtree were discarded as burn-in, and every 50 cycles were sampled to obtain a total of 100,000 individuals.
To make sense of the population dynamic expansion history, the DNASP was used to calculate the mismatch distribution (Rozas et al., 2003). It is supposed that multimodal and rough distributions would characterize populations that have been stable for a long time, whereas populations that have experienced a sudden demographic expansion should display a unimodal and smooth distribution (Rogers and Harpending, 1992). Meanwhile, the neutrality test statistics was carried for population genetic data. Tajima’s D (Tajima, 1989), Fu and Li’s D, Fu and Li’s F (Li and Fu, 1993) statistics were used to estimate the population expansion and direct selection that the population underwent during the evolution process and the result can be complete by DNASP software (Rozas et al., 2003). These statistics are sensitive to the past population expansion and contraction. It is generally expected to be negative in the case of population expansion, and positive in the case of population contraction or the recent bottleneck events (Tajima, 1989; Li and Fu, 1993). To infer that E. sibiricus experienced bottleneck effect in the recent years, the SMC++ software was used to explore population size changes over time (Terhorst et al., 2016).
Analysis of Ecological Niche Model
The maximum entropy modeling (MaxEnt) (Phillips et al., 2006) and ArcGIS were used to construct the suitable area distribution of E. sibiricus during the Last Inter-Glacial period (LIG ∼120,000–140,000 years ago), Last Glacial Maximum (LGM ∼22,000 years ago), Mid-Holocene (6,000 years ago), present (∼1950–2000), and projected future (2070) (Phillips et al., 2006). The species distribution data included 88 sites come from field sampling and the Global Biodiversity Information Facility (GBIF)2. Environmental data (Supplementary Table 2) was downloaded from the WorldClim database3 and the map of QTP (1:400 million) from the National Geographic Center of China. The receiver operating characteristic (ROC) evaluated the accuracy of the model by the area under the curve (AUC) (Ritchie et al., 2017). The evaluation criteria of AUC (Greiner et al., 2000) was liated as follows: Less accurate (0.5 < AUC ≤ 0.7), moderately accurate (0.7 < AUC ≤ 0.9), highly accurate (0.9 < AUC < 1), and perfect tests (AUC = 1). In addition, the key environmental factors were selected by Jackknife method (Peterson and Cohoon, 1999). Finally, we used the natural discontinuity method to divide the suitable area of E. sibiricus into the following four grades: Not suitable, low suitable, moderately suitable, and highly suitable.
Results
Sequencing and SNP Calling
A total of 1,221.47 Mb reads were obtained from the Illumina HiSeq 2500 platform. The reads with a length of 125 bp were selected for the further evaluation and analysis to ensure the quality of the data. The average Q30 was 93.13% and the average GC content was 45.91%, which indicated that the reliability of the data. In addition, a total of 6,747,462 SLAF tags were explored, which 3,890,215 were polymorphic and the average depth was 5.10x. Finally, 30,159,737 SNPs were obtained in total (Supplementary Table 3).
Genetic Diversity
The average values of observed heterozygous (Ho) and expected heterozygous (He) were 0.4167 and 0.3750, respectively. The maximum and minimum values of Nei’s diversity index (H), Shannon–Wiener index (I), and polymorphism information content (PIC) appeared correspondingly in population QH09 and SC15, and these average values at species level were 0.4196, 0.5543, and 0.2971, respectively (Table 1). Additionally, the regression analysis showed that there was a significant decline in Ho with increasing longitude (Supplementary Figure 1A: R2 = 0.314, p < 0.001) and latitude (Supplementary Figure 1B: R2 = 0.158, p = 0.008), and a significant ascent in Ho with increasing altitude (Supplementary Figure 1C: R2 = 0.22, p = 0.001). However, there were no significant differences in other genetic diversity indexes except Ho. The spatial distribution of genetic diversity of E. sibiricus across the whole sampling range was shown as that genetic diversity values were lower in northeastern region and western Hengduan Mountains of QTP and higher in the center and southeast of QTP (Figure 2). However, from the whole QTP, multiple genetic diversity centers were obviously observed across all E. sibiricus populations.
Table 1. Intra-population genetic diversity of 44 E. sibiricus populations.
Figure 2. Genetic diversity map of 44 E. sibiricus populations: IDW interpolation of expected heterozygosity values (He).
Genetic Differentiation and Spatial Structure Among Populations
To explore the divergence among populations, in this study, 216 individuals of 44 E. sibiricus populations were divided into three clusters based on the NJ tree as indicated with different colors in Figure 3. Notably, PCA further proved the discovery (Figure 4). Based on the population allele frequency, the PC1 and PC2 accounted for 52.81 and 5.6% of total variation respectively. Consequently, it was reasonable for us to divide 44 populations into three groups (Figure 1). Group I was located in the northeastern QTP, which contained 13 populations. Group II was located in South-Tibet river basin, which included nine populations mainly from the XZ region. Group III was located in Hengduan Mountains, which covered a wider range of individuals including 22 populations of the QH, GS, XZ, and SC regions.
Figure 3. Phylogenetic tree constructed with polymorphic SNPs between three groups.
Figure 4. Principal components analysis of 44 E. sibiricus populations. Each circle represents an individual, and the circle color represents the population attribute of each individual.
Analyses of molecular variance of 44 E. sibiricus populations indicated that most of total variations mainly existed among groups and populations with the differentiation coefficient of 66.39 and 64.44%, respectively (Table 2). The AMOVA analysis against each group further was conducted. It is interestingly observed that the low inter-population genetic differentiation among populations within three groups, and the fixation index (FST) of the three groups were 0.12, 0.397, and 0.453, respectively. It indicated the notable local gene flow among the populations within each group. In addition, the FST value among three groups (Table 3) showed Groups I and II > Groups I and III > Groups II and III. The Mantel test showed a significant positive relationship between FST and pairwise geographic distance, indicating a significant pattern of isolation by distance across all 44 E. sibiricus populations (Supplementary Figure 2; r = 0.307, p < 0.001).
Table 2. Results of hierarchical AMOVA for 44 populations of E. sibiricus sampled within four regions.
Table 3. The FST value between three groups.
From the unrooted maximum likelihood tree inferred by TreeMix (m = 10, adding six migration events, Figure 5), the low weight of these proposed migration events was related to the low proportion of alleles in the offspring from the ancestral populations. This result may reflect the past gene flow. Notably, a small scale of gene flow occurred among E. sibiricus populations throughout the QTP. The TreeMix tree showed the gene flow mainly existed between Groups I and III. There were limited gene flow between Groups II and III, and no gene flow between Groups I and II. Additionally, the divergence time among 216 E. sibiricus individuals on the ML tree (Figure 6) showed that Group I diverged earliest compared with the other two groups, and the divergence time was 8.8 Ma between Groups I and II (With Aegilops tauschii as an outgroup, the divergence time between them was about 16.5 Ma). Whereas the late divergence between Groups II and III happened 5.92 Ma ago. Thus, it was probably to infer that the northeast of QTP was the origin area of E. sibiricus population 16.08 Ma ago.
Figure 5. Phylogenetic network of the genetic relationships among the E. sibiricus populations with different possible migration events (m = 10, six migration events) inferred from TreeMix analysis. The graph shows the topology and branch lengths according to the drift parameter. Migration arrows were colored according to their weight. The scale bar shows 10 times the average standard error (s.e.) of the entries in the sample covariance matrix W.
Figure 6. Phylogenetic tree constructed with polymorphic SNPs and divergence time scale between three groups.
The population structure was analyzed using the Admixture software to further understand the degree of Admixture in the populations. As shown in Figure 7A, the best k-value was 10. However, when k = 3, CV errors showed a significant decline, and CV errors value was already very low. So, the best k-value can be chosen as 2 or 3, indicating a strong genetic cluster among populations. When K = 2, a total of 216 individuals were divided into two clusters and Group I was separated to other groups. When K = 3, the obvious grouping (Groups I, II, and III) was observed, and the cluster result was consistent with PCA and phylogenetic relationship. The result indicated that E. sibiricus populations did not originate independently from one location on the QTP. All individuals with a low degree of Admixture were seen from all studied populations. From the dynamics of K = 2–10 (Figure 7B), Group I showed the evidence of some Admixture which may be due to gene flow and the recent expansion events. Groups II and III was stable and almost not admixed.
Figure 7. (A) Admixture estimation of the number of groups for K values ranging from 1 to 10. (B) Population structure analysis of the 216 individuals based on 30,159,737 SNPs. The bars in the x-axis indicate different individuals. Colors in each row represent structural components.
Linkage Disequilibrium
The decay of LD (r2) was calculated for the three groups and the maximum r2 values of the three groups were 0.5584, 0.7174, and 0.5646, respectively (Figure 8). The degree of LD decreased with the increasing marker distance. When the mean value of r2 declined to half from the initial value, the average distances of Groups I, II, and III were 10, 0.8121, and 0.9302 kb, respectively. The LD in Group I presented a less dramatic drop than in the two other groups, whereas LD extent was very similar between Groups II and III. This linkage disequilibrium decay pattern supported the results of the stronger bottleneck in Group I and relatively lighter selection pressure in Groups II and III.
Figure 8. Decay of average linkage disequilibrium (LD) over distance among SNPs in E. sibiricus populations.
Demography Statistics
To understand the demographic history of E. sibiricus, Tajima’s D, Fu, and Li’s D, Fu and Li’s F statistics were calculated to test for evidence of range expansion (Table 4). In all populations across the QTP, Fu and Li’s D, Fu and Li’s F-values were all significantly positive, indicating that the E. sibiricus possibly experience the obvious population contraction on the QTP. Likewise in Group I, significantly positive Fu and Li’s D, Fu and Li’s F-values suggested that these populations experienced the obvious contraction. In contrast, Tajima’s D, Fu and Li’s D, Fu and Li’s F-values for Groups II and III were all not significant, so it can infer that Groups II and III possibly did not experience the population expansion in the recent. On the other hand, the mismatch distribution for all populations showed a bimodal curve (Figure 9), which indicated the 44 populations of E. sibiricus on the QTP did not experience a large-scale population expansion. In addition, a bimodal curve of the mismatch distribution in Group I also indicated there no obvious population expansion. However, the unimodal Poisson distribution in Groups II and III suggested that they experienced one large-scale expansion recently. Therefore, the results of the comprehensive neutral test and the mismatch distribution suggested that Groups II and III possibly experienced populations expansion recently, but Group I as well as the whole E. sibiricus populations on the QTP likely experienced slight contraction. The results of demographic inference indicated over the last 10,000 years, the whole QTP E. sibiricus populations experienced decline trend (Figure 10). Groups I and II had similar population size between 10,000 and 1,000 years ago. Groups I and III had similar population size between 1,000 and 100 years ago. Starting from around 500 years ago, Group I experienced a decline which persisted until about 100 years ago, but about 70 years ago, Group I declined sharply. About 10,000 years ago, Group II experienced decline until about 600 years ago, the population size reached at 107 to those day. Starting from about 500 years ago, Group III experienced decline until about 100 years ago, and then experienced a rapid expansion and slight population fluctuations. Eventually, the population size reached at about 106. In addition, the populations of the three groups have expanded from 10,000 to 6,000 years ago (Figure 10), consistent with the MaxEnt events exactly at the period of LGM-Mid-Holocene transition.
Table 4. Neutrality tests of E. sibiricus populations based on SNPs data.
Figure 9. Mismatch distribution analysis for SNPs data of Group I (A), Group II (B), Group III (C), and overall (D).
Figure 10. Results of effective population size inference for three groups.
Potential Distribution
In this study, we modeled the LIG, LGM, Mid-Holocene, current, and the future distributions for E. sibiricus using MaxEnt version 3.4.1 as shown in Figure 11. MaxEnt model had high predictive ability and five periods AUC were acquired (LIG: AUC = 0.942, LGM: AUC = 0.933, Mid-Holocene: AUC = 0.931, current: AUC = 0.944, future: AUC = 0.938), which demonstrated the models had good model performances. According to the MaxEnt result, E. sibiricus distributed mainly at the northeast and southeast edge of QTP during the LIG period. During the LGM period, highly suitable area underwent a slight expansion, but the whole distribution persisted a slight contraction. The highly suitable area had a slight movement toward to the edge of QTP (Hengduan Mountains region). Populations exhibited expansion from the LGM to Mid-Holocene period. However, the highly suitable area had a slight contraction. Currently, potential suitable distribution further contracted in overall QTP range. Considering the changes of the potential distribution range in five periods, E. sibiricus had undergone the repeated historical process of contraction and expansion. Interestingly in the future, the species will potentially experience an obvious expansion. In addition, it was also found that the key environmental factors affecting the distribution of E. sibiricus were mean annual temperature (Bio1), mean temperature of coldest quarter (Bio11), annual precipitation (Bio12), and elevation (Elev) (Supplementary Figure 3).
Figure 11. Suitable distribution maps of E. sibiricus during different periods on Qinghai–Tibet Plateau (A) the Last Inter Glacial period (LIG ∼120,000–140,000 years ago); (B) Last Glacial Maximum (LGM ∼22,000 ya); (C) Mid-Holocene (6,000 years ago); (D) present (∼1950–2000); (E) projected future (2070); (F) changes of suitable areas at different levels in five periods.
Discussion
Genetic Diversity and Population Structure
Genetic diversity is the basis of species diversity and an important precursor to the study of any species, because the quantity and distribution of genetic diversity affect the evolution and reproductive potential of species or populations (Futuyma and Fairbairn, 1979). There is no obvious regular pattern of genetic diversity across all QTP E. sibiricus populations in this study. The QTP is characteristic of complex geographical environment, which may be the cause of variation in within-population genetic diversity. Mountain barriers may play a key role in species formation and diversity as their topographic complexity may lead to ecological stratification and environmental heterogeneity (Fjelds et al., 2012). The wide geographical distribution of E. sibiricus may also be one of the reasons for the difference in genetic diversity (Xie et al., 2015). The work of Hamrick and Godt (1996) has shown that geographic ranges were significantly correlated with plant variation at both species and population levels. The wide-ranging species maintain more variability than narrow ranging species. This finding will cause the difficult strategy for conservation and utilization of E. sibiricus germplasm, and prove QTP as its resourceful deposit. However, the spatial distribution of genetic diversity of E. sibiricus within the whole sampling range was shown that higher genetic diversity in the central and southeast of QTP (Figure 2). The result also found several populations had higher genetic diversity, speculating that these places may be the possible resourceful refuges. Shelter populations have a higher level of genetic diversity because they have a longer history of population dynamics as compared to post-glacial populations (Comes and Kadereit, 1998; Taberlet et al., 1998), which can be mirrored by the genetic diversity pattern in this study. In respect of the correlation between genetic diversity and altitude, the result of regression analysis shown that there was no significant correlation (Supplementary Figure 1). However, there were no significant differences in other genetic diversity indexes except Ho. Yan et al. (2009) indicated that a higher genetic diversity was observed in populations occurring at the medium altitudes (3,200–3,600 m) than those at the low and high altitudes for the two Elymus species (Elymus nutans and Elymus burchan-buddae) on the QTP. However, the results of our study were not completely consistent with their conclusion even in different Elymus species. The reasons for this difference were likely to be that the complex geological environment of a large-scale sampling area conceals the influence of altitude on genetic diversity, or the sampling sites higher than 3,800 m and lower than 3,000 m in this study are relatively scaring, and without covering the entire altitude range. Genetic structure analysis at the intra-specific level is important for future adaptive change or evolution (Schaal et al., 1991), which is informative for germplasm conservation. A PCA analysis, phylogenetic tree analysis, and Admixture analysis demonstrated that E. sibiricus showed a certain regional cluster (Figures 3, 4, 7). Populations from the same or similar geographical origin tend to cluster together, and 44 populations was divided into three genetic groups. These three groups were located in three separate locations with obvious geographical differences on the QTP. These findings will make it relatively easy for conservation and utilization of the wheatgrass germplasm on the QTP.
Genetic Differentiation Among Populations and Groups
According to AMOVA analysis of all populations, the results of this study were consistent with previous genetic studies of E. sibiricus, which found that most of the variations were distributed within populations or geographical areas (Ma et al., 2008; Yan et al., 2010). Genetic variation and gene flow among populations are important factors to explain population structure and dynamics history (Behura, 2010). Analyses of the genetic differentiation of E. sibiricus populations showed obvious evidence of major geographical partitions and barriers to historical gene flow. A significant negative relationship between FST and pairwise geographic distance across all 44 E. sibiricus populations (Supplementary Figure 2: r = 0.307, p < 0.001), indicating that the isolation-by-distance model seems to describe the pattern of differentiation among the populations in this study. There was a significant correlation between genetic distance and geographical distance of E. sibiricus populations due to the main self-pollination breeding system of E. sibiricus. In general, selfing species commonly have lower levels of genetic diversity and higher differentiation between populations than outcrossing plants (Schoen and Brown, 1991; Hamrick and Godt, 1996). In addition to its self-pollination system, some geographical isolation via barriers, such as rivers and mountains, can prevent E. sibiricus from spreading to a very large scale geographical level. Bockelmann et al. (2010) hold that within a 100 m habitat range, E. athericus populations with unique habitats have undergone significant genetic variation. When the distance is greater than 60 km, the impact of geographic distance isolation becomes weaker, which confirms the impact of habitat conditions on the genetic variation of Elymus genus plants. Genetic differentiation among intra-specific populations was promoted by geographical barriers of a series of high mountains on the QTP, while the formation of mountain systems, river diversions and “Sky island” restricted intra-specific gene flow, disrupting genetic patterns and speciation (Matuszak et al., 2016). Therefore, the self-pollination, geographic barriers, and habitat effect may have led to current general population structure of E. sibiricus on the QTP. In addition, a low genetic differentiation was observed among populations within each of the three groups, suggesting that local gene flow existed three groups. Local gene flow not only caused by self-pollination system but also by extra power, such as grazing activity (Yan et al., 2010). In addition, the study also found the gene flow in the whole region seemed to be restricted by high mountains on the map (Table 2 and Figure 5). The higher within-group variation of the three geographical groups (66.39%) may be attributed to isolation by environment of sampling at a larger spatial scale (Nybom, 2004).
The uplift of the mountain triggered the differentiation of species, changed the genetic structure of alpine plants, and affected the evolution of alpine plants (Liu et al., 2012; Wen et al., 2014). The reasons for this genetic structure were geographic isolation or adverse environmental changes (Avise et al., 1987). To display the fine spatial differentiation pattern of E. sibiricus, three genetic groups was divided by the PCA, Admixture analysis and phylogenetic tree. Based on the topography of the QTP, it was possible that the Tanggula Mountains and Tongtian rivers formed a geographical barrier between the two groups (Groups I and II). Group II was distributed in the South-Tibetan river basin, which had almost no gene flow compared with the other two groups. The degree of Groups II and III was low with the cause of the river obstacle, and there were no gene flow among Groups I and II with the geographical isolation of Tanggula Mountains. Interestingly, the gene flow mainly existed between Groups I and III because two groups were not isolated by very huge physical barriers (Figure 5). In addition, due to the complex topography of the QTP, the separation of the high mountains led to the isolation from populations of each other, hindering the long distance spread of seeds and pollen, further promoting genetic differentiation among populations. By the FST results inferred that the Tanggula and Bayan Har Mountains as the geographical barriers of the three groups, probably imposed divergence of them. Therefore, it can conclude that QTP specific mountainous environment will play an important role in shaping the genetic structure of E. sibiricus by preventing the gene flow. In summary, it is speculated that the reasons affecting the current phylogeographical pattern were as follows: (1) The environment changes due to the uplift of the QTP; (2) The geographic distance between the populations (Groups I and III are close in geographic distance, and gene flow are frequent); (3) Geographical barriers (the Tanggula and Bayan Har Mountains between Groups I and II).
Species Divergence and Demographic History
In this study, the divergence time of the 44 E. sibiricus populations showed Group I in the northeast of QTP first split compare with other groups, and the divergence time was 16.08 Ma. The divergence time between Group III and Group II was 5.92 Ma (late Miocene-early pliocene). The lasting climatic aridification and other events following the Himalayan Motion during 12∼6 Ma possibly played a direct role in the diversification of Group II and Group III. In addition, the abiotic event during 16∼5.92 Ma that led to the differentiation of these groups may be QTP uplift in 13∼7 Ma (Harrison et al., 1992; Willis and McElwain, 2014). During the uplift of the QTP and adjacent Himalayan mountains, long-term geological events generated great environmental differences, which might have triggered the diversification of species in the genus Elymus. Therefore, it can be speculated that the divergence time of the 44 E. sibiricus populations was significantly related to the uplift and climatic change of the QTP, which had given evidence that India–Asia collision and QTP uplift has progressed northeastward in a stepwise manner. Understanding LD patterns and mismatch distribution enhances the knowledge of the demographic processes within the population (Barría et al., 2013). The LD decay indicated Group I decay rate was low, which illustrated the strong natural selection. Natural selection has an important impact on shaping the genetic variation of a population, therefore promotes local adaptation (Kawecki and Ebert, 2004). Group I located in the northeast of the QTP, where human activities are frequent and the intensity of grassland utilization is high. It probably due to the increasing of human activities in the recent years (Group I is located in low altitude areas, and the country has made great efforts to develop tourism in the recent years, resulting in a substantial increase in the flow of people in this area). Thus, Group I was more selective, more adaptable, more diverse, and differentiated earlier. Group II and Group III possibly experienced expansion recently, but Group I likely experienced slight contraction. Group II experienced population contraction at 10 ka, the reason was the South-Tibetan river basin was affected by the Quaternary Glaciation. Regarding of the special geographical location with many land forms and micro-climates, it was be speculated that there had many small refuges of E. sibiricus on the QTP. Based on the analysis of genetic diversity and MaxEnt model results, theorizing that the QTP has multiple refuges in E. sibiricus, such as eastern Himalayas, Hengduan Mountains. Such area as the edge of the QTP (toward the northeastern, southeastern) might have an important role as an escape region.
The uplift of the plateau may lead to the change of habitat, and the adverse environment may lead to the contraction of plains in other areas (Liu et al., 2019). The uplift of the QTP has a great impact on China’s climate even the whole Asia (An, 2001), including the temperature decrease in some areas. Some cold-tolerant species with a wide range of habitats and vegetation zones certainly survived in multiple refugia on the QTP throughout glacial/inter-glacial periods. There were many Ice Age shelters on the edge of the QTP, which caused the alpine plants to retreat in the shelters during the Ice Age, and the population range of these plants began to expand during the inter-glacial or post-glacial warming. The MaxEnt model result of E. sibiricus also proved this historical process. The LGM is the period when the global temperature is the lowest, and it is also the coldest period recently. After the LGM, the temperature rose, the climate became warmer and humid, the suitable living area expanded, and the populations began to expand. These new climatic conditions may have been conducive to the expansion of populations of cold-resistant plants. During the glacial periods, most species experienced adverse climatic conditions in high-altitude mountain areas. Then contracted into refuges at low latitudes, and their ranges expanded again after the Ice Age, which leads to species divergence or secondary contact evolution (Hewitt, 2004; Nybom, 2004; Ohsawa and Ide, 2010). In the past century, the general trend of glacier change on the QTP retreated, and the marginal mountain area was more sensitive than the central area. The ecological niche modeling (ENM) results of grassland dominant species fit this pattern perfectly, partially mirroring the QTP historical change in grassland function. Elymus sibiricus most likely experienced slight population size change since the LGM. From the LIG to the LGM period, E. sibiricus exhibited significant range contraction. However, in the Mid-Holocene, the population range was expanded again. This is first report on the divergence time and demography of the most ecologically and economically important wheatgrass on the QTP high-cold grassland, which help us understand the grassland history.
Currently, the population distribution has decreased. The reason for this distribution pattern may be first that glaciers covered the QTP during the LGM, causing the populations retreat into the northeastern edge of the QTP and to the refuge. After the glacial period, the species of the Mid-Holocene re-expanded from the refuges. However, in the recent years, with the change of global climate and the increase of human activities, the suitable habitats for E. sibiricus were destroyed. It led to the decrease in E. sibiricus populations, especially in the northeast of QTP, but human activities were more intensive and the most significant climate warming. Additionally, according to the suitable distribution map, there was a refuge for E. sibiricus on the QTP during the glacial period, mainly distribute in the eastern Himalayas. As the climate warms, glaciers continue to retreat, accelerating the melting of permafrost. QTP climate warming causes regional and local climate change: soil exposure, severe desertification, grassland productivity decline. The future result of ENM showed the overall suitable area of E. sibiricus decreased compared to the Mid-Holocene, especially the moderate and high suitable areas. Therefore, the future of the QTP is not optimistic, it will keep warming trend, climate and ecological risk of disasters increased. It also funds that the distribution of E. sibiricus has shifted to the edge of the QTP. The change caused by climate warming, and would influence the distribution of E. sibiricus in the future. Combined with the genetic diversity and geographical distribution pattern of E. sibiricus, the regions with high genetic diversity were selected as the priority areas for protection and utilization, and scientific and accurate protection strategies were formulated. Thus, the reasonable and scientific preventions should be made to maintain the grassland function and protect the genetic germplasm of QTP E. sibiricus based on these evidences.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: National Center for Biotechnology Information (NCBI) BioProject database under accession number PRJNA811778.
Author Contributions
XY, YXG, and YZG designed the research. DL, CBZ, LY, LC, MY, WG, XL, JY, and SB provided individuals. MH and JZ performed the research and analyzed the data. MH wrote the manuscripts. XY, SS, and CJZ revised the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Foundation Science Commission of China (Nos. 31872995 and 31372358).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.882601/full#supplementary-material
Footnotes
References
Abbott, R. J. (2000). Molecular analysis of plant migration and refugia in the arctic. Science 289, 1343–1346. doi: 10.1126/science.289.5483.1343
Alexander, D. H., Novembre, J., and Lange, K. (2009). Fast model-based estimation of ancestry in unrelated individuals. Gen. Res. 19, 1655–1664. doi: 10.1101/gr.094052.109
An, Z. S. (2001). Evolution of asian monsoon and phased uplift of the himalaya-tibetan plateau since late miocene times. Nature 411:5035. doi: 10.1038/35075035
Avise, J. C., Arnold, J., Ball, R. M., Bermingham, E., Lamb, T., Neigel, J. E., et al. (1987). Intraspecific phylogeography:the mitochondrial DNA bridge between population genetics and systematics. Ann. Rev. Ecol. Syst. 18, 489–522. doi: 10.1016/j.ympev.2007.01.005
Barkworth, M. E., and Dewey, D. R. (1984). The genomic system of classification as a guide to intergeneric hybridization with the perennial triticeae. Gene Manipulat. Plant Breed. 35:202. doi: 10.2307/1221077
Barría, A., Christensen, K. A., Yoshida, G., Jedlicki, A., Leong, J. S., Rondeau, E. B., et al. (2013). Whole genome linkage disequilibrium and effective population size in a coho salmon (Oncorhynchus kisutch) breeding population using a high-density SNP array. Front. Genet. 10:498. doi: 10.3389/fgene.2019.00498
Beger, M., Mcgowan, J., Treml, E. A., Green, A. L., White, A. T., Wolff, N. H., et al. (2015). Integrating regional conservation priorities for multiple objectives into national policy. Nat. Commun. 6:8208. doi: 10.1038/ncomms9208
Behura, S. K. (2010). Molecular marker systems in insects: current trends and future avenues. Mol. Ecol. 15, 3087–3113. doi: 10.1111/j.1365-294X.2006.03014.x
Bockelmann, A. C., Reusch, T., Bijlsma, R., and Bakker, J. P. (2010). Habitat differentiation vs. isolation-by-distance: the genetic population structure of Elymus athericus in european salt marshes. Mol. Ecol. 12, 505–515. doi: 10.1046/j.1365-294x.2003.01706.x
Comes, H. P., and Kadereit, J. W. (1998). The effect of quaternary climatic changes on plant distribution and evolution. Trends Plant Sci. 3, 432–438.
Dong, Q. M., Zhao, X. Q., Wu, G. L., Shi, J. J., and Ren, G. H. (2013). A review of formation mechanism and restoration measures of “black-soil-type” degraded grassland in the qinghai-tibetan plateau. Environ. Earth Sci. 70, 2359–2370. doi: 10.1007/s12665-013-2338-7
Excoffier, L., Laval, G., and Schneider, S. (2005). Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evol. Bioinform. Online 1, 47–50. doi: 10.1177/117693430500100003
Fjelds, J., Bowie, R., and Rahbek, C. (2012). The role of mountain ranges in the diversification of birds. Ann. Rev. Ecol. Evol. Syst. 43, 249–265. doi: 10.1146/annurev-ecolsys-102710-145113
Friesen, N., Fritsch, R. M., Pollner, S., and Blattner, F. R. (2000). Molecular and morphological evidence for an origin of the aberrant genus milula within himalayan species of allium (alliacae). Mol. Phylogenet. Evol. 17, 209–218. doi: 10.1006/mpev.2000.0844
Geng, X., Jiang, C., Yang, J., Wang, L., and Wei, W. (2016). Rapid identification of candidate genes for seed weight using the SLAF-seq method in brassica napus. PLoS One 11:e0147580. doi: 10.1371/journal.pone.0147580
Greiner, M., Pfeiffer, D., and Smith, R. D. (2000). Principles and practical application of the receiver-operating characteristic analysis for diagnostic tests. Prevent. Vet. Med. 45, 23–41. doi: 10.1016/S0167-5877(00)00115-X
Guo, J., Liu, Y., Wang, Y., Chen, J., Li, Y., Huang, H., et al. (2012). Population structure of the wild soybean (Glycine soja) in China: implications from microsatellite analyses. Ann. Bot. 110, 777–785. doi: 10.1093/aob/mcs142
Hamrick, J. L., and Godt, M. (1996). Conservation genetics of endemic plant species. Conservat. Genet. 1995, 231–304.
Harrison, T. M., Copeland, P., Kidd, W., and Yin, A. N. (1992). Raising tibet. Science 255, 1663–1670. doi: 10.1126/science.255.5052.1663
Hewitt, G. (2000). The genetic legacy of the quaternary ice ages. Nature 405, 907–913. doi: 10.1038/35016000
Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the quaternary. Philos. Trans. R. Soc. B Biol. Sci. 359, 183–195.
Hill, W. G., and Robertson, A. (1968). Linkage disequilibrium in finite populations. Theor. Appl. Genet. 38, 226–231. doi: 10.1007/BF01245622
Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A., and Cresko, W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genet. 6:1000862. doi: 10.1371/journal.pgen.1000862
Hoskin, C. J., Higgie, M., McDonald, K. R., and Moritz, C. (2005). Reinforcement drives rapid allopatric speciation. Nature 437, 1353–1356. doi: 10.1038/nature04004
Kawecki, T. J., and Ebert, D. (2004). Conceptual issues in local adaptation. Ecol. Lett. 7:684. doi: 10.1111/j.1461-0248.2004.00684.x
Li, H., and Durbin, R. (2010). Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/btp324
Li, H., and Fu, X. (1993). Statistical tests of neutrality of mutations. Genetics 133, 693–709. doi: 10.1093/genetics/133.3.693
Li, X. L., Gao, J., Brierley, G., Qiao, Y. M., Zhang, J., and Yang, Y. W. (2013). Rangeland degradation on the qinghai-tibet plateau: implications for rehabilitation. land degradation. Development 24, 72–80. doi: 10.1002/ldr.1108
Li, Z., Zou, J., Mao, K., Lin, K., Li, H., Liu, J., et al. (2012). Population genetic evidence for complex evolutionary histories of four high altitude juniper species in the qinghai–tibetan plateau. Evolutio 66, 831–845. doi: 10.1111/j.1558-5646.2011.01466.x
Liu, J. Q., Sun, Y. S., Ge, X. J., Gao, L. M., and Qiu, Y. X. (2012). Phylogeographic studies of plants in China: advances in the past and directions in the future. J. Syst. Evol. 50, 267–275. doi: 10.1111/j.1759-6831.2012.00214.x
Liu, M. L., He, Y. L., López-Pujol, J., Jia, Y., and Li, Z. H. (2019). Complex population evolutionary history of four cold-tolerant Notopterygium herb species in the qinghai-tibetan plateau and adjacent areas. Heredity 123, 242–263. doi: 10.1038/s41437-019-0186-2
Ma, X., Chen, S., Zhang, X., Bai, S., and Zhang, C. (2012). Assessment of worldwide genetic diversity of siberian wild rye (Elymus sibiricus L.) germplasm based on gliadin analysis. Molecules 17:424. doi: 10.3390/molecules17044424
Ma, X., Zhang, X. Q., Zhou, Y. H., Bai, S. Q., and Liu, W. (2008). Assessing genetic diversity of Elymus sibiricus (poaceae: triticeae) populations from qinghai-tibet plateau by ISSR markers. Biochem. Syst. Ecol. 36, 514–522. doi: 10.1016/j.bse.2008.03.003
Mao, K., Hao, G., Liu, J., Adams, R. P., and Milne, R. I. (2010). Diversification and biogeography of juniperus (cupressaceae): variable diversification rates and multiple intercontinental dispersals. New Phytol. 188, 254–272. doi: 10.1111/j.1469-8137.2010.03351.x
Mattioni, C., Angela Martin, M., Chiocchini, F., Cherubini, M., Gaudet, M., Pollegioni, P., et al. (2017). Landscape genetics structure of european sweet chestnut (castanea sativa mill): indications for conservation priorities. Tree Genet. Gen. 13:39. doi: 10.1007/s11295-017-1123-2
Matuszak, S., Muellner-Riehl, A. N., Sun, H., and Favre, A. (2016). Dispersal routes between biodiversity hotspots in Asia: the case of the mountain genus tripterospermum (gentianinae, gentianaceae) and its close relatives. J. Biogeogr. 43, 580–590. doi: 10.1111/jbi.12617
Mckenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A. M., et al. (2010). The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Gen. Res. 20, 1297–1303. doi: 10.1101/gr.107524.110
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Fonseca, G., and Kent, J. (2000). Bio-diversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.1038/35002501
Nie, Z. L., Funk, V., Sun, H., Deng, T., Meng, Y., and Wen, J. (2013). Molecular phylo-geny of Anaphalis (asteraceae, gnaphalieae) with biogeographic implications in the Northern hemisphere. J. Plant Res. 126, 17–32. doi: 10.1007/s10265-012-0506-6
Nybom, H. (2004). Comparison of different nuclear DNA markers for estimating intraspecific genetic diversity in plants. Mol. Ecol. 13, 1143–1155. doi: 10.1111/j.1365-294X.2004.02141.x
Ohsawa, T., and Ide, Y. (2010). Global patterns of genetic variation in plant species along vertical and horizontal gradients on mountains. Global Ecol. Biogeogr. 17, 152–163.
Peakall, R., and Smouse, P. E. (2010). GENALEX 6: genetic analysis in excel. population genetic software for teaching and research. Mol. Ecol. Notes 6, 288–295. doi: 10.1111/j.1471-8286.2005.01155.x
Peterson, A. T., and Cohoon, K. P. (1999). Sensitivity of distributional prediction algorithms to geographic data completeness. Ecol. Mod. 117, 159–164. doi: 10.1016/S0304-3800(99)00023-X
Petit, R. J., Hampe, A., and Cheddadi, R. (2005). Climate changes and tree phylogeogr-aphy in the mediterranean. Taxon 54, 877–885. doi: 10.2307/25065568
Phillips, S. J., Anderson, R. P., and Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecol. Mod. 190, 231–259. doi: 10.1016/j.ecolmodel.2005.03.026
Pickrell, J. K., Pritchard, J. K., and Tang, H. (2012). Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8:e1002967. doi: 10.1371/journal.pgen.1002967
Pradhan, S. K., Pandit, E., Pawar, S., Naveenkumar, R., and Patnaik, S. (2020). Linkage disequilibrium mapping for grain Fe and Zn enhancing QTLs useful for nutrient dense rice breeding. BMC Plant Biol. 20:2262. doi: 10.1186/s12870-020-2262-4
Price, A. L., Patterson, N. J., Plenge, R. M., Weinblatt, M. E., Shadick, N. A., and Reich, D. (2006). Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet. 38, 904–909. doi: 10.1038/ng1847
Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M., Bender, D., et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795
Ritchie, E., Boyd, P., Lawson-Halasz, A., Hawari, J., Saucier, S., Scroggins, R., et al. (2017). The ecotoxicity of zinc and zinc-containing substances in soil with cons-ideration of metal-moiety approaches and organometal complexes. Environ. Toxicol. Chem. 36, 3324–3332. doi: 10.1002/etc.3910
Rogers, A. R., and Harpending, H. (1992). Population growth makes waves in the distr-ibution of pairwise genetic differences. Mol. Biol. Evol. 9, 552–569. doi: 10.1093/oxfordjournals.molbev.a040727
Rozas, J., SÁnchez-DelBarrio, J. C., Messeguer, X., and Rozas, R. (2003). DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19, 2496–2497. doi: 10.1079/9780851994758.0139
Schaal, B. A., Leverich, W. J., and Rogstad, S. H. (1991). Comparison of methods for assessing genetic variation in plant conservation biology. Genet. Conservat. Rare Plants 1991, 123–134.
Schoen, D. J., and Brown, A. H. (1991). Intraspecific variation in population gene diversity and effective population size correlates with the mating system in plants. Proc. Natl. Acad. Sci. U.S.A. 88, 4494–4497. doi: 10.1073/pnas.88.10.4494
Sudhir, K., Glen, S., Li, M., Christina, K., and Koichiro, T. (2018). MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms. Molecular Biology & Evolution. 35(6):1547–1549. doi: 10.1093/molbev/msy096
Sun, H., Yang, N., Chen, Y. S., Song, B., Liu, C. Q., Peng, D. L., et al. (2014). Survival and reproduction of plant species in the qinghai–tibet plateau. J. Syst. Evol. 52, 378–396. doi: 10.1111/jse.12092
Sun, X., Liu, D., Zhang, X., Li, W., Liu, H., Hong, W., et al. (2013). SLAF-seq: an efficient method of large-scale de novo snp discovery and genotyping using high-throughput sequencing. PLoS One 8:e58700. doi: 10.1371/journal.pone.0058700
Taberlet, P., Fumagalli, L., Wust-Saucy, A., and Cosson, J. (1998). Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7, 453–464. doi: 10.1046/j.1365-294x.1998.00289.x
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595. doi: 10.1093/genetics/123.3.585
Terhorst, J., Kamm, J. A., and Song, Y. S. (2016). Robust and scalable inference of population history from hundreds of unphased whole genomes. Nat. Genet. 49, 303–309. doi: 10.1038/ng.3748
Tian, B., Liu, R., Wang, L., Qiu, Q., Chen, K., and Liu, J. (2010). Phylogeographic analyses suggest that a deciduous species (ostryopsis davidiana decne., betulaceae) survived in northern China during the last glacial maximum. J. Biogeogr. 36, 2148–2155. doi: 10.1111/j.1365-2699.2009.02157.x
Tian, X., Luo, J., Wang, A., Mao, K., and Liu, J. (2011). On the origin of the woody buckwheat Fagopyrum tibeticum (=Parapteropyrum tibeticum) in the qinghai-tibetan plateau. Mol. Phylogenet. Evol. 61, 515–520. doi: 10.1016/j.ympev.2011.07.001
Upadhyaya, N. M. (2007). Rice Functional Genomics: Challenges, Progress And Prospects. New York, NY: Springer.
Wang, L., Wu, Z. Q., Bystriakova, N., Ansell, S. W., and Zhang, X. C. (2012). Phylogeography of the sino-himalayan fern lepisorus clathratus on “the roof of the world”. PLoS One 6:25896. doi: 10.1371/journal.pone.0025896
Wang, X. D., Wang, Z. P., and Zou, Y. P. (1996). An improved procedure for the isolation of nuclear DNA from leaves of wild grapevine dried with silica gel. Plant Mol. Biol. Rep. 14, 369–373. doi: 10.1007/BF02673369
Wen, J., Zhang, J. Q., Nie, Z. L., Zhong, Y., and Sun, H. (2014). Evolutionary diversifications of plants on the qinghai-tibetan plateau. Front. Genet. 5:4. doi: 10.3389/fgene.2014.00004
Xie, W., Zhao, X., Zhang, J., Wang, Y., and Liu, W. (2015). Assessment of genetic diversity of siberian wild rye (Elymus sibiricus L.) germplasms with variation of seed shattering and implication for future genetic improvement. Biochem. Syst. Ecol. 58, 211–218. doi: 10.1016/j.bse.2014.12.006
Xu, T., Abbott, R. J., Milne, R. I., Mao, K., Fang, K. D., Wu, G., et al. (2010). Phylogeography and allopatric divergence of cypress species (Cupressus L.) in the qinghai-tibetan plateau and adjacent regions. BMC Evol. Biol. 10:1–10. doi: 10.1186/1471-2148-10-194
Yan, J. J., Bai, S. Q., Chang, D., You, M. H., Zhang, C., and Li, D. X. (2010). Genetic diversity of wild Elymus sibiricus L. germplasm from qinghai-tibetan plateau in china detected by SSR markers. Acta Pratacult. Sin. 19, 173–183. doi: 10.3724/SP.J.1142.2010.40491
Yan, X. B., Guo, Y. X., Zhao, C., Liu, F. Y., and Lu, B. R. (2009). Intra-population genetic diversity of two wheatgrass species along altitude gradients on the qinghai-tibetan plateau: its implication for conservation and utilization. Conservat. Genet. 10:359. doi: 10.1007/s10592-008-9596-z
Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum Likelihood. Mol. Evol. 24, 1586–1591. doi: 10.1093/molbev/msm088
Zeng, B., Yan, H., Liu, X., Zang, W., Zhang, A., and Zhou, S. (2017). Genome-wide association study of rust traits in orchardgrass using SLAF-seq technology. Hereditas 154, 1–8. doi: 10.1186/s41065-017-0027-3
Zheng, B. X., Xu, Q. Q., and Shen, Y. P. (2002). The relationship between climate change and quaternary glacial cycles on the qinghai-tibetan plateau: review and speculation. Quat. Int. 97, 93–101.
Keywords: MaxEnt, diversity, differentiation, expansion, phylogeography, wheatgrass
Citation: Han M, Zhang J, Li D, Sun S, Zhang C, Zhang C, Yan L, Guan Y, Chen L, Guo Y, You M, Gou W, Li X, Yan J, Bai S and Yan X (2022) Phylogeographical Pattern and Population Evolution History of Indigenous Elymus sibiricus L. on Qinghai-Tibetan Plateau. Front. Plant Sci. 13:882601. doi: 10.3389/fpls.2022.882601
Received: 25 February 2022; Accepted: 06 June 2022;
Published: 29 June 2022.
Edited by:
Genlou Sun, Saint Mary’s University, CanadaReviewed by:
Xing Fan, Sichuan Agricultural University, ChinaWengang Xie, Lanzhou University, China
Copyright © 2022 Han, Zhang, Li, Sun, Zhang, Zhang, Yan, Guan, Chen, Guo, You, Gou, Li, Yan, Bai and Yan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xuebing Yan, yxbbjzz@163.com
†These authors have contributed equally to this work and share first authorship