Population genetic structure and evolutionary history of Psammochloa villosa (Trin.) Bor (Poaceae) revealed by AFLP marker

Abstract Psammochloa villosa is an ecologically important desert grass that occurs in the Inner Mongolian Plateau where it is frequently the dominant species and is involved in sand stabilization and wind breaking. We sought to generate a preliminary demographic framework for P. villosa to support the future studies of this species, its conservation, and sustainable utilization. To accomplish this, we characterized the genetic diversity and structure of 210 individuals from 43 natural populations of P. villosa using amplified fragment length polymorphism (AFLP) markers. We obtained 1,728 well‐defined amplified bands from eight pairs of primers, of which 1,654 bands (95.7%) were polymorphic. Results obtained from the AFLPs suggested effective alleles among populations of 1.32, a Nei's standard genetic distance value of 0.206, a Shannon index of 0.332, a coefficient of gene differentiation (G ST) of 0.469, and a gene flow parameter (Nm) of 0.576. All these values indicate that there is abundant genetic diversity in P. villosa, but limited gene flow. An analysis of molecular variance (AMOVA) showed that genetic variation mainly exists within populations (64.2%), and we found that the most genetically similar populations were often not geographically adjacent. Thus, this suggests that the mechanisms of gene flow are surprisingly complex in this species and may occur over long distances. In addition, we predicted the distribution dynamics of P. villosa based on the spatial distribution modeling and found that its range has contracted continuously since the last interglacial period. We speculate that dry, cold climates have been critical in determining the geographic distribution of P. villosa during the Quaternary period. Our study provides new insights into the population genetics and evolutionary history of P. villosa in the Inner Mongolian Plateau and provides a resource that can be used to design in situ conservation actions and prioritize sustainable utilization.


| INTRODUC TI ON
The Quaternary period, comprising the Holocene and Pleistocene Epochs, spanned the last ~2.6 million years (Myr) and has been characterized by distinct climatic oscillations, especially alternating glacial and interglacial cycles in the Northern Hemisphere (Elias, 2013).
The glacial cycles covaried with, and probably profoundly affected, other aspects of the climate, including the intensity of the Asian monsoon, even in unglaciated regions Liu et al., 2018). Climate fluctuations during the Quaternary glaciations led to dramatic changes in the geographic distribution, genetic structure, and population demography of plant species (Jia et al., 2012;Liu et al., 2018).
Quaternary climate change is known to have strongly affected the distributions of plants in Northern China .
Northern China was dominated by deserts, which likely developed during the Quaternary based on geological interpretations of deposits of silt-sized sediment (i.e., loess) (Sun et al., 2018). The process of desert formation in Northern China was once understood to be a result of sustained orogenesis of the Qinghai-Tibetan Plateau and surrounding areas during the Quaternary (Meng & Zhang, 2011;Wu, 2002). However, new data showed that the plateau reached its height earlier than originally thought in the Miocene (Hu et al., 2020;Staisch et al., 2020) and suggested that the deserts might have originated during glacial periods, when global ice volume was high and, thus, liquid water was more limited (Sun et al., 2018). The same climatic processes that gave rise to the deserts also appear to have profoundly shaped regional plant diversity and yielded highly complex demographic histories of native species (Ge et al., 2005). In particular, many plants within the deserts of Northern China underwent an adaptive and demographic change in response to cold, arid conditions, such that the intermittent glacial periods may be the primary mechanism explaining modern distributions (e.g., El-Tayeh et al., 2020;Su & Zhang, 2013;Xu & Zhang, 2015). For example, Su and Zhang (2013) proposed that the onset of aridity during Quaternary glacial periods was a primary driver of population processes and structures in Nitraria sphaerocarpa Maxim. (Nitrariaceae), and Xu and Zhang (2015) revealed that periods of cold, arid conditions during the Pleistocene glaciations resulted in genetic differentiation and demographic structuring in Atraphaxis frutescens (L.) K. Koch (Polygonaceae).
Nevertheless, these prior studies on population histories of desert plants of Northern China have focused primarily on woody species, and studies on herbs of the region are largely lacking. To our knowledge, the only such study on an herb is on Delphinium naviculare W. T. Wang (Ranunculaceae) , which is endemic at midelevations within the Tianshan Mountains of Xinjiang Province of China (Wang & Warnock, 2001). Herbaceous plants of the deserts of Northern China merit further study because they comprise a vital component of desert plant communities, and studies focusing on dominant grass species are especially warranted (Meng & Zhang, 2011). Herbaceous plants may have been more sensitive to Quaternary climatic oscillations because they differ markedly from woody species in their responses to cold, often via the death of aboveground biomass as part of either an annual or perennial life cycle.
Moreover, dominant species likely achieved their present abundances due to their responses to the glacial cycles.
In modern times, deserts and semideserts, such as in Northern China, are extremely fragile ecosystems, the stability of which impacts global environmental conditions and influences climate change (Su, 2013). Although deserts typically have sparse vegetation, plants are critical for maintaining their integrity. In Northern China, the desert grassland ecosystems in particular are becoming rapidly degraded due to long-term overgrazing and desertification and, simultaneously, the desert is encroaching on arable land within the region (Deng et al., 2014;Li et al., 2004;Wei et al., 2020). These desert grasslands represent a large area within China and adjacent countries and occur at both low and high elevations. Dominant plant species within the grasslands are often psammophytes, which have special adaptations to resist being buried by sand and to tolerate having periodically exposed roots. At the same time, these plants help to anchor sands in place and prevent wind erosion. Therefore, they are critical for promoting environmental stability within the desert grasslands and preventing desert encroachment (Pan, 2006;Zhou et al., 2011).
In this study, we focused on one psammophyte, Psammochloa villosa (Trin.) Bor, which was treated in a monotypic genus of tribe Stipeae in Poaceae. This species, commonly called sand whip, is a perennial rhizomatous herb that is primarily distributed in the desert grasslands of northwestern China, especially in the Inner Mongolian Plateau, the Hexi Corridor of Gansu, central and northern Ningxia, and Northern Shanxi (Ma, 1994). It also occurs in the Gobi Desert of Mongolia, especially in Ömnögovǐ and Bayankhongor Provinces in the south and southwest of the country (Hilbig, 1995). Psammochloa villosa is ecologically widespread at low and high elevations (900-2,900 m) (Wu & Phillips, 2006). Its flowering and fruiting period is from September to November, and the seeds are 5-7 mm long with an average weight of 5.507 ± 0.053 mg (mean ± standard error; Huang, 2003). These lightweight seeds are potentially dispersed by the high winds that occur throughout much of its natural desert habitat. Nevertheless, seedlings are rarely observed (Zhu et al., 2005).
Psammochloa villosa is known to have high resistance to drought, cold, alkaline soils, disease, wind, and burial by sand, all of which likely represent evolutionary adaptations that facilitate its survival in grassland and dune areas (Lu, 1987;Wu & Phillips, 2006).
Previous research on P. villosa has been mainly focused on studying K E Y W O R D S desert grasslands, ecological niche modeling, Inner Mongolian Plateau, population genetics, SAMOVA its anatomy, embryology, and microbiology (e.g., Huang et al., 2004;Lv et al., 2018;Wang et al., 2011), with only a few focused on molecular markers (e.g., Li & Ge, 2001).
In this study, we investigated the influence of aridification and climatic oscillations on the genetic structure and evolutionary processes of P. villosa during the Quaternary in northwestern China using a population genetics approach based on amplified fragment length polymorphism (AFLP) combined with ecological niche modeling (ENM) to compare past, present, and future environmentally suitable habitats for the species. We used AFLPs, which were multilocus markers, and their mode of inheritance was dominant, because they remained extremely efficient for investigating genetic diversity, genetic structure, and population demography due to their high levels of polymorphism, their reproducible, reliable results that were unaffected by the developmental stage of plant materials, and their universality among plant species (Wang et al., 2008).
In addition, they have been used to resolve genetic structures and population demography in many diverse grass species such as Oryza sativa, Leymus racemosus, Orinus thoroldii, and O. kokonoricus (Cai et al., 2017;Liu et al., 2019;Zhang & Jia, 2002 We believe that, taken together, our results can provide a scientific basis for improved protection and sustainable utilization (e.g., as forage) of P. villosa within the fragile desert grassland ecosystems where the species occurs.

| Population sampling
We randomly sampled five to ten individuals from 43 populations of P. villosa in the field throughout its natural range in China and obtained a total of 210 individuals (Table S1 & Figure S1). We sampled individuals spaced at least 20 m apart in order to avoid sampling a single clone more than once. In the field, we immediately put the fresh leaves into sealed bags filled with silica gel and then stored them in the laboratory in a −20℃ freezer until processing. At each population collection locality, we obtained geocoordinates using a GPS measuring instrument (Garmin eTrex201x). From each population, we obtained and deposited one representative voucher specimen in the herbarium of Qinghai-Tibet Plateau Museum of Biology (QTPMB), Northwest Institute of Plateau Biology, Chinese Academy of Sciences, China.

| DNA extraction and AFLP scoring
We extracted total DNA from each sample according to a modification of the CTAB procedure (Doyle & Doyle, 1987) and accessed DNA quality using 1.0% agarose gel electrophoresis and the A260/ cally were included in the analysis (Rocha et al., 2015). In assessing the recovered fragments, we did not account for polyploidy because P. villosa was known to be diploid (2n = 40; Li et al., 2012). In total, we assessed 1,728 AFLP markers for the 210 individuals, and all interpretations were performed randomly. Besides, we estimated the error rate with the ratio between the observed number of phenotypic differences and the total number of phenotypic comparisons in order to track and assess genotyping errors (Bonin et al., 2010).

| Genetic diversity and population genetic structure
For each population, we evaluated genetic diversity and population genetic structure according to standard metrics GenAlEx 6.5 (Peakall & Smouse, 2012) and AFLP-SURV v1.0 under the assumption of Hardy-Weinberg equilibrium (HWE, Vekemans et al., 2002). These metrics included the number of individuals (N), percentage of polymorphic loci (PPL), observed number of alleles (Na), effective number of alleles (Ne), Shannon's information index (I; Lewontin, 1972), Nei's genetic diversity (h), expected heterozygosity (He), Nei's standard genetic distance (GD), total population diversity (Ht), genetic diversity within populations (Hs), genetic diversity between populations (Hb), and the population differentiation (F ST ). Meanwhile, we inferred the correlation between I and h, and I and He using spearman ranking in R 4.04 (http://www.r-proje ct.org/), and calculated the degree of genetic differentiation between populations (G ST ) as (Ht − Hs)/Ht (Nei, 1973), the parameter of gene exchange as Nm = 0.5(1 − G ST )/G ST (McDermott & McDonald, 1993), genetic diversity, coefficients of gene differentiation, and gene flow for eight pairs of AFLP primers in POPGENE 1.32 (Yeh et al., 1999). In addition, to eliminate the influence of codominance from AFLP molecular markers, we estimated the population differentiation (θ B ) using the Bayesian method in HICKORY v1.1 (Holsinger & Lewis, 2007;Holsinger et al., 2002), whose advantages of this method were that it did not assume HWE. We performed the Full, f = 0, and θ B = 0 models with default parameters (burn-in = 5,000, sample = 25,000, thin = 5) and determined the most suitable model based on the deviation information criterion (DIC) (Holsinger & Wallace, 2004).
Identifying the genetic structure of species or clusters of genetically associated populations could facilitate the detection of finerscale geographic structures (i.e., within groups) juxtaposed with broader, regional patterns (i.e., between groups) (Li et al., 2020). To assess structuring among the 43 populations of P. villosa, we generated a UPGMA tree from the genetic distance matrix derived from the binary AFLP dataset and performed bootstrapping of this tree using a custom R script (Supplemental File 1). We used the UPGMA tree to identify strongly supported clusters of populations.
As complementary to the UPGMA approach, we also constructed a similarity-based network in SplitsTree 4.13 (Huson & Bryant, 2006) to infer the relationships between individuals and populations by applying the Neighbor-Net algorithm with Jaccard's measure of distance. We further examined clusters of populations using a principal coordinate analysis (PCoA), from which we determined the optimal number of clusters by calculating the gap statistic (Tibshirani & Hastie, 2001) for axes one and two. The gap statistic represents a mathematically tractable method compared to the examination of a scree plot (Thorndike, 1953). In addition, we evaluated clusters of populations using SAMOVA, which took geographic adjacency into account. Within SAMOVA, we used a K-means method to select the best clustering scheme based on genetic variation coefficients (F CT ) (Li et al., 2020). For possible numbers of clusters, K, in the range two to ten, we performed 100 heuristic searches with 10,000 steps each, and we selected the value of K that minimized within-cluster F CT without over-partitioning. Subsequently, we inferred clusters of populations of P. villosa using STRUCTURE V2.2 (Hubisz et al., 2009), which differed from SAMOVA by not requiring that groupings be geographically adjacent. In STRUCTURE, we performed the analyses using an admixture model with independent allele frequencies for 90 independent runs for the number of clusters (K) ranging from one to ten. We applied 2 × 10 5 repetitions of the Markov chain Monte Carlo with a burn-in of 25%. To determine the best value of K for the STRUCTURE analyses, we used the ΔK statistical method (Evanno et al., 2005).
Based on all assessments on the optimal number of clusters of populations and the constituent populations of the optimal grouping, we further evaluated the genetic variation between populations within groups and between groups in SAMOVA 1.0 via an analysis of molecular variance (AMOVA, Excoffier et al., 1992) in ARLEQUIN v3.01 (Excoffier et al., 2005). Moreover, we performed the tests of neutrality with Tajima's D and Fu's Fs in ARLEQUIN v3.01 (Excoffier et al., 2005) and determined the correlation between F ST inferred from the binary matrix of scored AFLPs and geographic distance of the populations via a Mantel test (Mantel, 1967) in GenAlEx 6.5 (Peakall & Smouse, 2012) with 9,999 permutations to evaluate significance.

| Distribution modeling of P. villosa
In order to predict the effects of Quaternary climatic oscillations on the geographic distributions of P. villosa, we used ecological niche modeling (ENM) to compare the potential distributions of P. villosa at the Last Inter-Glacial (LIG, ~120,000-140,000 years before present), the Last Glacial Maximum (LGM, ~21,000 years before present), and the present. As input for the ENMs, we used the geocoordinates of the 43 populations of P. villosa that we sampled for this study as well as geocoordinates of specimens recorded in the Chinese Virtual Herbarium (CVH, http://www.cvh.ac.cn), the Global Biodiversity Information Facility (http://www.gbif.org), the China National Specimen Information Infrastructure (http://www.nsii.org. cn), and the Specimen Resources Sharing Platform for Education (http://mnh.scu.edu.cn/main.aspx). In total, after removing duplicate and ambiguous records, we obtained 155 georeferenced data points for using in the ENMs (Table S2).
Initially, to perform ENMs, we obtained 19 bioclimatic variables and three geographic factors (altitude, slope, and aspect) as global information system (GIS) layers from the WorldClim database (Hijmans et al., 2005, www.world clim.org) at 2.5 arc-min resolution.
This represented 22 total environmental variables for modeling. We followed Peterson and Nakazawa (2008) and applied the Spearman's correlation test to exclude highly correlated variables, a correlation of <0.75 compared to other variables. We performed preliminary modeling in MaxEnt 3.3.3k (Phillips & Dudík, 2008) with 75% of localities randomly selected for training and 25% selected for testing 500 times independently to ensure reliable results. We used these models to assess the relative contributions of the 22 environmental variables and excluded those exhibiting a relative contribution score ≥0.8.
Based on the outcome of Spearman's and the preliminary ENMs, we retained 11 variables to generate the final models. We preformed modeling using the same procedure as above with these 11 variables and the occurrence data and subsequently evaluated model performance using the area under the curve (AUC) of the receiver operating characteristic (ROC). The value of AUC ranges between 0 (randomness) and 1 (exact match), and values above 0.9 are generally regarded as indicating good model performance (Swets, 1988).
In order to obtain the geographic distribution, we projected the We used the ENMs to assess niche similarity between populations occurring in the clusters, or groups, that we inferred. To accomplish this, we calculated Schoener's D (Schoener, 1968) and standardized Hellinger distance (calculated as I) in ENMTools 1.3 (Warren et al., 2010). We obtained the null distribution of niche models for the identity test based on 1,000 pseudo-replicates generated by random sampling from the data points pooled for each pair of clusters. We determined D and I by comparing with null distributions drawn from pooled occurrences, retaining the original cluster sizes, and we generated histograms of frequency distributions in R 4.04.

| Conservation assessment of P. villosa
We performed a conservation assessment for P. villosa using the extent of occurrence (EOO) (Moat, 2007;Velzen & Wieringa, 2014) and following guidelines for interpretation from IUCN (2001). EOO comprises the minimum convex polygon covering all known or predicted sites for the species. It is frequently used as a preliminary assessment tool, such as when a new species is described or when populations of a species are found in new places or were locally extirpated (e.g., Lachenaud et al., 2013;Velzen & Wieringa, 2014). In the case of P. villosa, there has been no prior conservation assessment for the species, and its status is not presently included in the IUCN Red List of Threatened Species (IUCN, 2020). Therefore, we analyzed the EOO based on the 155 occurrence data points also used for ENM (Table S3).

| Polymorphism of AFLP markers
The 43 populations of P. villosa were all assigned correctly to the corresponding individual after amplification and scoring. The eight primer pairs yielded 1,728 clearly identifiable amplified bands, of which 1,654 (95.7%) were polymorphic (Table 1)

| Clusters of populations and genetic structuring
The UPGMA analysis on 43 populations of P. villosa showed that these populations could be assigned to two clusters that have high support (100% bootstrap support (BS), Figure 1). The gap statistical analysis based on all 210 individuals yielded five groups in which individuals generally clustered within their populations (Figure 2), and groups of populations represented subgroups on the UPGMA tree. Based on SAMOVA, we found that the optimal number of clusters of populations was three (K = 3) ( Table S3). As with the PCoA analysis (Figure 3), the three groups represented subgroupings of the two highly supported ones according to UPGMA. Specifically, separate groups of populations (P) 1-6 and P7-12 were resolved in SAMOVA. However, in the field, we observed striking similarities in habitat among P6-8, which were also geographically proximal.
Meanwhile, both SplitsTree and STRUCTURE revealed two clusters of populations (Figures 4 and 5). However, in the SplitsTree analysis, some individuals from populations 36, 37, 38, 39 did not cluster with others from their populations and were resolved in the opposing cluster (Figure 4). These same individuals also showed similar phenomenon in the analysis of PCoA. Overall, we chose to treat the populations as belonging to the two groups identified based on UPGMA, though we acknowledged that results from other analyses, such as the gap statistic and SAMOVA, suggested that additional structures with weaker signal might also exist among the sampled populations.
F I G U R E 1 Dendrogram of P. villosa generated by unweighted pair group method analysis (UPGMA) cluster analysis from the genetic similarity matrix obtained using amplified fragment length polymorphism genetic distance (see Figure S1 for population codes) Hereafter, we refer to the two groups as Groups 1 and 2, which comprise P1-12 and P13-43, respectively. Populations of Group 1 occurred mainly in the central and eastern regions of the Inner Mongolian Plateau, while populations of Group 2 were distributed throughout the range of the species in China. Notably, the populations of Group 1 tended to be found at lower elevations compared to those of Group 2 (Table S1 & Figure S1). showed that the populations of Group 1 exhibited minimally greater genetic diversity than those of Group 2 (except the value of Na). In addition, we noted a strong correlation between I and h (Table 1: Spearman ranking correlation, R = 0.976, p = 0; Table S1: Spearman ranking correlation, R = 0.973, p = 0) and between I and He (Table S1:  Table 2).

| Genetic diversity of populations
The Mantel test revealed that there was a significant positive correlation between geographic distance and F ST for the 43 populations (r = .282, p < .05) (Figure 6). Similarly, we detected a strong, significant, positive correlation between geographic distance and F ST for Group 1 (r = .622, p < .05) and a weak but significant positive correlation for Group 2 (r = .372, p < .05).

| Distributional change of P. villosa
ENMs for P. villosa yielded relatively high AUC, demonstrating reliable model performance (AUC = 0.969, Figure S2). For the eleven variables used for modeling, the most significant factor for the spatial distribution pattern of P. villosa was the altitude (Alt), followed by temperature annual range (Bio 7) and precipitation of warmest quarter (Bio 18), whose contribution rates were 40.0%, 17.2%, and 16.7%, respectively (Table S5). Based on model projections, we observed a contraction in a highly suitable habitat during the LGM compared with the LIG (Table 4 & Figure 7). Nevertheless, there was less suitable highly suitable habitat in the present compared to during either the LIG or LGM (Figure 7). Highly suitable habitat projected for the present day was largely congruent with the actual geographic distribution of P. villosa, within the Inner Mongolia Plateau.

F I G U R E 3
A two-dimensional plot of the principal coordinate analysis (PCoA) based on variation of amplified fragment length polymorphism markers for P. villosa (see Figure S1 for population codes; ellipse, 95% confidence interval) Simultaneously, we estimated the future changes in the potential spatial distribution under the RCP 2.6 and RCP 8.5 scenarios for the 2050s and 2070s. According to the future model predictions, our projections of the models based on future climates showed that, in general, the areas of suitable habitat for P. villosa would remain stable under the climatic scenario of RCP 2.6 for the 2050s and 2070s, whereas there was an increase in highly suitable areas based on RCP 8.5 (Table 4 & Figure 8).
When we compared the niches of Groups 1 and 2, we found that D and I were significantly lower than the values expected from the pseudo-replicated datasets. Thus, there is distinct niche differentiation between the two groups (p < .01) (Figure 9). The niches of the two groups differed mainly in that Group 2 occurred at higher elevations and under temperature annual range.

| Genetic diversity of P. villosa
Genetic diversity is closely linked to the evolutionary potential of a species to adapt to adverse environments . In the present study, we observed high genetic diversity at the species level in P. villosa (h = 0.21) and at the population level (h = 0.13).
The underlying drivers of genetic diversity within species are generally a combination of biological factors, such as dispersal abilities and life history, and environmental factors, such as climate and anthropogenic activities (Gerzabek et al., 2020;Prazeres et al., 2020).
The life history of P. villosa frequently involves clonal reproduction via its rhizomes under harsh environmental conditions, although the species also reproduce sexually by seed following wind pollination (Li & Ge, 2001;Wang et al., 1999). In comparison with L. chinensis, C. bulbosum, and D. glomerata, the relatively high genetic diversity of P. villosa might be explained by one or more factors. Among these, our study design comprised more populations, which might lead to greater accuracy in inferring genetic diversity. However, biological explanations are more likely and include possible higher clonal F I G U R E 4 Neighbor-Net split network of P. villosa based on amplified fragment length polymorphism datasets using Jaccard's distances. Lines of red and yellow represent Group 1 and Group 2, respectively fitness of P. villosa as it has extremely robust, hardy rhizomes, and high lifetime rates of seed production and regeneration via seedlings (although annual regeneration via seedling is at a low rate) (Eriksson & Bremer, 1993;Shimizu et al., 1998).

| Genetic differentiation and genetic structure
The genetic structure of a species is effectively the sum of genetic differentiation among and within populations (Hamrick & Godt, 1989).
Overall, genetic structure occurring among populations results from the evolutionary history of the species in question; natural selection; genomic factors (e.g., mutations, reorganization, and genetic drift); and biological characteristics, including gene flow, mating system, mode of reproduction, and seed dispersal mechanisms (Slatkin, 1987;Zhen, 2010). Genetic differentiation is primarily controlled by aspects of gene flow, such as its rate and directionality (Hamrick & Godt, 1989). In plants, gene flow occurs most often via the transmission of pollen and seeds during sexual reproduction. However, for clonal species, such as P. villosa, asexual propagules may be more common than seeds, but these usually have limited dispersal distance and, thus, restrict gene flow among populations (Xia et al., 2002).

F I G U R E 5
Results of the Bayesian clustering analysis in STRUCTURE of 210 individuals representing P. villosa. (a) ΔK values from the mean log-likelihood probabilities through STRUCTURE runs where inferred cluster (K) ranged from one to ten; (b) Estimated genetic clustering for K = 2, where unique colors correspond to assignment at different clusters; (c) Geographic origin from 43 populations of P. villosa and their colorcoded grouping according to the structure analysis for the model with K = 2 For P. villosa, we inferred that more than 56% of the genetic vari-   (Takhtajan, 1986) and contained a vastly richer flora than the broader, surrounding area (Jiang et al., 2007). The transitional nature of this area may be more likely to support fixation of genes introduced from outside populations; that is, among other populations outside of the Helan Mountains, dispersal of propagules may occur commonly, but selection favors local ecotypes. Overall, with limited effective gene flow, P. villosa has undergone considerable genetic divergence and exhibits a high level of genetic structure.
Although long-distance dispersals might be a one critical aspect of genetic structure in P. villosa, genetic distance was also significantly correlated with geographic distance based on a Mantel test. Therefore, genetic structure in P. villosa might primarily result from geographic isolation imposed by mountains (e.g., Yin Mountains; Helan Mountains) and large deserts in northwestern China (e.g., Tengger Desert; Mu Us Sandy Land) as well as range contraction and population fragmentation induced by climatic oscillations as also observed for Gymnocarpos przewalskii Maxim. and Helianthemum songaricum Schrenk (Ma et al., 2012;Meng et al., 2014;Su et al., 2011). In addition, founder effects and population bottlenecks might have also contributed to the genetic structures of the species (Birky et al., 1989;Liu et al., 2015).

| Demographic history of P. villosa
The genetic diversity within Group 2 was slightly lower than that of Group 1, despite that Group 1 comprises a smaller number of populations (12 vs. 31). Based on this, extant populations of this species might originate from the genetic stock of Group 1, as geographic areas with both high genetic diversity and frequency of dominant genes usually represent centers of origins for source populations (Vavilov, 1926). However, our study design and results cannot discern the exact center of origin for the species nor the main migrational patterns of P. villosa, and accomplishing this will require additional molecular data and informatics approaches.

F I G U R E 7
Potentially suitable climatic distribution of P. villosa under different climate change scenarios in the Inner Mongolian Plateau Climate oscillation during the Quaternary has often been hypothesized to be an important factor in influencing the current geographic distribution and demographic history of plant species (Hewitt, 2004;Su & Zhang, 2013). One widely utilized approach to comparing past and future distributions of plant species and determining the primary environmental factors driving them is ENM (e.g., Bai et al., 2017;Nabout et al., 2016;Wei et al., 2018). Specifically, our models showed that the range of P. villosa was the most Northern Hemisphere during the early-Middle Pleistocene (Meng et al., 2014;Xu et al., 2010;Yi et al., 2004). Nevertheless, it is surprising that the species range did not rebound as temperatures grew warmer following the LGM. This may be because of the onset of extreme aridity within the region during the Quaternary period, as this is widely known to have played a significant role in determining the geographic distribution and evolutionary history of many plant species (Meng & Zhang, 2011;Su & Zhang, 2013;Su et al., 2011). For example, in a previous study of Helianthemum songaricum (Cistaceae), which occured in Northern China and adjacent desert areas of central Asia (Yang & Gilbert, 2007), the worsening of the dry climate constrained the distributional range, and acceptable habitats for the species gradually became reduced and fragmented (Su et al., 2011).
Future studies of P. villosa may utilize genomic data and seek to understand the evolution of genes involved in adaptation to aridity.

| Germplasm conservation of P. villosa
Psammochloa villosa is a dominant species in its desert habitat, and sometimes it is the only herbaceous species occurring within its plant community. The species helps to maintain a fragile desert ecosystem by preventing wind erosion, the development of quicksand, and further desertification (Cai, 2016). After observing populations at 43 sampling locations during our field work, we noted that some populations of the species presently grow in severely degraded habitats. While we found that P. villosa is a species of least concern (LC) based on EOO (2,064,370 km²), habitat degradation may be an impending threat to the species and jeopardize its vital ecological role.
Therefore, we advocate for continued ecological monitoring of this dominant, keystone desert grass.
Psammochloa villosa may have great potential for sustainable utilization as a forage plant for livestock. The sand whips have relatively long inflorescences with large spikes that make it suitable for forage. Moreover, its adaptations to drought may make this species a valuable source of genetic resources for molecular breeding of other crop and forage species as, presently, it is one of few forage species that can withstand the intensifying long-term drought conditions in northwest China. Developing a sustainable use strategy for P. villosa will also help to ensure its continued availability as a keystone species within desert communities of the Inner Mongolian Plateau and adjacent areas.

ACK N OWLED G M ENTS
We greatly appreciate Dr. Jiabin Zou helping to improve an earlier draft of our manuscript. This work was financially supported by the National Natural Science Foundation of China (Grant Nos. 41761009 and 31800310), the Natural Science Foundation of Qinghai Province F I G U R E 9 Results of the nicheidentity test. (a) Schoener's D; (b) Warren's I. The arrow in each panel represents the observed niche similarity between occurrence points for the corresponding pair of clusters. The histograms represent the distribution of niche similarities obtained from pairs of pseudo-niches constructed by random resampling of occurrence points of the two clusters (Grant No. 2019-ZJ-7011), "The Dawn of West China" Talent Training Program of the Chinese Academy of Sciences (2019-1-4).

CO N FLI C T O F I NTE R E S T
None declared.