Population Genetic Structure and Biodiversity Conservation of a Relict and Medicinal Subshrub Capparis spinosa in Arid Central Asia

: As a Tertiary Tethyan relict, Capparis spinosa is a typical wind-preventing and sand-ﬁxing deciduous subshrub in arid central Asia. Due to its medicinal and energy value, this species is at risk of potential threat from human overexploitation, habitat destruction and resource depletion. In this study, our purpose was to evaluate the conservation strategies of C. spinosa according to its genetic structure characteristics and genetic diversity pattern among 37 natural distributional populations. Based on genomic SNP data generated from dd-RAD sequencing, genetic diversity analysis, principal component analysis, maximum likelihood phylogenetic trees and ADMIXTURE clustering, the signiﬁcant population structure and differentiation were explored. The results showed the following: (1) Six distinct lineages were identiﬁed corresponding to geographic locations, and various levels of genetic diversity existed among the lineages for the natural habitat heterogeneity or human interferences; (2) The lineage divergences were inﬂuenced by isolation by distances, vicariance and restricted gene ﬂow under complex topographic and climatic conditions. Finally, for the preservation of the genetic integrity of C. spinosa , we suggest that conservation units should be established corresponding to different geographic groups, and that attention should be paid to isolated and peripheral populations that are experiencing biodiversity loss. Simultaneously, monitoring and reducing anthropogenic disturbances in addition to rationally and sustainably utilizing wild resources would be beneﬁcial to guarantee population resilience and evolutionary potential of this xerophyte in response to future environmental changes.


Introduction
Arid central Asia is deemed to be the largest arid region in the temperate zones of the northern hemisphere and even the world [1]. The floristic compositions here are mostly descended from Tethyan xerophytic vegetation flora and have been phenoecological plant taxa shaped by the long-term continental-arid desert climate [2,3]. Because of the specific geomorphologic landscapes comprising the vertical and horizontal convergence of several major Asian tectonic mountains (Himalayas, Karakoram, Kunlun and Tianshan Mountains) and plateaus (Pamirs and Tibet Plateau), as well as the mosaic distribution of many large basins (Tarim, Turpan-Hami and Junggar basins) and deserts (Taklimakan, Kumtag and Gurbantunggut deserts), the monsoon current and moisture from the Indian and Pacific Oceans have been blocked and intercepted, and the amount of precipitation is low and unevenly distributed, which has constituted the complexity of the extreme geographic and climatic conditions, consequently forming a very fragile ecological environment [1]. The types of vegetation are relatively sparse and simple in this region. To adapt to environment variability, desert species from oligotypic genera and a monotypic genus with drought resistance and barren tolerance have mainly occurred and been discontinuously A number of previous studies have hitherto been carried out on the genetic variation and diversity of C. spinosa [15][16][17][18][19][20][21]. Earlier on, scholars mainly devoted their efforts to morphological taxonomy and subspecies or variety identification [15][16][17]. In recent years, with the development of molecular marker technology, researchers have achieved certain advancements in intraspecific diversity and infraspecific differentiation by virtue of methodologies such as RAPD, AFLP, ISSR, IRAP and EST-SSR [18][19][20][21]. Nonetheless, most of these studies have focused on Mediterranean coastal regions with a subtropical dry summer climate. However, the relict wild populations remaining in arid central Asia, dominated by a temperate continental desert climate, have been shown a lack of concern. Valuable information on the genetic structure and diversity conservation of C. spinosa in this region is relatively scarce and is of significant biogeographical and conservation interest to phytologists in the fields of plant genetic diversity and biodiversity conservation.
The conservation of genetic diversity is of great importance to the maintenance of biological species and ecosystem diversification, which provides a crucial foundation for the survival and development of species and enhances their long-term adaptive evolution to environmental changes, reducing the risk of extinction [22][23][24]. The evaluation of levels of genetic diversity based on integrated population structures is a critical prerequisite for species protection [25]. For source-sink metapopulations of terrestrial plants, the genetic structures and divergences are not only dominated by their own species and population characteristics, such as growth rates, effective population sizes and dispersal capacities, but are also closely related to their demographic responses to external spatial and temporal variability of environment [26,27]. Drivers of isolation by distances (IBD), vicariance and landscape heterogeneity influence the species' evolutionary processes (genetic drift effects, gene exchanges, natural selection and local adaptation), thereby determining the genetic diversity patterns of metapopulations [28,29]. Human disturbance is also a contributory external factor influencing biodiversity. Since the beginning of the Anthropocene in the latter part of the 18th century, hazards of wild habitat destruction and loss caused by expansions in the range of anthropogenic activities have been becoming increasingly serious [30,31]. Under the background of global environmental and climatic changes, mankind's excessive pursuits of economic benefits are causing the growing reduction in biological resources, the extinction of species and populations, as well as the modification of natural landscapes [32]. Notably for small-scale isolated or peripheral populations, due to lack of connectivity with source populations, they are particularly prone to genetic depletion through genetic drift and inbreeding [33], whereas anthropogenic interference undoubtedly further weakens their fitness, decreasing survivorship and reproduction, and increasing extinction risks [34]. Therefore, research on and the protection of these sink populations are considered to be essential for maintaining the genetic integrity of large-scale metapopulations. Hitherto, the conservation of genetic resources has mainly focused on native or endemic rare and endangered species in global biodiversity hotspots, such as the Andes Mountains, Mediterranean Basin, eastern Himalayas and Hengduan Mountains in the tropics and subtropics [35][36][37][38]. However, some widespread natural resources, typically with medicinal or economic value, experiencing diversity reduction in temperate arid regions have not been sufficiently recognized and effectively protected [39]. If this continues, these species would not escape from the hidden dangers of resource exhaustion or even extinction.
Compared with traditional molecular genetic markers, genome-wide single-nucleotide polymorphisms (SNPs) have the characteristics of higher genotyping efficiency and data quality, analytical simplicity, broader genome coverage and sheer abundance, as well as advantages in further analyzing complex population structures and distinguishing the evolutionary relationship of genetic polymorphisms [40,41]. In recent years, they have been increasingly widely applied in research fields of phylogenetics, molecular ecology, conservation genetics and biogeography [41][42][43]. As a high-efficiency next-generation highthroughput sequencing technology (NGS) that has been developed and commonly used, restriction-site-associated DNA sequencing (RAD-Seq) can rapidly identify and score highresolution SNP data from thousands of orthologous sites, while reducing the complexity of genomic sequences, even in the absence of reference genomes [44,45], becoming one of the main effective methods to address critical issues in the fields of evolution and genetics at the genomic level [46], especially for studies of species groups with ancient origins and complex evolutionary and genetic structures [47,48].
In the present study, for the first time, we used the SNP marker from dd-RAD sequencing, ADMIXTURE clustering, principal component analysis (PCoA) and maximum likelihood (ML) phylogenetic analysis methods to fully explore the integrated genetic diversity mechanisms and conservation significance of C. spinosa across large-scale environmental gradients in arid central Asia, based on spatial source-sink theory of metapopulation dynamics. Thus, we aimed to: (1) elucidate the population structure and intraspecific divergence of the lineages of this species; (2) to adequately reveal the spatial pattern of genetic diversity among identified geographical groups and, accordingly, formulate feasible conservation strategies for maintaining the genetic integrity and population resilience of this xerophytic subshrub; and (3) to provide a referential understanding for future research on the diversity conservation of relict plants in arid regions.

Sampling, DNA Extraction and dd-RAD Library Construction
We located and collected a total of 37 natural populations (5-12 individuals for each population) of C. spinosa in arid central Asia during 2019-2020 ( Figure 2, Table S1), according to the distribution described in Flora of China, Flora Xinjiangensis and Flora of Tibet as references [12][13][14]. The sampling spacing was at least 100 m between the neighboring individuals within each population. We also collected Capparis erythrocarpos Isert (northeastern Kenya, Africa) and Capparis bodinieri Lévl. (Xishuangbanna, China) as outgroups for the subsequent phylogenetic analysis. Total genomic DNA was extracted from silica-dried leaf tissues using a DNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA), following the manufacturer's protocol. The cleaned DNA had previously conducted simulated double digestion with 15 combinations among seven common or rare restriction enzymes to predict enzyme cutting sites (Table S2). Considering the yielded tag numbers and average sequencing depths, we finally selected the optimal combination of BfaI (5 -C↓TAG-3 ) and DpnII (5 -↓GATC-3 ) (New England Biolabs, Ipswich, MA, USA) to construct the dd-RAD library, following a modified approach of Peterson et al. [45]. The library with an inserted fragment between 200 and 600 bp was recovered for sequencing. All the qualified samples' sequencing was performed on the Illumina Hiseq X Ten platform (Illumina, San Diego, CA, USA) with paired-end reads of 2 × 150 bp in length. Library preparation and sequencing were performed by Shanghai Personalbio Biotechnology Co., Ltd.

Data Processing and SNP Calling
Raw data produced by Illumina were processed to yield RAD tags and to obtain SNPs using the Stacks v1.4.8 program [49]. Initially, paired-end reads were cleaned and filtered for quality control and then trimmed to 140 bp × 2 in length using the process_radtags module. Individuals with sufficient RAD tags and coverage depth were used for SNP calling. Since there were no available reference genome data for C. spinosa, we next chose the ustack module to assemble reads per individual into RAD tags following the protocol of de novo workflow [49]. To remove high-repetitive stacks via a deleveraging algorithm, we set coverage depth to at least m = 3, the maximum distance allowed between stacks as M = 2 and the maximum number of mismatches allowed between loci as n = 1. Then, we combined the consensus loci among all 263 samples to build a catalog file in a cstacks module and aligned each loci sequence with the catalog to generate alleles in an sstacks module. Finally, filtered SNP matrices were identified and exported using a populations module, where the optimal operation parameters were set as follows: complete degree > 0.5 (individuals where the same RAD locus was detected accounted for at least 50% of each population), minor allele frequency (MAF) > 0.05 and minimum stack depth m = 4, p = 100 (one locus appeared in at least 100 populations). Based on the above calculations, the resulting high-quality SNP dataset was used for further analysis.

Genetic Structure Analysis
Population genetic structure, phylogenetic analysis and principal component analysis (PCoA) were comprehensively considered to evaluate the spatial distribution pattern of genetic variation. We performed a maximum likelihood estimation to analyze the genetic structure of C. spinosa using the ADMIXTURE program, version 1.3.0 [50], in order to identify the group clustering of each individual based on unlinked SNP datasets via a block relaxation algorithm. When the numbers of grouping (K value) were preset as 1 ≤ K ≤ 10, ten replicated runs were implemented in ADMIXTURE to compute the mean cross-validation (CV) error for each K. The optimal K value was determined in conditions of the lowest CV error. For the analysis of phylogenetic relationships, we then reconstructed a maximum likelihood (ML) tree based on the unlinked genomic SNPs among these 263 individuals using the PhyML program, version 3.0 [51]. Capparis erythrocarpos Isert and Capparis bodinieri Lévl were chosen as the outgroups to root the tree. After the tree topology with maximum likelihood was established, the reliability of the branches was verified through a bootstrap (BP) with 1000 replications. Lastly, principal component analysis (PCoA) was performed to exhibit the first two axes of population variation in the species using GCTA software, version 1.93.2 [52]. After removing SNPs with MAF values of less than 0.05, genomic SNP data of different individuals were formed into a genetic matrix where the feature vectors that made a major contribution to variance were extracted, and then a scatter plot was drawn using the first two of the eigenvectors. According to the distribution characteristics of the principal components, the clustering relationship of individuals could be inferred.

Genetic Diversity and IBD Analysis
Based on the results of genetic grouping, population genetic statistics, including observed heterozygosity (H o ), expected heterozygosity (H e ), nucleotide diversity (P i ) and inbreeding coefficient (F IS ), were computed using the populations command in the Stacks v1.4.8 package [49] at the group and population levels of the species. Degrees of deviations from the Hardy-Weinberg equilibrium were measured for all 37 populations. Hierarchical AMOVA (analysis of molecular variance) in the ARLEQUIN version 3.5.2.2 program was used to classify the total genetic variation partitioning within and across groups by calculating the evolutionary distance between alleles or genotypes [53]. Pairwise genetic differentiation indices (F ST ) among populations were also calculated in the populations module in the Stacks. Additionally, the significance of isolation by distance (IBD) was further tested to compare Pearson's correlation between the genetic distance (F ST ) matrix and geographic distances matrix among these 37 populations using the cor.test function in the R v3.3.3 package [54].

SNPs from dd-RAD Analysis of C. spinosa
A total of 263 C. spinosa individuals from 37 populations were sequenced using the dd-RAD library construction methodology. Both the base content and quality distribution across all the sequences were uniform and optimal, considering that the average proportion of bases whose recognition accuracy exceeded 99.9% (Q30) was greater than 90% and that the peak value was around Q36. The numbers of raw reads for each sample ranged from 6,455,458 to 26,065,222. After filtering out low-quality reads, the numbers of high-quality reads ranged from 5,995,376 to 24,640,428, and the mapping rate ranged from 90.24% to 97.98%. According to the clustering algorithm through the Stacks procedure, the obtained tag number per sample ranged from 247,350 to 823,581, and the average depth of sequencing was from 10× to 31×. Via mutation spectrum analysis, C:G > T:A and T:A > C:G account for the overwhelming majority of mutation types of SNPs ( Figure S1). The final dataset, which contained 559,586 high-integrity SNPs, could be used for the subsequent population genetics analyses.

Population Genetic Structure of C. spinosa
The ADMIXTURE analysis showed that, when K = 6, the minimal value of the crossvalidation (CV) error was detected ( Figure S2), indicating that the best-fit number of genetic groupings of C. spinoca was six: (1) Table S1). This pattern of population genetic structure was also supported by the topological relationship of the ML phylogenetic tree, with the main nodal bootstrap values exceeding 89.3% at the intraspecific level ( Figure 4). Populations from group WH (red) always diverged first from all other populations, whether in the ADMIXTURE clustering or in the phylogenetic tree (Figures 3 and 4). Then, group EP (blue) branched off, followed by group WT (yellow) and group ST (green). Group NT (azure) was nested within group ST, and populations from group ET (purple) formed one clade ( Figure 4).  The first two coordinate axes of PCoA identified three major clusters, and corresponding principal components explained 16.55% (PC1) and 6.99% (PC2) of the total genetic variation, respectively (Figure 5a). Considering the preliminary results, group WH formed the first cluster, which was completely separated from the rest of the groups; group EP and group WT were generally merged into the second cluster; and group ST, group NT and group ET were broadly integrated into the third cluster. Then, we conducted further PCoA analyses on the latter two clusters, separately. The substructure within the second cluster showed that group EP was gradually split from group WT, with 16.15% (PC1) and 6.99% (PC2) of the total variation (Figure 5b). The substructure within the third cluster showed that populations in groups ST, NT and ET were mostly divided from each other (PC1, 16.15% vs. PC2, 6.99%), whereas admixed populations crossed among them (Figure 5c), which was also evident in the ADMIXTURE diagram ( Figure 3).

Genetic Diversity Pattern and IBD
The results of the genetic diversity analysis showed that the overall genetic diversity among geographical groups was ranked as follows: EP > ST > NT > WT > ET > WH, which exhibited distinct differences in various degrees (Table 1). At the population level, the maximum observed heterozygosity (H o ), expected heterozygosity (H e ) and nucleotide diversity indices (P i ) were in Manas (C26) and Jiashi (C10), while the minimum H o , H e and H e values were in Zabrang (C01) ( Table 1). Inbreeding coefficients (F IS ) that measured degrees of deviation from the Hardy-Weinberg equilibrium were positive for all 37 populations. The results of the AMOVA showed that the percentage of variation among the six geographical groups was relatively low; the variation was mainly from individuals within the populations, and the variation rate across populations within groups was the lowest ( Table 2). Significant genetic differentiation existed among geographic popula-tions ( Figure S3, Table S3), with F ST values on a scale from 0.076 to 0.707 (Table S3). The maximum genetic differentiation was detected between C02 (Sarang population) and K01 (Bishkek population) (F ST = 0.707), while the minimum genetic differentiation existed between C27 (Shihezi population) and C29 (Wusu population) (F ST = 0.076) (Table S3). Nei's genetic distance was roughly synchronous with pairwise F ST between populations ( Figure S3). Furthermore, the result of isolation by distance (IBD) analysis displayed a considerable correlated pattern between genetic divergence and geographical distance among these 37 populations, given that the two-tailed test of the Pearson correlation coefficient was of significance (r = 0.499, p < 0.001) ( Figure 6).

Relationships between Genetic Structure, Lineage Differentiation and Geographic Distribution
In this study, we obtained a clear genetic structure of C. spinosa through RAD-seq and identified six clustering results corresponding to geographic locations (Figures 2 and 3, Table S1). For this Tertiary relic plant, the effects of isolation by distance and vicariance has caused long-term shaping of local ecological environments, which had a profound influence on population structure and lineage differentiation among different geographic units. It is speculated that the formation of this disjunction is the result of intense and frequent orogeny during the geological period, coupled with extensive Quaternary glaciation, as well as the long-term process of aridification in central Asia [5][6][7]. The allopatric divergence and local adaptive evolution among distant geographical groups, especially for the isolated populations, are probably attributed to the incapability of long-distance genetic exchange and being hindered by mountains and immense deserts. Similar disjunctive or fragmented population structure patterns have also been reported in other Tethyan relicts that occur in this area, such as Gymnocarpos przewalskii, Amygdalus mongolica and Ammopiptanthus, owing to the influence of long-term aridification [9][10][11]. The constantly extended deserts and climatic changes have led to the fragmentation and separation of suitable habitats for those xerophytic taxa, thereby restricting the migration among distant geographic populations [55].
The relatively isolated populations distributed in the northern piedmont of the western Himalayas (Group WH, C01-C03) have the maximum genetic distance from other populations. Plant taxon located in this desert plateau region is related to the sub-frigid arid climate, precipitous terrain and topography, as well as abnormally barren soil conditions (Figure 1a). Due to the complex and discontinuous natural environment, the sink populations here, to a certain extent, are differentiated by hereditary and habitat variability. Compared to populations distributed at high altitudes, populations in medium-altitudinal valleys and low-altitudinal desert basins have different habitat requirements (Figure 1). Among the eastern Pamir region and the southern piedmonts of the Tianshan Mountains, we found that population structures were associated with spatial geographic distribution, as well as heterogeneous environments [29]. For instance, the Jiashi population (C10) is geographically nearby the Kalpin (C12) and Tumxuk (C11) populations, but due to the hereditary differences and landscape heterogeneity (gravel Gobi vs. dried clay desert) they belong, respectively, to group EP and group ST, although some of the mixed genotypes were found in these divergent populations (Figure 3). On the contrary, according to research results from Mims et al. (2016), among populations within the same geographic unit, the habitat preferences are of continuity and similarity, where high degree of gene flow could exist among them in a relatively flatter and more expansive range [56]. For example, source populations that inhabit piedmont proluvial fans on the southern side of the Tianshan Mountains (Group ST, C11-C17) are more closely distributed, and the connectivity of high-quality breeding habitats may play a role in facilitating closer reduced genetic distance between these populations.

Genetic Diversity Pattern of Species Metapopulations
The genetic diversity of C. spinosa (Pi = 0.2644) ( Table 1) through RAD-seq is higher than that of other relic species (such as Euptelea pleiosperma and Gymnocarpos przewalskii) [11,48] because widespread species have a higher genetic diversity than species with endemic distribution. The total diversity index in arid central Asia is also different from that obtained by previous research using AFLP, RAPD and EST-SSR for C. spinosa in the Mediterranean coast (southern Europe, north Africa and west Asia) [18,20,21]. The primary reasons may be the different systematic molecular markers and the varying regional environments and weather conditions. In the arid region of central Asia, the distribution of C. spinosa covers most of the mountain piedmont, warm dry valley and Gobi desert terrains. Based on the results of SNPs, the range of this species is divided into six geographic groups with various levels of genetic diversity, indicating the levels of geographic isolation and environment gradients. According to the observed heterozygosity, the genetic diversity of source populations in groups ST and NT remained at relatively high and stable levels ( Table 1), indicating that populations within these two groups are both likely to have experienced long-term migration of gene flow. In these geographic units, the high-quality piedmont alluvial or diluvial fan habitats are usually continuous, which is likely conducive to maintaining a high level of genetic diversity throughout these regions. In contrast, sink populations in group WH (western Himalayas) inhabit low-quality desert plateaus with high altitudes and low temperatures, especially the isolated Zabrang population (C01) with lowest genetic diversity (Table 1) adjacent to the middle reaches of the Sutlej River Basin, which grows on steep rocky slopes and seems a bit far away from the water source ( Figure 1a); even worse, natural disasters, such as landslides and mudslides, occur frequently during the annual flowering and fruiting phases. In addition, because of insufficient visits by pollinators, inbreeding or self-breeding dominate the reproduction of C. spinosa under this situation, resulting in significant population depression (current effective size of <50 individuals). Many previous studies have demonstrated that, in nature, isolated populations with smaller effective sizes often have more genetic depletion and extinction risks due to genetic drift and inbreeding, which are, to a large extent, disadvantageous to the persistence of population diversity [57,58]. In the process of evolutionary history, the Zabrang population survives in a separate habitat but has not adapted to the extreme environment, which may be one of the reasons why its population scale is declining and even endangered. Nevertheless, it is noteworthy that the Taxkorgan (C04, 2855 m alt) and Akto (C05, 2845 m alt) populations, which are also located at high altitudes, could maintain certain high levels of genetic diversity because they are not completely isolated in the eastern Pamirs region (group EP). These populations have continuous extension and buffer along the descending elevation gradients, and they also have a certain degree of gene exchange with the Yegisar population (C06) in the piedmont belt.
Researchers have proved that, compared with diffusion within short distances, the long-distance dispersal of gene flow plays a more pivotal role in the persistence of the entire connectivity and genetic diversity of plant metapopulations [59,60]. Simultaneously, the contributions of wind, pollinators and river water help the pollen flow to spread farther. However, within some geographic groups (group ET, group WT), moderate degrees of genetic differentiation (F ST > 0.1) (Table S3) have appeared in neighboring populations, indicating that anthropogenic activities may have encroached into these natural territories, leading to the discontinuous gene flow. Driven by economic interests and commercial purposes, human over-exploitation has weakened the reproductive capacity of C. spinosa populations and has arbitrarily cut off the long-distance dispersal chains of gene flow, making short-distance spreading predominant, increasing the rates of inbreeding depression and resulting in population decline, thereby making it difficult for them to maintain their genetic diversity and evolutionary potential. Some populations near human activity areas, such as Toksun (C18), Hami (C21), Dunhuang (C22), Yining (C31) and Gongliu (C32), are frequently threatened by natural resource plundering and other disturbances, such as riverway-or highway-widening, agricultural reclamation and industrial expansion, etc. The relatively low genetic heterozygosity and diversity levels ( Table 1) found in these populations may be primarily caused by the results of the above depredations.

Implications for Conservation
Evaluating the spatial genetic structure and survival situation of metapopulations across regional environment variability could efficiently and effectively support and improve the management measures of biodiversity conservation, thereby conducing to maintain the evolutionary potential of widespread species in response to changing environments [33,56]. In this study, most of the central source populations with large effective sizes are located in the expansive piedmonts of north and south Tianshan Mountains, where multiple hotspots of biodiversity have been shaped in their long evolution history in the arid region. In contrast, isolated populations on the desert plateau of western Himalayas, as well as peripheral populations in the desert basin adjacent to eastern Tianshan, are experiencing a decrease in or even disappearance of biodiversity. As a result of the adverse effects of hereditary factors or human disturbance, these sink populations are extremely vulnerable to habitat fragmentation and resource exhaustion [34], so it is particularly important to monitor population trajectories, minimize human disturbances and restore natural landscapes [27]. Given the limited manpower and material resources, in addition to relevant policies for protecting xeromorphic vegetations in arid regions, we recommend that conservation management units should be established that are consistent with the six effective evolutionary groups and that the emphasis of preservation work should be laid on the genetic rescue of isolated populations and habitat restoration of peripheral populations in order to persist the genetic integrity of this plant species.
Conservation of sink populations has been described as an important part of biodiversity persistence of integrated metapopulations because they often contain certain rare alleles or genotypes in the process of historic immigration [61]. The Zabrang population (C01) is distributed in more sparse and localized low-quality habitats, resulting in its greater genetic differentiation with other populations within group WH (F ST = 0.1357-0.1878) (Table S3). However, due to a lack of corresponding attention and conservation, the genetic diversity and uniqueness of such isolated population are being threatened and almost disappearing. Empirical studies have illustrated that inputting a certain level of stable gene flow from large source populations to small sinks may counter the negative effects of inbreeding depression and genetic drift [26,33]. Currently, we suggest that genetic rescue such as immigration, exchange of pollen flow and hybridization, might reduce the extinction risk of isolated populations of C. spinosa. If conditions permit, ex situ preservation and reintroduction instead of natural dispersal may be more effective protection strategies to increase local effective population sizes.
The disturbance of human activities, such as the expansion of farmland, makes this subshrub lose its competitive advantage, which could consequently lead to its habitat shrinking. Additionally, the industrial waste generated by newly built factories in the suburbs has caused serious pollution to the soil environment and increased the salinealkali stresses on its habitat, whereupon the mortality rates of plant individuals have also risen. If these trends persist, the relic species would scarcely escape from the hidden danger of resource exhaustion or even the verge of extinction. Populations in the Tuha Basin-Hexi Corridor (Group EP) are located on the edge of the Kumtag Desert, where water shortages and soil salinization are becoming increasingly serious. These findings are the consequences of continuous increasing human pressure on the barren habitats of xerophytes. According to our field investigation conducted from 2011 to 2020, the population sizes within group EP in the desert basin were originally large (>300 individuals per population in 2011). Nevertheless, over the past decade, due to humans' immoderate harvesting and unceasing expansion of agriculture and industry, the living space of C. spinosa in this region has decreased sharply. In particular, the genetic diversity of the peripheral Hami and Dunhuang populations (C21 and C22) has been gradually decreasing, and the plant morphology and population sizes (<100 individuals per population in 2020) here have also been shrinking under these dual pressures, which weaken the plant's ability to prevent wind and sand and preserve water and soil (Figure 1c). This region has been trapped in a vicious circle of natural desert vegetation destruction, water loss and soil erosion, as well as aggravated saline-alkali desertification. In order to realize the sustainable development of the ecological economy in arid regions, we should rationally conduct agricultural production on the basis of protecting natural desert landscapes, reduce the pollution of industrial waste to primitive habitats and limit the predatory exploitation of germplasm resources, therefore protecting the biodiversity of wild drought-resistant species to ultimately slow down the processes of desertification and salinization.

Conclusions
In this study, we analyzed the population genetic structure and biodiversity conservation management of C. spinosa, a Tethyan relic medicinal subshrub in arid central Asia. Six geographical units with differences in genetic diversity and lineage differentiation were distinguished. This excellent wind-preventing, sand-fixing and soil-and water-conserving xerophyte is at potential risk of natural resource depletion. Populations at high altitudes of the western Himalayas have been experiencing a gradual decline in genetic diversity and effective scales as a result of inbreeding depression by IBD and habitat isolation. On the other hand, peripheral populations at low altitudes in eastern desert basins with original rich genetic diversity are suffering from the plundering of germplasm resources and the destruction of survival environment. Therefore, to restore the genetic integrity of the metapopulations, we suggest that central source populations should be preserved in situ corresponding to different geographical units and that conservation priority should be focused on the genetic rescue of isolated populations and the habitat restoration of peripheral populations. Meanwhile, minimizing human activities in addition to the rational and sustainable utilization of natural resources will be of significance to the population resilience and evolutionary potential of C. spinosa in response to the long-term aridification and future changing environment, and it will eventually maintain the ecosystem balance, as well as slow the desertification process in arid regions.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/d14020146/s1, Table S1: Details of geographical locations, altitude and effective sizes of natural populations of C. spinosa. Table S2: Digestion sites of several common or rare restriction enzymes, as well as their combinations for screening schemes in this study. Table S3: Pairwise genetic differentiation indices (F ST ) among the 37 populations. Figure S1: Mutation spectrum of genome-wide single-nucleotide polymorphism (SNP) dataset of C. spinosa. Mutation types of SNPs are distinguished by different colors. T:A > C:G and C:G > T:A account for the overwhelming majority. Figure S2: Distribution of cross-validation (CV) error in the ADMIXTURE analysis. K = 6 is the optimal K value in condition of the lowest CV error. Figure

Institutional Review Board Statement: Not applicable.
Data Availability Statement: All RAD-seq data of C. spinosa were submitted to the NCBI (SRA accession: SAMN24698411-SAMN24698673).