In the present study, 361 silica-dried leaf samples from 35 populations of E. nutans were collected in the southeastern QTP from August to September 2019. In order to reveal the genetic characteristics of this species in response to environmental change, populations sampling in this study was carried out on a longitude gradient. For each population, we sampled 7 to 15 individuals, ensuring that the sampling locations were at least 5 meters apart. The seeds of these specimens were planted at the experimental field of Sichuan Academy of Grassland Science, Hongyuan, China (32°46.61’N, 102°32.63’E). Supplementary Table 1 provides detailed information regarding the collection sites, geographical coordinates, and elevations for each of the source populations. Genomic DNA extraction was carried out using the DP350 Plant DNA Kit (Tiangen Biotechnology, Beijing, China), following the manufacturer’s protocol. Subsequently, the quality and quantity of the extracted DNA samples were assessed using a NanoDrop-Lite instrument (Thermo Scientific, Waltham, MA, USA) and 1% agarose gel electrophoresis, respectively. Only the DNA samples meeting the required quality criteria were considered for further analysis. Subsequently, the qualified DNA samples were diluted to a concentration of 10 ng/μL for use.

Four cpDNA fragments (trnH-psbA, trnL-F, matK, and rbcL) were amplified and sequenced across all E. nutans samples. These fragments were amplified using primers previously reported in the literature, known for showcasing polymorphisms in other Elymus species (Xiong et al., 2022). Amplification was conducted in a 30 μL reaction volume, comprising of 30 ng of genomic DNA, 1.5 μL of each primer, and 15 μL of 2× Es Taq MasterMix (CoWin Biosciences, Beijing, China), using a C1000 Touch Thermal Cycler (BIO-RAD, Foster City, CA, USA). The polymerase chain reaction (PCR) conditions throughout the experiment were consistent with those described by Shaw et al (Shaw et al., 2005), and the resulting products were sequenced by Tsingke Biotech (Beijing, China). The genotypes of all 361 samples were determined using nine pairs of EST-SSR primers developed for E. nutans in our previous study based on transcriptomic data, and the PCR procedure followed the same protocol as we described previously (Li et al., 2023). The PCR products were subjected to capillary electrophoresis in an ABI 3730xl DNA analyzer (Applied Biosystems, Foster, CA, USA) and analyzed using GeneMarker v. 2.2 software (SoftGenetics, Pennsylvania, USA).

All cpDNA data were edited and adjusted using DNAMAN to obtain consensus sequences. The sequences were then aligned using the ClustalW algorithm in MEGA 6.0 (Tamura et al., 2013). PhyloSuite v1.2.3 was utilized to concatenate the sequences serially (Zhang et al., 2020). DnaSP v5.0 was employed to calculate haplotype (H), haplotype diversity (Hd), nucleotide diversity (π), and population genetic distance based on the concatenated sequences (Librado and Rozas, 2009). A haplotype network was constructed using the median-joining method in PopART software (Leigh and Bryant, 2015). The haplotype phylogenetic tree was constructed by the maximum likelihood (ML) method in MEGA 6.0 software, Elymus repens and Elymus ciliaris were chosen as outgroups (Bouckaert et al., 2014). The maximum clade credibility tree and additional summary statistics were visualized using FigTree 1.3.1. This enabled the assessment of potential spatial expansion within the populations, Tajima’s D and Fu’s Fs were calculated using DnaSP version 5.0. Mismatch distribution analysis and neutrality tests were employed (Tajima, 1989; Fu and Li, 1993). Analysis of molecular variance (AMOVA) was conducted using Arlequin v3.5 to estimate genetic variance within and between populations (Excoffier and Lischer, 2010).

In the analysis, the raw data matrix generated by EST-SSR markers was examined using the GenAlEx 6.5102 project (Peakall and Smouse, 2012). This encompassed parameters such as Shannon’s information index (I), observed number of alleles (Na), and the number of effective alleles (Ne). The genetic structure of E. nutans populations was investigated via the Bayesian clustering method in Structure v.2.3.4 (Falush et al., 2007). Determination of the optimal number of clusters (K) involved 20 independent runs for each K value ranging from 1 to 18. Structure Harvester v0.6.94 (Earl and Vonholdt, 2012) facilitated this analysis, each run comprising 500,000 Monte Carlo Markov Chain (MCMC) replicates. An admixture model was utilized, incorporating a burn-in period of 10,000 replicates. The ΔK method (Evanno et al., 2005) was employed to estimate the most likely number of clusters. Subsequently, the 20 replicates underwent clustering and permutation using CLUMPP v1.1 with the LargeK Greedy algorithm (Jakobsson and Rosenberg, 2007). To explore genetic differentiation at the population level, an analysis of molecular variance (AMOVA) was conducted within the GenAlEx 6.5102 project. Additionally, Barrier v2.2, based on Monmonier’s maximum difference algorithm (Manni et al., 2004), was employed to predict major genetic barriers’ geographical locations between populations.

In investigating the relationships between geographical and environmental factors and population genetic differentiation, population genetic distance indices for SSRs and cpDNA were computed using GenAlEx 6.5102 and DnaSP v5.0, respectively. Pairwise geographical distances were calculated utilizing the “geosphere” package in R, while environmental Euclidean distances were derived from nineteen bioclimatic variables obtained from the WorldClim website (https://www.worldclim.org/). Regression analyses were employed to evaluate isolation by distance (IBD) and isolation by environment (IBE) through the plotting of genetic distance against geographical and environmental Euclidean distances, respectively. These analyses facilitated an exploration of potential influences of geographical and environmental factors on population genetic differentiation.

To delve further into potential local adaptation due to environmental fluctuations, the correlation between climate variables and molecular markers was examined. Considering the high ploidy of E. nutans, variations in the detected allele numbers among different molecular markers across populations were anticipated. Assessing adaptive trends at specific sites involved establishing a linear regression relationship between the number of effective alleles and climate variables. This analysis offered insights into potential adaptive alterations in response to environmental factors at the identified loci.

The total length of the four cpDNA sequences was 4689 bp, and the lengths of the trnH-psbA, trnL-trnF, matK, and rbcL regions were 667, 1052, 1538, and 1432 bp, respectively. All sequences were deposited in GenBank under accession numbers: OR421574-OR423017. As cpDNA regions are uniparentally inherited markers, we used these four chloroplast fragments in our subsequent population genetics analysis. Among the 361 individuals sampled from 35 populations of E. nutans, a total of 20 variable sites were detected in the combined cpDNA. Genetic diversity indicators for various regions of E. nutans can be found in Table 1 . The DQ4 population exhibited the highest values for number of haplotypes (h = 6) and haplotype diversity (Hd = 0.911), followed by population DQ3 (h = 5, Hd = 0.833). The lowest values were observed in populations RQ1, JD2, and SN1, with a single haplotype each (h = 1, Hd = 0). On the other hand, the GG4 population showed the highest nucleotide diversity (π = 1.74 ×10^-3), followed by population DX2 (π = 1.06 ×10^-3), these comprehensive information regarding sample locations, sample sizes, and descriptive statistics of genetic variation for the populations is available in Supplementary Table 1 . Analyses of molecular variance (AMOVAs) based on cpDNA sequence data revealed that the variation among populations of E. nutans was equivalent to the variation within populations, accounting for 49.76% and 50.24%, respectively ( Table 2 ).

A total of 19 distinct haplotypes (H1–H19) were identified, and the spatial distribution of each population, along with its corresponding haplotype composition, is illustrated in Figure 1 . The differentiation among haplotype categories is visually discernible through color distinctions, as demonstrated in the haplotype network depicted in Figure 2 . Both the ML tree and the haplotype network consistently exhibit correspondence in Figure 2 . Based on the distribution patterns in the ML tree and haplotype network, the haplotypes can be systematically grouped into two categories. Group1 comprises 13 haplotypes, prominently featuring the most widely distributed and frequently occurring haplotype H4, alongside its subsequent haplotype H5. These two extensively distributed haplotypes hold central positions within the network diagram, suggesting their potential status as ancestral haplotypes. Additionally, Group1 encompasses 7 unique haplotypes, with H13, H14, H16, and H18 exclusive to Qamdo, while H11 exclusively appears in Nagqu. In contrast, Group2 consists of 6 haplotypes, none of which are observed in Nagqu. Notably, Group2 includes 3 haplotypes from Qamdo (H2, H15, H17) and 1 haplotype H8 from Nyingchi. To gain insight into the historical dynamics of E. nutans, mismatch distribution analysis revealed a multimodal pattern ( Figure 3 ), and neutrality tests showed that both Fu and Li’s F* (Fu and Li’s F* = 1.22427, p > 0.10) and Tajima’s D value (Tajima’s D = –0.07339, p>0.10) were not statistically significant. These results do not support population expansions in E. nutans.

Nine EST-SSR markers used to assess the diversity of E. nutans populations exhibited polymorphism ( Supplementary Table 2 ). The observed allele number (Na) detected by each primer ranged from 1.37 (EN91) to 3.57 (EN62), with an average of 2.65. The effective number of alleles (Ne) per locus ranged from 1.14 (EN91) to 3.30 (EN62), with a mean of 2.21. The average Shannon’s information index (I) was 0.46, varying from 0.15 (EN91) to 0.74 (EN5). The observed heterozygosity (Ho) and expected heterozygosity (He) varied across primers, with the lowest values observed in EN91 (0.06 and 0.11, respectively) and the highest in EN67 (0.99 and 0.55, respectively). The overall average Ho and He were 0.32 and 0.29, respectively. Table 1 provides an overview of the genetic variability at the population level. The Na and Ne per population ranged from 1.56 (MK2) to 3.56 (MG) and from 1.55 (MK2) to 2.82 (ZG3), with mean values of 2.65 and 2.21, respectively. Shannon’s information index (I) and expected heterozygosity (He) detected in each population ranged from MK2 (I = 0.12; He = 0.08) to ZG3 (I = 0.74; He = 0.46), with averages of 0.46 and 0.32, respectively.

The software Barrier was employed to identify genetic barriers among the 35 populations, represented by red lines in Figure 4A . The study confirmed that the Salween River serves as a significant genetic barrier in the region. Additionally, crisscrossing genetic barriers were observed near the Mekong River in the eastern part of the study area. The Bayesian cluster analysis (Structure) based on microsatellite data indicated the presence of three optimal clusters ( Supplementary Figure S1 ), suggesting that the 361 individuals likely belong to two main genetic clusters. Using CLUMPP to determine the most optimal of the 20 replicates, a plot of the structure of the 35 populations was constructed ( Figure 4B ). Despite a considerable degree of hybridization between populations, their genetic backgrounds can still be differentiated through genetic barriers. Furthermore, the AMOVA results demonstrated that a significant proportion (62%) of the genetic variation exists within the 35 E. nutans populations, while 38% of the genetic variation is found among populations ( Table 2 ).

Based on the Mantel test using pairwise distance values, our study revealed that geographic distance had a weak effect on the genetic differentiation detected by microsatellites ( Figure 5A , r=0.181; p=0.008), while environmental distance did not show a significant impact ( Figure 5B , r=0.064; p=0.128). However, there was no significant correlation between either geographic distance or environmental distance and the population genetic differentiation detected by cpDNA ( Figures 5C, D ). Nevertheless, we did observe a significant correlation between geographic distance and environmental distance among sampling points ( Figure 5E , r = 0.495; p = 0.001), as well as a significant correlation between genetic distances among populations detected by cpDNA and microsatellites ( Figure 5F , r = 0.209; p = 0.001).

The variation of alleles is likely attributed to adaptation to external environmental factors (Ramírez-Valiente et al., 2010). Here, we examined the correlation between the number of effective alleles (Ne) in each population and the environmental factors. Three of the markers detected a significant correlation between the number of effective alleles detected and ecogeographical data ( Figure 6 ). Specifically, the Ne detected by marker EN57 and EN80 exhibited a negative correlation with annual precipitation ( Figures 6A, B ), while the Ne detected by marker EN80 showed a significant positive correlation with the seasonal variability of precipitation ( Figure 6C ). Additionally, the Ne detected by marker EN5 displayed a significant positive correlation with the annual temperature range ( Figure 6D ). Upon annotating these loci with gene functions within the database, it was revealed that the marker EN57 locus is linked to glutathione S-transferase T1. In contrast, the marker EN5 and EN80 loci exhibited homology with hypothetical proteins in species such as Aegilops tauschii, Hordeum vulgare, and Setaria italica, among others. Notably, the precise functions of these hypothetical proteins remain uncharacterized and warrant further accurate annotation.

Genetic diversity, encompassing the total genetic variation among individuals within different populations of a species, serves as a crucial indicator of population adaptability to changing environments (Bell and Collins, 2008). Given the climate sensitivity of the QTP region, investigating the genetic diversity of species in this area holds significant importance (Wambulwa et al., 2021). E. nutans, a prominent herbaceous plant widely distributed in the QTP, has demonstrated exceptionally high levels of genetic variation through various approaches, including phenotype, ISSR, and SSR analyses (Chen et al., 2009; Li et al., 2023). However, prior studies on genetic variation primarily relied on individual or pooled genotyping resources, which might not sufficiently capture the spatial variation of genetic diversity among populations within this species.

To bridge the existing knowledge gap, a comprehensive analysis was undertaken, examining genetic variability among 361 specimens drawn from 35 distinct populations throughout the QTP region, with a specific focus on variations in nuclear and plastid fragments. The study unveiled substantial genetic diversity in both nuclear and plastid genes. Notably, the number of alleles (Na = 2.65) and effective alleles (Ne = 2.21) identified through EST-SSR markers were in close alignment with those reported in a prior investigation of the closely associated and sympatric species, E. breviaristatus (Li et al., 2022; Na = 2.71; Ne = 2.28). Nevertheless, it is imperative to acknowledge that the Shannon information index (I) recorded for the E. nutans population (I = 0.46) was significantly elevated in comparison to that of E. breviaristatus (I = 0.13), underscoring a distinctive and intriguing pattern of genetic diversity. As a species with an extensive distribution, E. nutans demonstrates a more homogeneous allele distribution, potentially attributed to its propensity to sustain larger population sizes and engage in more substantial gene flow across various populations (Magota et al., 2021). Conversely, the endangered status of E. breviaristatus may exert limitations on gene flow and genetic exchange, culminating in a diminished Shannon’s information index (I). This stark contrast elucidates the impact of species distribution and conservation status on genetic diversity and gene flow within these species.

A high Hd (0.805) and π (0.67×10^-3) were revealed in our study, indicating values lower than those observed in the related and widespread species of Elymus sibiricus in the QTP region detected based on cpDNA (Hd = 0.834, π = 1.08×10^-3). Conversely, these values surpassed those documented for E. sibiricus populations in Xinjiang and northern China, which exhibit respective Hd values of 0.527 and 0.543, and π values of 0.042×10^-3 and 0.53×10^-3, respectively (Xiong et al., 2022). Similarly, this high degree of genetic diversity has also been confirmed in other plants on the Tibetan Plateau, such as Stellera chamaejasme (Hd = 0.834) (Zhang et al., 2010), Iris lozyi (Hd = 0.820) (Zhang et al., 2021). This observation corroborates the notion that the QTP serves as a focal point for the differentiation in most species, including Elymus species. Further, on a population scale, specific groups displaying elevated haplotype counts and notable haplotype diversity, as exemplified by populations such as DQ3, DQ4, RQ2, and others, are strategically positioned within the geographical span delimited by the Salween River and the Mekong River. These observations substantiate the hypothesis that this geographic region holds substantial significance as a sanctuary during ice ages for E. nutans. Furthermore, the locale has historically held a pivotal role as a biogeographic demarcation within the East Asian floral context, commonly recognized as the Ward Line-Mekong-Salween Divide (MSD) (Ward, 1921; Luo et al., 2017), which is highly suggestive of its potential as a hub for E. nutans diversification. In addition, this hypothesis is supported by the previous findings of Yu et al (Yu et al., 2019), who identified this field as one of the nine evolutionary hotspots in the QTP, bolsters the proposition of its paramount significance in the evolutionary trajectory of E. nutans, namely the eastern part of Nyenchen Tanglha Mountains.

The analysis of cpDNA haplotypes and microsatellite structures in this study suggests that the genetic differentiation among the 35 populations of E. nutans lacks distinct definition. This conclusion is supported by the presence of individuals displaying mixed genetic backgrounds within the majority of populations. Moreover, the AMOVA results emphasize that genetic variation primarily arises within populations, with comparatively lower levels of genetic divergence observed among populations. We attribute this observed divergence to E. nutans’ wider ecological niche and larger effective population size, which likely enhance individuals’ migration and dispersion across diverse habitats, facilitating gene flow among populations. Such occurrences are common among widely distributed plant species and are pivotal in maintaining genetic diversity and adaptability within populations (Broadhurst et al., 2018; Kahl et al., 2021; Veto et al., 2023).

Nonetheless, there still exist pronounced genetic barriers among certain adjacent populations. These barriers could be attributed to specific factors such as geographic isolation, environmental differences, or ecological niche divergence (Ryan et al., 2017; Yu et al., 2021), which collectively restrict gene flow among populations. As indicated by the results obtained from the Barrier analysis ( Figure 3 ), our study reveals that the genetic barriers detected are primarily aligned with the Salween River and the Mekong River, impeding gene flow between populations situated on opposite sides of these rivers. The Ward Line-MSD, formed by these two rivers, holds significance as a major biogeographic boundary within the East Asian plant region (Ward, 1921). It has been extensively studied and confirmed to exist in various species, such as Marmoritis complanatum, Koenigia forrestii, and Sinopodophyllum hexandrum (Li et al., 2011; Luo et al., 2017; Rana et al., 2023). Our research, focusing on the intraspecific genetic differentiation of populations and operating at a finer scale, further underscores the importance of MSD in the dynamic process of species differentiation. The elucidation of these intricacies significantly enriches our comprehension of species evolution and ecological dynamics. Such an enhanced understanding, in turn, facilitates the implementation of more refined and effective strategies for the collection, conservation, and management of the germplasm resources pertaining to this species.

Although our study highlights the extensive effect of the MSD concerning the genetic separation of E. nutans populations, which restricts gene flow and leads to prompt genetic divergence, a noteworthy phenomenon presents itself in the form of a weak yet statistically significant pattern of isolation by distance (IBD) among populations only demonstrated via microsatellite data. In contrast to the robust IBD pattern uncovered in our earlier study of E. breviaristatus populations (Li et al., 2022), the relatively weaker or even absent IBD pattern discerned in E. nutans populations further reinforces the notion of heightened dispersal capabilities within these populations. However, the insights gleaned from the cpDNA analysis do not corroborate this phenomenon. It is posited that this disparity may originate from the fundamental differences in transmission modes, rates of genetic drift, and migration patterns inherent to nuclear and plastid DNA, culminating in the divergent genetic configurations manifested in these two distinct genetic substrates (Soltis and Kuzoff, 1995).

Perennial plants strategically accumulate diverse allelic variants in response to various environmental conditions, reflecting their adaptive mechanism to external selection pressures (Castillo et al., 2010; Sork et al., 2010; Raschke et al., 2015). In our investigation, we identified significant divergence within three analyzed loci among populations situated in regions marked by more extreme climatic conditions, such as low precipitation or imbalances in precipitation and temperature. Of particular significance, notable correlations were observed between the number of effective alleles and specific environmental factors, indicating a discernible influence of natural selection on these genetic markers, thereby facilitating localized genetic differentiation. While two of these loci lack prior annotations, it is noteworthy that marker EN67 is unequivocally associated with a glutathione S-transferase T1 in wheat and Arabidopsis. Glutathione S-transferase, an essential component of the glutathione antioxidant system, plays a pivotal role in managing oxidative stress and detoxifying harmful compounds within plants (Sappl et al., 2009; Nianiou-Obeidat et al., 2017). We posit that this functional attribute may contribute to allelic divergence observed in arid regions. Additionally, although the marker EN5 and EN80 motifs lack established annotations, we hypothesize that they are linked to mechanisms enabling adaptation to climatic extremes. These motifs hold promise as potential candidate genes for genetic breeding in wheat plants, particularly to enhance resilience against environmental challenges. Further research is warranted to unravel the precise functional implications of these loci and their roles in plant adaptation.

Collectively, our findings, coupled with our prior results, implies that concerning gene loci affected by environmental selection, the lack of a clear-cut isolation by environment (IBE) pattern signifies a more intricate association between genetic differentiation and environmental adaptation than previously presumed. This complexity likely arises from a confluence of interacting factors, highlighting the multifaceted nature of the evolutionary processes governing population differentiation and adaptation.