A total of 633 individuals were sampled from 52 natural populations of N. sphaerocarpa and N. tangutorum in the study, including 119 individuals from 14 populations of N. sphaerocarpa, and 514 individuals from 38 populations of N. tangutorum. These populations almost cover the whole geographical ranges of the two species except a few populations of N. sphaerocarpa occurring in Mongolia ( Table 1 ; Figure 1A ). However, this will not bring a significant influence on the inference of the population structure of the two species, since these populations in Mongolia just occur on the fringe of the whole distribution of N. sphaerocarpa and represent only a small part of the gene pool in this study. In the sampled populations, five population sites are overlapped for the two species, for example populations 11 of N. sphaerocarpa and 12 of N. tangutorum share the same site, and 15 and 16 share another site ( Table 1 ; Figure 1A ). In each population, fresh leaves were collected from 5-30 individuals at least 100 meters apart and rapidly dried with silica gel in the field. The site information of each population, including latitude, longitude and altitude, were recorded using an Etrex GIS monitor ( Table 1 ). In addition, we collected a voucher specimen for each population, and some seeds and flowers for further studies on morphology and germination.

Total genomic DNA was extracted using DNeasy Plant Mini Kits (Qiangen, Valencia, CA, USA). Ten samples from 10 geographically distant populations were selected for primer scanning with primer pairs of six cpDNA fragments, i.e. trnL-trnF, rpl32-trnL (Taberlet et al., 1991), rpS-trnK (Shaw et al., 2007), trnH-psbA, atpH-atpI, and ndhC-trnV (Dong et al., 2012; Qian et al., 2016). Two cpDNA fragments, trnH-psbA and atpH-atpI, were found to contain useful and appropriate polymorphic loci, and further used for surveying genetic variations of the two species in the present study. Polymerase chain reaction (PCR) amplification was performed in a Bio-Rad PCR system S1000 Thermal Cycler (USA). The reaction volume was 25 µL, containing 1.0 µL genomic DNA extract (about 20-40 ng DNA), 2.5 µL 10×PCR buffer, 0.2µL Taq DNA polymerase (5U/µL, TakaRa Biotech Co., Dalian, China), 0.5 mmol/L dNTPs, 1.5 mmol/L MgCl₂, and 2µmol/L each primer. The thermal profile for trnH-psbA included the following phases: 94°C for 5 min, followed by 30 cycles at 94°C for 45 s, 56°C for 45 s, and 72°C for 105 s, then a final extension phase of 7 min at 72°C. The thermal profile for atpH-atpI was as follows: 94°C for 3 min, followed by 30 cycles at 94°C for 30 s, 58°C for 30 s, and 72°C for 60 s, then ended with an extension step of 5 min at 72°C. The PCR products were directly used for sequencing on an ABI 3130xl Genetic Analyzer platform (Applied Biosystems, USA), with forward and/or reverse primers. The obtained DNA sequences were aligned with Clustal_X (Thompson et al., 1997), and carefully revised manually. All the newly obtained sequences of N. sphaerocarpa and N. tangutorum were deposited in EMBL GenBank (https://www.ncbi.nlm.nih.gov/) under the accession numbers of OR801339-OR801365.

Nucleotide diversity parameters, including number of haplotypes, nucleotide diversity (Pi), number of polymorphic sites, were calculated using DnaSP 5.10 software (Librado and Rozas, 2009). We used NETWORK version 4.2.0.1 (Bandelt et al., 1999; available at http://www.fluxus-engineering.com) to construct Median-joining network of cpDNA haplotypes based on the data generated through DnaSP. We also reconstructed the phylogenetic tree among these haplotypes, using one individual of Zygophyllum xanthoxylum as outgroup. The neighbor-joining (NJ) tree, maximum-parsimony (MP) tree and maximum-likelihood (ML) tree, with bootstrap values of 1000 replicates (Felsenstein, 1985), were obtained respectively through MEGA 4 software (Tamura et al., 2007).

Unbiased genetic diversity (H _E) was estimated for each population based on cpDNA haplotype composition (Nei, 1987). According to the methods described by Pons and Petit (1996), average gene diversity within populations (Hs), total gene diversity (H _T) and population differentiation between populations (G _ST and N _ST) were calculated for the two species using the PERMUT program (available at http://www.pierroton.inra.fr/genetics/labo/Software/Permut/). G _ST depends only on the frequencies of the haplotypes while N _ST takes into account the similarity between haplotypes as well as haplotype frequencies (Pons and Petit, 1996). G _ST and N _ST were compared using the u-statistic to determine the presence of phylogeographic structure. A higher N _ST than G _ST denotes a significant phylogeographic structure, indicating that some closely related haplotypes usually occur in the same or nearby geographical regions (Pons and Petit, 1996). Analyses of AMOVA were carried out with Arlequin version 3.0 to estimate genetic differentiation within populations, between populations within groups and between groups, with significant tests of 1000 permutations (Excoffier et al., 2005). The correlation between geographic distance and genetic differentiation among populations was tested using the Mantel test through TFPGA version 1.3 (Mantel, 1967).

To further reveal spatial genetic structure of the two species, we performed a spatial analysis of molecular variance using SAMOVA version 1.0 (Dupanloup et al., 2002, http://web.unife.It/progetti/genetica/Isabelle/samova.html) based on the cpDNA haplotypes and geographical localities of all populations sampled in the study. In the analysis, the number of population groups (K) was defined using a simulated annealing approach, and K values were set between 2 and 12 with each simulation starting from 100 random initial conditions. A genetic differentiation index (F _CT) among groups was calculated, and the optimal configuration of groups was determined using an iterative simulated annealing process (Wright, 1978). In addition, biogeographical boundaries or potential gene flow barriers among populations were examined with Monmonier’s maximum-difference algorithm in BARRIER v2.2 (Manni et al., 2004). The robustness of these barriers was evaluated through running BARRIER program on the basis of 100 replicates of population average pairwise difference matrices. The difference matrices in the test were obtained through bootstrapping of haplotype sequences in SEQBOOT (Felsenstein, 2005).

Firstly, we carried out mismatch distribution analyses (Schneider and Excoffier, 1999) in DnaSP 5.10 version (Librado and Rozas, 2009) to test sudden demographic expansions of the two species. In general, populations with a sudden demographic expansion in the past should display a smooth and unimodal distribution (Slatkin and Hudson, 1991; Rogers and Harpending, 1992). Secondly, we performed Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) statistical tests using Arlequin version 3.0 (Excoffier et al., 2005) to infer historical demographic expansions of N. sphaerocarpa and N. tangutorum respectively. The validity of the expansion model was determined with the sum of squared deviation (SSD) between the observed and expected mismatches. The significance of expansion model (Harpending, 1994) was tested using Harpending’s raggedness index and associated P values. We used the formula τ = 2ut (Rogers and Harpending, 1992) to estimate the expansion time (t). In the formula, τ is the mode of mismatch distribution defined with units of evolutionary time, and u is the mutation rate per generation, which was estimated through the relationship u = 2 μkg. Here μ refers to the substitution rate per nucleotide site per year (s/s/y), k is the sequence length, and g represents the generation time in years. The average sequence mutation rates of chloroplast DNA in angiosperms (1.0-3.0×10^-9s/s/y) were used to estimate the expansion time of the species (Wolfe et al., 1987), since there was no fossil to calibrate cpDNA mutation rates of Nitraria species. According to our field observations and some cultivation experiments, the generation time of the two species was about 3 years. Thirdly, we performed ecological niche modeling (ENM) on the basis of maximum entropy approach to infer the potential distribution ranges of the two species in the LIG (ca. ~140Kya), the LGM (ca. ~21 Kya), the present, and the future (2050), using the MAXENT version 3.3.3k (Phillips et al., 2004; 2006; Peterson and Soberón, 2012). In the simulating analysis, a “maximum training sensitivity plus specificity” threshold was used to determine suitable/unsuitable habitat (Liu et al., 2010). A jackknife test was carried out to estimate the contributions of different bioclimatic variables to the prediction of the distributional models. The goodness of models was evaluated with the area under the receiver operating characteristic (ROC) curve (AUC scores) (Fawcett, 2006). In order to reduce effects of spatial autocorrelation, we deleted duplicates of records. A total of 128 distribution sites, including 84 sites for N. tangutorum and 44 sites for N. sphaerocarpa, were used in the ENM analyses. These sites consist of 52 sampling sites in the study and 76 specimen records from the Chinese Virtual Herbarium (CVH, http://www.cvh.org.cn/). All bioclimatic layers with 19 bioclimatic variables at a resolution of 30 arc seconds were obtained from the WorldClim database (available at http://www.worldclim.org/; Hijmans et al., 2005). In order to avoid multivariate collinearity of environmental variables, which could lead to model over-fitting, we only retained distinct sets of variables that contributed most to models, and eliminated one variable per pair with correlations of r ≥ 0.8 according to Pearson correlation value (Wan et al., 2016). Seven bioclimatic variables, i.e. precipitation of warmest quarter (Bio18), temperature seasonality (Bio4), precipitation of driest quarter (Bio17), annual mean temperature (Bio1), mean temperature of coldest quarter (Bio11), annual precipitation (Bio12), and precipitation seasonality (Bio15), were selected for modeling the potential range of the two Nitraria species. We used 30 replicates with 80% of the geographical sites for training and 20% for testing. Graphics for distribution models in different periods of the two species were drawn using DIVA-GIS 7.5.

A total of 23 single nucleotide substitutions and 2 indels of 3 or 6 nucleotides within trnH-trnA sequence were detected in 633 individuals sampled from 52 populations and the length of the cpDNA fragment was 834 bp. However, only 8 nucleotide substitutions were detected in the atpH-atpI cpDNA fragment. The total length of combined the two chloroplast DNA sequences was 1794 bp, comprising 24 parsimony informative sites. Haplotype diversity (H _d) and nucleotide diversity (Pi) calculated with DnaSP 5.10 software were 0.7819 and 0.00094, respectively. A total of twenty-four different chlorotypes (H1-H24) were identified among all individuals sampled in the study ( Supplementary Table S1 ). However, the relationships among these chlorotypes have not been resolved by phylogenic analysis ( Supplementary Figure S1 ). Half of the populations sampled in this study were fixed for a single chlorotype, and the other half of the populations were polymorphic ( Table 1 , Figure 1A ). The most common chlorotypes were H4, H6 and 19 ( Figures 1A, B ). H4 widely occurred in most populations of N. tangutorum at high frequency and a few populations of N. sphareocarpa (13 and 32). H6 only occurred in some populations of N. tangutorum as a single chlorotype (24, 26, 27, 29 and 31) or at different frequencies (4, 15, 22, 23, 30 and 41). However, H19 was only fixed for most populations of N. sphaerocarpa as a single chlorotype (8, 10, 11, 16, 17, and 19), and for one population (32) at lower frequency. The other haplotypes occurred in a few populations at different frequencies. Especially, there were eight rare haplotypes (H2, H7, H9, H11, H12, H21, H22 and H24) which occurred only in one or two individuals of the two species, and the other two rare haplotypes, H13 and H18, were present only in two populations of N. tangutorum at low frequencies, respectively (H13: 39, 42; H18: 50, 51). In addition, H4 and H14 were shared by two species while all the other haplotypes identified in this study were fixed by only one of the two species ( Figure 1B ). The haplotype diversity (H _E) calculated for each population ranged from 0 to 0.758, with a mean value of 0.208 ( Table 1 ). Population 4, which had the greatest haplotype diversity (H _E = 0.758), was polymorphic for chlorotypes H1, H4, H5, and H6. Furthermore, the other four populations (21, 34, 42 and 50) had high haplotype diversities (H _E = 0.714, 0.600, 0.692 and 0.724, respectively) for three or four haplotypes. Although the two populations (5 and 6) contained high haplotype numbers, their values of H _E were low (H _E = 0.490, and 0.466, respectively) because of low frequencies of some chlorotypes occurring in these populations (e.g., H1, H3). Notably, population pairs of N. tangutorum and N. sphareocarpa at the same sampling site contained completely different chlorotype composition ( Figure 1A ; population pairs 11/12, 15/16, 18/19, 24/25, and 27/28), implying that a complete reproductive isolation has been established between the two related species.

The SAMOVA analysis revealed that the differentiation among groups (F _CT) of N. tangutorum reached a plateau when the value of K was 6 ( Supplementary Tables S2 , Supplementary S3 ; Supplementary Figure S2 ), and all populations sampled for this plant clustered into six groups. The group 6 included most of these populations (25 populations) which occurred everywhere in the whole distribution area of this species, while the group 5 consisted of 7 populations (22, 24, 26, 27, and 29-31) which occupied almost the same geographical region ( Figure 1A ). However, all the other groups (group 1-4) only included one or two populations. For the species N. sphaerocarpa, the SAMOVA divided all sampled populations into five groups (K = 5) ( Supplementary Tables S2 , Supplementary S3 ; Supplementary Figure S2 ). Except for groups 4 and 5 which contained most populations of the species (10 populations), the other groups only consisted of one or two populations. The populations of group 4 or 5 also occupied approximately the same geographical region ( Figure 1A ). The group structures based on the SAMOVA analysis were not congruent with geographical distributions of all populations of the two species. However, most populations in each close geographical region obviously contained homologous or unique chlorotype, for example, populations 33-38 in the eastern region of sampled area (Inner Mongolia) fixed unique H10, and populations 42-44 in North Xinjiang region fixed unique H15 ( Figure 1A ). Especially, almost all the populations of N. sphaerocarpa in the western range of Gansu province fixed the same chlorotype (H19).

To further reveal genetic structure of these populations both in the whole range and in different geographical regions, therefore, all the populations of N. tangutorum in this study were divided into six geographical groups (Gg1-Gg6) according to different geographical regions and climatic habitats ( Table 2 ). Gg1 consisted of populations 1-7 which mainly occurred in the Qinhai-Tibet Plateau (QTP) with an average altitude of more than 4000m. Gg2 included populations 9, 12, 15, 18, and 20, which were distributed in the west of Gansu province where belongs to extreme arid area with little rainfall of less than 100 mm. The range of Gg3 (22-24, 26, 27, and 29-31) lied in arid area where the rainfall is more than 100mm but less than 200mm. Gg4 populations (33-38) are largely distributed near the Yellow River where belongs to semi-arid regions. Gg5 and Gg6 populations mainly occurred in the northern and western regions of Xinjiang, respectively, which are geographically far from the populations of other geographical groups. However, the populations of N. sphaerocarpa were not divided according to geographical regions, because most populations of the species in a close geographical rang (in the west of Gansu province) fixed a single chlorotype and other geographical regions contained only a few sampled populations.

The PERMUT analysis indicated that the two species, N. tangutorum and N. sphaerocarpa, have low average gene diversities within populations (H _S = 0.249 and 0.099, respectively), while the total gene diversities across all populations were high (H _T = 0.654 and 0.758, respectively) ( Table 2 ). The two species also showed high level of differentiation among all the populations (G _ST = 0.620 and 0.870, N _ST = 0.624 and 0.932, respectively). The total gene diversity (H _T) and interpopulation differentiation (G _ST and N _ST) of N. sphaerocarpa were obviously higher than that of N. tangutorum. A permutation test detected a significant phylogeographic structure (N _ST > G _ST; P < 0.01) among populations of N. sphaerocarpa ( Table 2 ), while no clear phylogeographic structure was detected among populations of N. tangutorum. For the six geographical groups of N. tangutorum ( Table 2 : Gg1-Gg6), the total gene diversity of Gg1 in QTP (H _T = 0.818) was apparently higher than that of the other geographical groups, and only this geographical group showed a clear phylogeographic structure. The Mantel test also identified an obvious correlation between genetic distance and geographical distance for N. tangutorum (R = - 0.103, P = 0.063), while the correlation between the two distance matrixes was not clear for N. sphaerocarpa (R = 0.406, P = 0.019) ( Figure 2 ).

AMOVA analyses showed that approximately 58% of the total chloroplast DNA variations occurred among populations of N. tangutorum whereas about 87% of variation occurred among populations of N. sphaerocarpa, also indicating that the level of genetic differentiation among populations of N. sphaerocarpa was higher than that of N. tangutorum ( Table 3 ). For N. tangutorum and its six geographical groups, about 33% of variation occurred among groups (F _CT = 0.3287), revealing that there was a significant differentiation among these geographical groups. Actually, the BARRIER analysis also detected a few potential biogeographical boundaries among the populations of the two species. Two robust boundaries with high bootstrap values (> 95%) were identified among populations of N. tangutorum and N. sphaerocarpa, respectively ( Figures 1C, D ). AMOVA analysis also revealed that the molecular variation within populations of the four geographical groups (Gg2, 93.98%; Gg4, 85.76%; Gg5, 64.90%; and Gg6, 76.61%) were notably higher than that of among population within groups of N. tangutorum ( Table 3 ).

The mismatch distribution analysis using DnaSP 5.10 revealed that the distribution of pairwise differences for N. tangutorum populations displayed a smooth and unimodal curve ( Figure 3A ), implying this species experienced a sudden demographic expansion in the past. However, the distribution of pairwise differences of N. sphaerocarpa was a bimodal curve ( Figure 3B ). The higher P-values of Harpending’s raggedness index (RAG) and the sum of squared deviation (SSD) for N. tangutorum ( Table 4 ) further revealed that the species experienced a rapid range expansion in its whole geographical distribution. Furthermore, the obviously negative values of Tajima’s D and Fu’s Fs in the whole populations of N. tangutorum also supported the sudden demographic expansion of the species ( Table 4 ). The demographic expansion of N. tangutorum approximately happened between 7 and 21 Kya, according to the average mutation rates of cpDNA sequences in angiosperms (Wolfe et al., 1987), the two chloroplast DNA sequences length of 1794bp, and the generation time of 3 years. However, a rapid demographic expansion was not identified in the whole geographical distribution of the species N. sphaerocarpa according to statistics for neutrality tests and mismatch distribution analysis ( Table 4 ).

The distribution models of the two species in ENM obtained the high average AUC scores (0.875 and 0.953 for N. tangutorum and N. sphaerocarpa, respectively), based on the 30 replicates of the MAXENT runs ( Supplementary Figure S3 ). Estimates of relative contributions of the environmental variables according to the jackknife tests showed that seven environmental variables, including Bio1, Bio4, Bio11, Bio12, Bio15, Bio17, and Bio18, principally influenced the geographical ranges of N. tangutorum and N. sphaerocarpa ( Supplementary Table S4 ). However, contribution rates of these seven environmental factors to the distribution ranges were slightly different between the two species. Among these environmental variables used in the study, four variables, i.e. Bio1, Bio4, Bio17, and Bio18, played a primary role in determining the potential ranges of the two species, and their total contribution rate accounted for over 85% ( Supplementary Table S4 ). Based on the 128 sites dataset and the seven environmental factors above, the potential distribution ranges of the two species were modeled for the LIG, the LGM, the present day, and the future, respectively ( Figures 4A–F ). The results showed that the potential ranges of the two species on the basis of a high habitat suitability index (>0.50) obviously fluctuated during the LIG and LGM ( Figures 4A–D ; Supplementary Table S5 ), comparing to the present and future ranges ( Figures 4E, F ). Especially, their potential ranges contracted significantly during the LGM period ( Figures 4C, D ; Supplementary Table S5 ), and subsequently, expanded or recolonized during the present day ( Figures 4E, F ; Supplementary Table S5 ). Furthermore, the two species will experience different degrees of range contraction under the future climate scenario, and the range contraction of N. sphaerocarpa is significantly greater than that of N. tangutorum based on a higher habitat suitability index (>0.74) ( Figures 4G, H ; Supplementary Table S5 ). Notably, the potential niche maps based on the ENM showed that the two species occupied distinct fragmented habitats in different periods, implying that the two species had multiple geographically isolated refugia in northern China. In addition, N. tangutorum occupied a significantly broader potential distribution range than N. sphaerocarpa, indicating that N. tangutorum had a higher ecological plasticity or adaptability to different habitats than N. sphaerocarpa.

In the present study, we examined high total haplotype diversity (H _d = 0.7819) and nucleotide diversity (Pi=0.00094) of the two cpDNA sequences (trnH-psbA and atpH-atpI). A total of twenty four chlorotypes were identified from 52 populations of N. tangutorum and N. sphaerocarpa. Out of these haplotypes, there was a single dominant haplotype (H4) which widely occurred in the entire geographical distribution of the two species ( Figures 1A, B ), and a lot of rare haplotypes which were fixed by only one or two individuals of them ( Table 1 ; Figure 1A ). However, the bootstrap values of most clads in the phylogenetic tree were very low (<50%, Supplementary Figure S1 ), indicating that the relationships among these chlorotypes were not clear. The PERMUT analysis indicated that the two species had high levels of total genetic diversity ( Table 2 : N. tangutorum: H _T = 0.654; N. sphaerocarpa: H _T = 0.758). In addition, we also found that four N. tangutorum populations (4, 34, 42, and 50) and one N. sphaerocarpa population (21) had obviously higher haplotype diversities than other populations ( Table 1 ; Figure 1A ), and these populations scattered different geographical regions. For example, population 4 occurs on the QTP with a series of mountains and valleys, population 42 lies in northern Xinjiang, while population 34 is distributed in central Inner Mongolia. Notably, these populations with higher haplotype diversities are geographically far from each other, and occupy different geographical regions. These results together imply that the ancestor populations of the two species were widely distributed in northwestern China before the Quaternary (Hewitt, 2000). These ancient populations have been gradually isolated into different geographical or ecological groups owing to the habitat fragmentation and formation of geographical barriers caused by past geological events and/or climate oscillations (Abbott et al., 2000; Hewitt, 2004), and subsequently formed the current genetic diversity pattern (Li et al., 2018; Jia and Zhang, 2019; Hu et al., 2022). The aridification and desert expansion in northwestern China induced by uplifting of the QTP (Sun et al., 1998; Guo et al., 2002; Sun and Liu, 2006) probably accelerated the range fragmentation, geographical subdivision and diversification of these desert species in this region (Li et al., 2012; Meng et al., 2015; Zhang et al., 2017; Hu et al., 2022). The rare haplotypes identified in this study were highly likely to be randomly retained in fragmented populations due to genetic drift during the population contraction/expansion process of the two species (e.g. Avise, 2004; Chen et al., 2008).

Although the total gene diversity of N. sphaerocarpa was significantly higher than that of N. tangutorum, most populations of N. sphaerocarpa were only fixed for a single haplotype (e.g. populations 8, 10, 11, 16, 17, 19, 25; Figure 1A ). The higher G _ST value of P. sphaerocarpa than N. tangutorum indicated that the genetic differentiation among populations of N. sphaerocarpa was more severe than that of N. tangutorum ( Table 2 ). AMOVA analysis also showed that approximately 87% of the total genetic variations occurred among populations of N. sphaerocarpa, while only about 57% variations occurred among populations of N. tangutorum ( Table 3 ). This different genetic structure between the two related species with co-distributed ranges was probably triggered by gene flow between populations (e.g. Zink et al., 2001; Polihronakis and Caterino, 2010; Zhang et al., 2012). Habitat dependence/preference of species with similar dispersal ability probably affects the gene flow among their populations. Species with high ecological plasticity and wide ranges have more opportunities to survive and spread in adverse conditions than other species with narrow niches and small ranges (Michaux et al., 2005; Hamer and McDonnell, 2010; Lange et al., 2010; Morris-Pocock et al., 2010). N. tangutorum is widely distributed across large parts of northern China, covering almost all desert areas in northwestern China. Its habitats include shifting or semi-fixed sandy land, and gravel or salinized sandy land. However, N. sphaerocarpa is scattered in desert areas discontinuously, and only occupies some narrow desert regions in northwestern China, including Midwest of Inner Mongolia, Hexi corridor of Gansu province, and a few areas of Xinjiang. This specie prefers to grow on fixed or gravel sandy land. This habitat preferences and habitat fragmentation probably affected the gene flow between populations of N. sphaerocarpa, and subsequently further contribute to the genetic divergence of the species.

Notably, among the twenty-four chlorotypes detected in this study, only two chlorotypes (H4 and H14) were shared by N. tangutorum and N. sphaerocarpa ( Figure 1B ). According to species-specific chlorotype composition of two populations in the same site (e.g. population pairs: 11 and 12, 18 and 19), we inferred that the two related plants were reproductively isolated by some biological mechanisms. Therefore, these shared chlorotypes between the two species were probably derived from retention of ancestor polymorphisms, but not from hybridization or introgression. These retained ancestor haplotypes had probably experienced an incomplete lineage sorting before a complete reproductive barrier was established between the two species (Wendel and Doyle, 1998).

The climate change induced by uplift of the QTP since the Pliocene accelerated the desert formation and expansion in northwestern China (Guo et al., 2002), and further caused range fragmentation and subsequent diversification of plants in this region (Sun and Li, 2003; Al-Shehbaz et al., 2006). Nitraria is an ancient Tertiray relic taxon, and its ancestor populations were highly likely to widely occur in northern China before these geological events (Temirbayeva and Zhang, 2015; Zhang et al., 2015). Our results revealed that the two species, N. tangutorum and N. sphaerocarpa, probably experienced habitat fragmentation caused by geological and climatic changes (aridification, desert formation and expansion/contraction) (Li et al., 2012; Meng et al., 2015). The hypothesis was supported by no clear phylogenic relationship among haplotypes and geographical distributions of populations with high haplotype diversity (e.g. populations 4, 21, 34, 42, 50). These populations with high diversity were likely to be the putative glacial refugia for the species (Abbott et al., 2000). We also did not detect any significant phylogeographical structure in the entire range of N. tangutorum ( Table 2 ), implying that all these populations of the species probably originated from different ancestor populations. The Mantel test also showed that there is no significant correlation between population genetic differentiation (Fst) and geographical distance of N. tangutorum in its whole distributional ranges ( Figure 1A ; 2). Therefore, it is highly likely that the climatic oscillations and associated environmental changes (e.g. desert formation and expansion/contraction) in the Quaternary accelerated the range fragmentation and population isolation, and subsequently providing chances for allopatric differentiation within species induced by selection and/or genetic drift (Meng et al., 2015). This inference has been further confirmed by a large number of rare haplotypes identified in the study. Actually, we also found some distinctive chlorotypes in different geographical regions, e.g. H15 in northern Xinjiang region, H16, H17 and H18 in western Xinjiang region, and H3 in Qinghai (QTP region) ( Figure 1A ; Table 1 ). These unique chlorotypes were also likely to result from genetic drift and randomly retained in fragmented geographical distributions.

The population groups divided by SAMOVA were not congruent with geographical subareas of the two species, for example, the group 6 of N. tangutorum included 25 populations which occurred everywhere in the entire ranges of this species ( Figure 1A ; Supplementary Table S3 ). However, we did find some geographically close populations which obviously contained some homologous or unique chlorotypes, such as populations 33-38 and 22-31, indicating that phylogeographical structure existed among populations in a few subareas of N. tangutorum ( Table 2 ). In the whole ranges of N. sphaerocarpa, we examined the significant phylogeographical structure (N _ST > G _ST, P < 0.01) and high interpopulation differentiation (G _ST = 0.870), which may be caused by limited gene flow among populations affected by habitat preference to narrow gravel sandy land. Actually, in the subarea of Hexi corridor both in Gansu and central Inner Mongolia, we discovered obvious lineage genetic differentiation among geographical populations of N. tangutorum and N. sphaerocarpa (21-38). A few genetic barriers based on Monmonier’s maximum-difference algorithm have also been detected among these populations in the subarea of the two species ( Figures 1C, D ). This genetic pattern of the two Nitraria species is partly congruent with the other two typical desert taxa we surveyed using cpDNA variations (Yu et al., 2013; Zhang et al., 2017). In addition, the geographical region Gg7 (populations 1-7) of N. tangutorum distributed in the fringe of QTP contained three distinctive chlorotypes (H1, H3, and H5) and one rare chlorotype (H2), and had the highest total gene diversity (H _T = 0.818) and interpopulation differentiation (G _ST = 0.570), compared to other geograpihical groups ( Table 2 ). A significant phylogeographical structure was also identified in this subarea group. All these results implied that a series of high mountains and deep valleys in this region played an important role in accelerating the formation of regional genetic differentiation and phylogeographical structure of the species (Dutech et al., 2004; Meng et al., 2015).

A few studies have revealed that recent population expansion probably results in shallow genetic divergence (Conroy and Cook, 2000; Hewitt, 2000; 2004; Polihronakis and Caterino, 2010). The results obtained by PERMUT revealed that genetic differentiation among populations in some subareas of N. tangutorum (e.g. Gg2, Gg4, Gg6; Table 2 ) were very low, and most of variations occurred within populations ( Table 3 ), implying that these populations in the subareas probably experienced recent regional expansion in their limited geographical regions. This scenario was further supported by mismatch analysis results ( Table 4 ; Figure 3A ). For N. sphaerocarpa species, we did not find a rapid demographic expansion in its whole geographical distribution ( Table 4 ; Figure 3B ), which was probably related to its habitat preference to narrow gravel sandy land. However, the ENM results revealed that the two species experienced apparent expansion and contraction during the LIG and LGM, respectively ( Figures 4A–D ; Supplementary Table S5 ), and their ranges have obviously expanded again in the present ( Figures 4E, F ; Supplementary Table S5 ). The mismatch analysis based on a range of possible cpDNA mutation rates indicated that the demographic expansion of N. tangutorum occurred approximately between 21 and 7 Kya before present, which was partly congruent with the ages of LGM (ca. 21 - 18 Kya before present). Numerous studies revealed that the deserts in central Asia, involving northwestern China, obviously enlarged during the Quaternary glacial ages (e.g. Bush et al., 2004; Sun and Liu, 2006; Qian et al., 2016). The range expansion of N. tangutorum probably followed the enlargement of the deserts during the ice ages, especially the LGM. Therefore, the recolonization of the species may have occurred after the LGM (ca. 18 - 7 Kya). This inference is basically consistent with the results of ENM and mismatch analysis. In addition, the jackknife tests in ENM revealed that two environmental variables, i.e. Bio18 and Bio4, played vital roles in modeling the potential ranges of the two species ( Supplementary Table S4 ), indicating that precipitation and temperature are the two key factors for the population survival and expansion of the two species. Some studies showed that anthropogenically induced climate change within the last decade is causing shifts in the distribution ranges of many species (Walther et al., 2002; Parmesan and Yohe, 2003; Parmesan, 2006), and these range shifts are likely to continue because of carbon emissions and global climate warming. Therefore, the global warming and aridification induced by human activities will also probably produce significant effects on potential distributions of N. tangutorum and N. sphaerocarpa. The ENM results in this study indicated that the suitable ranges of the two species will contract substantially in the future climate scenario ( Figures 4G, H ; Supplementary Table S5 ), when comparing to the present ( Figures 4E, F ). Moreover, the contraction degrees of distribution ranges are obviously different between the two species.

These different demographic patterns of the two species in responses to past climate changes are probably correlated with their contrasting habitat preferences and ecological plasticity (Michaux et al., 2005; Hamer and McDonnell, 2010; Lange et al., 2010; Morris-Pocock et al., 2010). In general, phenotypic plasticity can produce morphologies adapted to local conditions, and is very beneficial for the survival of plants in heterogeneous environments (Puijalon et al., 2008). Morphogenetic differences ensure resilience of species to new and changing ecological conditions (Schoelynck et al., 2015). N. tangutorum and N. sphaerocarpa had contrasting habitat preferences and different phenotypic plasticity (Pan et al., 1999). Furthermore, N. tangutorum had higher average gene diversity within populations (H _s) than N. sphaerocarpa ( Table 2 ), implying that the former had higher adaptability to heterogeneous environments than the latter.

Providing that the two species extensively occurred in the Northwest of China before the Quaternary, their current populations might have originated from a few separated glacial refugia following their habitat fragmentation induced by climate oscillations in the Quaternary (Chen et al., 2008; Li et al., 2012). Firstly, the range fluctuation of one species generally remains two evident genetic imprints: a wide distribution of a dominant haplotype, and multiple rare haplotypes (Hewitt, 1996, 2000; Comes and Kadereit, 1998). In the present study, we identified a single dominant haplotype (H4) which widely occurred in the entire geographical ranges of the two species ( Figures 1A, B ), and a lot of rare haplotypes which were fixed by only one or two individuals of them ( Table 1 ; Figure 1A ). Furthermore, we also found that some different geographical regions fixed a few distinctive chlorotypes, e. g. Gg1 fixed H3, Gg3 fixed H6, Gg5 fixed H15, and Gg6 fixed H16, H17 and H18, indicating that the species have experienced regional-scale range expansion/contraction in these geographical regions. Secondly, glacial refugium areas generally contained most of haplotypes and have high levels of genetic diversity (Abbott et al., 2000), while new recolonized areas usually have low haplotype diversities owing to founder effects (e.g., Petit et al., 2003; Stewart et al., 2010), and haplotypes should be decreased gradually from refugium (Hewitt, 2000; Heuertz et al., 2004). In the present study, we found four N. tangutorum core populations (4, 34, 42, and 50) and one N. sphaerocarpa core population (21) which had obviously high levels of genetic diversity ( Table 1 ; Figure 1A ), and these core populations are scattered across different geographical regions ( Figure 1A ). Moreover, we also found that other populations near these core populations had lower levels of genetic diversity and a few chlorotypes originated from their core populations. This genetic structure implies that the two species are highly likely to have independent refugia in these separated geographical regions, and they have experienced regional demographical expansion within different geographical regions.

Topographically heterogeneous areas are likely to act as refugia for species because they facilitate survival during regional climatic stress due to availability of a range of microenvironments (Byrne et al., 2022). Generally, mountainous regions can provide a few microhabitats for species survival when they facing adverse conditions or environments, while flat regions with few major geomorphologic features usually have a poorer refugial capacity (Trew and Maclean, 2021). Therefore, topographically complex regions usually retained higher levels of genetic diversity (Garrick, 2011; Byrne et al., 2017). In the present study, the four putative refugia (populations 4, 34, 42, 50) for N. tangutorum had high haplotype diversities, and H _E ranged from 0.600 to 0.758 ( Table 1 ). Apart from population 34 with a relatively lower level of haplotype diversity (0.600), the other three populations harbored diversity values with no significant difference among them, implying that these putative refugia had similar ecological capacities for preserving genetic diversity of the species. Furthermore, these populations with high haplotype diversities (putative refugia) are all located near some mountains ( Figure 1A ), for example, populations 34, 42 and 50, are located near the Helan Mountain, Altai Mountains, and Kunlun Mountains, respectively, although the species prefers to grow on flat sandy land in deserts.

Our results together suggest that these current populations of N. tangutorum and N. sphaerocarpa originated from a few separated glacial refugia following their habitat fragmentation in the Quarternary, and these populations have experienced regional range expansion after the LGM. Interestingly, we also found that these putative refugia in the study are mostly located near some rivers. For example, population 21 is located near the Shiyang river originating from the Qilian Mountain, 34 is located near the Yellow River, 42 is located near the Ertiz river originating from the Altai Mountains, and 50 is located near the Yarkant river originating from Kunlun Mountain. The jackknife test of the ENM revealed that Precipitation of Warmest Quarter (Bio18) is the most important factor which influenced the geographical ranges of the two species ( Supplementary Table S4 ). Therefore, these seasonal or permanent rivers probably provided the possibility for the survival of the two species in different refugia when the glacial climate became dry and cold (Xu et al., 2010; Ge et al., 2011; Meng and Zhang, 2013; Meng et al., 2015).

In the present study, we only used two maternally inherited chloroplast DNA fragments to highlight the demographic history and genetic structure of closely related species with overlapping ranges in response to past geological events and climatic oscillations. Our results together suggest that the two typical desert species, N. tanguotorum and N. sphaerocarpa, responded to the past climate fluctuations in different ways due to special habitat preferences. However, phylogenetic relationships among populations in the whole geographical distribution have not been well resolved. Intraspecific gene flow among and within sub-geographical regions are not clear. Further studies using more plastid DNA fragments and nuclear genes (e. g. unlinked low-copy nuclear genes) are needed to reveal the demographic history of these desert plants with overlapping ranges.