The role of East Asian monsoon system in shaping population divergence and dynamics of a constructive desert shrub Reaumuria soongarica

Both of the uplift of Qinghai-Tibet Plateau (QTP) and the development of East Asian monsoon system (EAMS) could have comprehensively impacted the formation and evolution of Arid Central Asia (ACA). To understand how desert plants endemic to ACA responded to these two factors, we profiled the historical population dynamics and distribution range shift of a constructive desert shrub Reaumuria soongarica (Tamaricaceae) based on species wide investigation of sequence variation of chloroplast DNA and nuclear ribosomal ITS. Phylogenetic analysis uncovered a deep divergence occurring at ca. 2.96 Mya between the western and eastern lineages of R. soongarica, and ecological niche modeling analysis strongly supported that the monsoonal climate could have fragmented its habitats in both glacial and interglacial periods and impelled its intraspecific divergence. Additionally, the population from the east monsoonal zone expanded rapidly, suggesting that the local monsoonal climate significantly impacted its population dynamics. The isolation by distance tests supported strong maternal gene flow along the direction of the East Asian winter monsoon, whose intensification induced the genetic admixture along the latitudinal populations of R. soongarica. Our results presented a new case that the development of EAMS had prominently impacted the intraspecific divergence and population dynamics of this desert plant.


Results
Genealogies and geographical distributions of haplotypes. Five cpDNA fragments and one nrITS fragments for a total of 272 in-group and 3 out-group individuals were sequenced. The sequences obtained from five cpDNA fragments were concatenated into a total length of 3,836 bps, with 41 chlorotypes defined by 44 polymorphic sites, including two indels (Supplementary Table S1). A total of 65 ribotypes (59 ribotypes from in-groups and 6 ribotypes from out-groups) were detected from the nrITS fragment (662 bps). Accordingly, the genealogical topology of chlorotypes and ribotypes were both clustered into two major lineages from the out-groups, and this clustering was strongly supported by the maximum likelihood (ML) trees ( Fig. 1; Supplementary Fig. S1a). All chlorotypes and ribotypes from the two lineages were then projected to the population location on the geographical map (Fig. 2a,b). Most of the western populations harbored the chlorotypes and ribotypes from the western lineages, whereas most of the eastern populations were dominated by the chlorotypes and ribotypes from the eastern lineages.
Scientific RepoRts | 5:15823 | DOi: 10.1038/srep15823 Chlorotypes C1 and C9, which coalesced to chlorotype C39 connecting with the out-group, could be the most recent common ancestral chlorotypes for the eastern and western clades, respectively ( Fig. 1a; Supplementary Fig. S1a). Chlorotype C1 was found in most of the populations from the eastern group BJ-TD (grouping refers to Materials and Methods), whereas C9 was dominant in the western group TaD. Interestingly, chlorotypes C1 and C9 were both found in the populations from the QB and KD-GG groups. Four of the five populations from the group QB and one marginal population HSW from the group BJ-TD near the QB harbored both the western and eastern chlorotypes ( Fig. 2a; Supplementary  Table S2), which suggested that the maternal genetic admixture occurred along the populations of the QB "corridor". In the KD-GG region, however, one population SSG was fixed with the western chlorotype C9, while the other three populations were fixed with the eastern chlorotypes ( Fig. 2a; Supplementary  Table S2), indicating that strong genetic structures were characteristics of chloroplast DNA in the group KD-GG.
Similar to the chlorotypes, most of the ribotypes could be divided into two lineages, except for the ancestral ribotype R4 being shared among all populations from the eastern and western regions (Fig. 1b), which showed incomplete divergence and/or imcomplete lineage sorting on nrITS within this species. However, we failed to detect the genetic admixture occurred on the nuclear genome in the group QB, as none of western ribotypes but ancestral ribotype R4 and the eastern ribotypes were found in all the populations from the group QB (Fig. 2b). Furthermore, we could not detect genetic structure on the nuclear genome within the group KD-GG, as the population SSG and the other three populations from the KD-GG region were all fixed with the eastern ribotypes and R4.
As noted above, all chlorotypes and ribotypes found in the group GuD were derived from both lineages (Fig. 1a,b). Interestingly, in this region, two populations (SW and FK) were dominated by the  59 nuclear ribotypes (R1-R59) in R. soongarica from the five groups in ACA, e.g. Taklamakan Desert (TaD), Gurbantunggut Desert (GuD), Qaidam Basin (QB), Kumtag Desert and Gashun Gobi (KD-GG), and Badain Jaran and Tengger deserts (BJ-TD). The size of the circles corresponds to the frequency of each haplotype. Small gray circles indicate unsampled or extinct haplotypes. Each number near the line represents the mutational steps interconnecting two haplotypes, and only the steps over two mutations are listed. There were three mutation steps between C39 and C9, as each mutation represents 5. 35  western chlorotypes and eastern ribotypes, whereas the other population, HSS, was dominated by the eastern chlorotypes and western ribotypes (Fig. 2a,b). These results indicated that R. soongarica recently colonized into the GuD region from two directions.
Divergence and diversification of chlorotype lineages. Due to insufficiency of pollen or fossil records of Reaumuria, we carefully chose an average mutation rate of 4.87 × 10 -10 substitutions per site per year (s/s/y) for chloroplast genome, which is suitable for desert shrub taxa 27 , to estimate the divergence time between the genealogical lineages. According to the BEAST analyses of cpDNA based on the above mutation rate, the intraspecific divergence time between the eastern and western lineages can be dated to ca. 2.96 Mya (95% HPD, 1.45-5.50 Mya, point "c" in Fig. 3), corresponding to the last peak of higher temperature during the Pliocene 28 and/or the prominent uplift of the QTP 10 . The eastern lineage coalesced at ca. 1.72 Mya (95% HPD, 0.83-3.13 Mya point "d" in Fig. 3), considerably earlier than the western lineage at ca. 1.26 Mya (95% HPD, 0.52-2.46 Mya point "e" in Fig. 3), suggesting the asynchronous lineage diversifications of R. soongarica between the two regions (Badain Jaran and Tengger deserts and Taklamakan desert) before the Middle Pleistocene. Interestingly, in the group GuD, the presence of (a) cpDNA chlorotype distribution; (b) ITS ribotype distribution. Population locations are presented as small black spots with white margins; slice sizes are proportional to the frequency of the haplotypes (eastern or western clade). The colors of the pie diagrams correspond to the clades of the network genealogies in Fig. 1. "Old" refers to C26, C39 and C24. The gene flow along the East Asian winter monsoon (EAWM) and East Asian summer monsoon (EASM) directions are displayed with red arrow lines in (a). In addition, genetic barriers are superimposed on the map between populations using a gray bold line, with the width of the barrier lines representing the bootstrap values. These figures were drawn by H.X.Y. using CorelDraw X6 (Corel Corporation, Ottawa, Canada) after being originated from the software packages Diva-GIS version 7.5.0 (http://www.diva-gis.org/) loading with the Chinese map. the west-like chlorotype C12 and the east-like C10 and C11 could be dated to later than ca. 0.37 Mya (point "f " in Fig. 3), indicating that the bi-directional colonization of R. soongarica could have occurred in the Gurbantunggut desert during the Late Pleistocene.
Genetic isolation and gene flow. The overall genetic diversity (H T = 0.934) of R. soongarica on cpDNA was considerably higher than the average diversity within populations (H S = 0.527), which was concordant with the significant phylogeographical structure found on cpDNA (N ST = 0.612 with P = 0.04, G ST = 0.436 with P = 0.04, Table 1) and on nrITS (N ST = 0.386 with P = 0.03, G ST = 0.279 with P = 0.02; Supplementary Table S3). A hierarchical AMOVA indicated that the molecular variation among five groups (48.04%) was nearly two-fold higher than that within populations (24.74%, Table 1). In addition, strong significant pairwise genetic differentiation (Φ ST ) was found between the group TaD and the other four groups (from 0.675 to 0.819, Table 2), which supported the geographical isolation of group TaD. This finding suggested that the gene flow from other groups couldn't have affected the population diversity within group TaD. In contrast, negligible pairwise Φ ST values were found between the three latitudinal groups (Φ ST (QB, KD-GG) = 0.005, Φ ST (QB, GuD) = 0.084, Φ ST (GuD, KD-GG) = 0.108; Table 2). This finding supported the hypothesis that frequent latitudinal gene flow on chloroplast genome occurred among these three groups during the evolutionary history of R. soongarica.
Significant population structure was also found within the groups GuD and KD-GG (Hs = 0.083 and 0.179, F ST = 0.993 and 0.961, respectively; Table 1). In contrast, the highest genetic diversity within populations (Hs = 0.629) and the lowest genetic differentiation (F ST = 0.191) were found within the group QB, supporting the frequent maternal gene flow among the populations within this group ( Fig. 2a; Table 1).
To determine whether genetic differentiation was significantly related to physical distance (isolation by distance), Mantel tests were conducted along the latitudinal and longitudinal directions both on maternally inherited cpDNA and bi-parentally inherited nrITS. Significant genetic isolation by physical distance was found along both latitudinal and longitudinal directions on nrITS (r = 0.430, P < 0.001; r = 0.794, p < 0.001, respectively), and along the longitudinal direction on cpDNA (r = 0.609, p < 0.001; Table 3). However, no significant genetic isolation was detected along the latitudinal direction on cpDNA (r = 0.262, P = 0.110; Table 3). This finding implied that the EAWM played a crucial role in enhancing the maternal gene flow, which minimized the genetic differentiation between populations along the latitudinal direction. (a-f) represent the coalescence time of each major clade, and the length of the blue bars represents 95% highest posterior density. The capital letters at the end of each branch representing the original groups, e.g. Taklamakan Desert (T), Gurbantunggut Desert (G), Qaidam Basin (Q), Kumtag Desert and Gashun Gobi (K), and Badain Jaran and Tengger Desert (B).

Regional population dynamics.
A marginally significant negative value of Tajima's D (-1.27, Table 1) and MDA results ( Supplementary Fig. S2) implied that the overall population of R. soongarica underwent a complex expansion of its population structure. Significant negative values of Tajima's D (− 1.55, P = 0.035) and Fu's Fs (− 12.13, P < 0.001) strongly confirmed that the eastern group BJ-TD underwent a population expansion, whereas the western group TaD remained comparatively stable, as only Tajima's D was significant negative (− 1.72, P = 0.021; Table 1). The skyline plots with MIGRATE-N analysis showed that the population size of the overall R. soongarica populations had been increasing significantly ( Fig. 4; Supplementary Fig. S3a). In the separate groups, the population size of group BJ-TD increasing more rapidly and reached a maximum at 80 Ka, whereas group TaD exhibited a slower population growth relatively recently ( Fig. 4; Supplementary Fig. S4b,d). However, population decreasing was found in the three other groups with genetic admixture/structure: GuD, QB and KD-GG ( Fig. 4; Supplementary Fig. S4c,e,f).
Distribution range shifts in the past. The ENM simulations based on Maxent supported the present distribution range of R. soongarica, which covers nearly all the desert regions across ACA (Fig. 5a). Our modeling data indicated that Bio19 and Bio11, namely, the precipitation and average temperature during the coldest season, together represented more than 60% of the niche contribution and over 80% of the permutation importance (Supplementary Table S4). These results indicated that the EAWM markedly influenced the distribution range of R. soongarica. Compared with the present distribution, the area of suitable habitat determining the range of R. soongarica was extremely contracted and fragmented by the QTP during LGM period, and this tendency further strengthened during LIG period (Fig. 5b,c). However, the local distribution range fluctuated asynchronously: the western region experienced range expansion, and the eastern region underwent more frequent range shifts, following LIG period (Fig. 5b,c). These results indicated that the paleoclimate change could have further driven the habitat fragmentation of R. soongarica and that the EAMS could have markedly affected vegetation distribution in different ways in the two local regions.

Discussion
Based on thorough investigation on genetic structure and distribution range shifts under monsoonal climate scenarios, we intend to understand the mechanism of population divergence and population dynamics of R. soongarica, a Tertiary relic constructive zonal shrub species widely distributed across all the desert ecosystems in ACA. As the formation and evolution of ACA had been affected by the development of EAMS since Pliocene 7 , the present work supported an idea model to discuss the role of monsoonal climatic changes in shaping the habitat fragmentation which triggered the intraspecific divergence/speciation of the dominant species of ACA. The cpDNA genealogy network and phyML tree ( Fig. 1a; Supplementary Fig. S1a) suggested that R. soongarica diverged into two contrasting lineages. This finding is consistent with our previous study 21 . The chlorotypes from the western lineage were mainly distributed in the group TaD, whereas those from the eastern lineage were dominant in the group BJ-TD (Fig. 2a). Marked genetic differentiation between these two groups was found in both nrITS and cpDNA (Table 2). Due to the effective population size of the biparentally inherited nuclear genome is twofold larger than that of the uniparentally inherited chloroplast genome, we found incomplete divergence on nuclear ribosome DNA. The most recent common ancestral ribotype, R4, was shared among all groups, excluding GuD (Fig. 1b). However, we detected significant genetic barriers on nrITS between the group TaD and the other groups (Fig. 2b). This marked genetic differentiation between the eastern and western populations in ACA is not specific to R. soongarica. It has also been found in some sympatric plant species such as Nitraria sphaerocarpa 18 and   To further clarify the cause of the intraspecific divergence between the two groups, molecular dating was conducted with a relaxed clock method via Bayesian simulation based on substitution rate of 4.87 × 10 −10 s/s/y for the chloroplast genome 27 . The divergence between the two out-group genera could be dated to ca. 43 Mya (95% HPD: 16.63-74.68 Mya, point "b" in Fig. 3b). This result did not conflict with the earliest fossil record of the genus Tamarix 29 . The crown age of genus Reaumuria and the out-groups could be dated to ca. 66 Mya (95% HPD: 111-27 Mya, point "a" in Fig. 3b). It is consistent with a previous report 30 , and corresponding to the period of Paleocene-Eocene Thermal Maximum (PETM, Fig. 3a), when a sharp increase in concentration of atmospheric CO 2 caused the global surface temperature to increase by approximately 5-6 °C 28 . The intraspecific divergence of R. soongarica between the eastern and western lineages could be dated to ca. 2.96 Mya (95% HPD: 1.45-5.5 Mya, point "c" in Fig. 3b). This time point is close to the strong tectonic uplift of QTP, especially, a rapid uplift of Qaidam Basin (the northeastern QTP) during the late Pliocene 7,14 . Such a rapid uplift is consistent with our molecular dating of the intraspecific divergence of R. soongarica to a period prior to the beginning of the dramatic glacial and interglacial oscillations of the global climate (Fig. 6a).
However, the abrupt uplift of the QTP had also caused the intensification of EAMS at ca. 3.6 to 2.6 Mya. Actually, based on ENM simulation of the potential distribution range of R. soongarica during the LGM (EAWM dominated) and LIG (EASM dominated), both the glacial and interglacial climates could have greatly contracted the habitat of this species and finally fragmented its distribution range (Fig. 5b,c). On the other hand, the present distribution of R. soongarica was not fragmented at all (Fig. 5a), indicating that the two lineages of R. soongarica should have contacted easily and the two reciprocal habitats could have merged together under the moderate climatic conditions. Therefore, we could not arbitrarily conclude that the geological event was the only condition that triggered the intraspecific divergence of R. soongarica. Based on the considerations above, we proposed that once the isolation occurred as a result of geological cause, the cycles of the monsoonal climatic oscillations was most likely to play a key role in the habitat fragmentation and intraspecific divergence of R. soongarica. After R. soongarica diverged into two regional lineages, further diversification of two lineages began ca. 1.72 Mya (Eastern clade, point "d" in Fig. 3b) and 1.26 Mya (Western clade, point "e" in Fig. 3b), which are fairly concordant with the two main turning points of the EAMS 31,32 . Stratigraphic data from loess deposits and NW-Pacific ocean deposits suggested that the EAWM steadily strengthened during Pleistocene, whereas the EASM underwent stepwise weakening until ca. 1.8 Mya, when the EASM was increasingly restrengthened. Considering the more fragmented habitat of R. soongarica in the interglacial period such as LIG (Fig. 5c), we speculated that the strengthening EASM of the warming stage ca.1.8 Mya could have fragmented the habitat of R. soongarica and accelerated the diversification of the eastern lineage. However, the western lineage  Fig. 2 as follows: the brown spot represents the ancestor of R. soongarica, the blue pie represents the haplotyps from the eastern clade, and the green pie represents ones from the western clade. The shadow represents the expanding arid regions responding to the habitat changes of R. soongarica during different geological periods. These figures were drawn by H.X.Y. using CorelDraw X6 (Corel Corporation, Ottawa, Canada).
could not start diversification until ca. 1.3 Mya, when the restrengthened EAWM effectively overwhelmed the EASM and induced the further aridification of ACA 32 . Therefore, we suspected that the EAWM played a crucial role in the diversification of the western lineage of R. soongarica.
Since the asynchronous diversification happened between the two lineages of R. soongarica, could the two regional groups have undergone the similar population dynamics? In fact, from the skyline plots obtained with MIGRATE-N, we found a significant increase of the population size in the eastern groups of R. soongarica (Fig. 4), concordant to significant negative values of the Tajima's D and Fu's Fs and star-like shape of the genealogy network ( Fig. 1; Table 1). Comparatively, the population size of the western group TaD was much more stable ( Fig. 4; Supplementary Fig. S4d). Considering the distribution range of R. soongarica expanding slightly in the western region from LIG to the present based on ENM analysis (Fig. 5b,c), the stable population dynamic of the western group TaD could be related to the geographical isolation caused by the Tianshan-Kunlun-Arjin Mountains and equable climate during the late Pleistocene 15,33 .
Given that the population size of the western group remained stable while its distribution range expanded, and that the eastern group underwent a remarkable population growth while its distribution range shifted significantly (Figs 4 and 5b,c), the population expansion of the eastern group could be attributed to the frequent range shifts rather than to habitat range expansion. The frequent range shifts during the Pleistocene exerted strong selective forces that facilitated the accumulation of new mutants in the local populations, and consequently induced the effective population to grow more rapidly in the eastern climate-sensitive group than in the western climate-stable group. Previous studies have suggested that during the Pleistocene climate oscillation, population growth was quite common for large-sized species, such as conifers endemic to the climatic-sensitive region, e.g., the southeast QTP 4,34 and other plants 16 and animals 35 endemic to ACA and adjacent regions. Unfortunately, these studies failed to interpret the potential effects of local climate change on the speed and scale of population growth. Here, based on the comparison between the regional groups of the population dynamics and range shifts of R. soongarica, we presented a novel example that the populations in the monsoonal climate zone tended to grow more rapidly than those far from it.
Additionally, an interesting geographical pattern of chlorotypes of R. soongarica showed that frequent genetic admixture of the two lineages occurred in most populations within the QB "corridor" (Fig. 2a). On the other hand, we detected negligible pairwise genetic differentiation (Φ ST ) among the groups GuD, Ku-GG and QB, which suggested that the latitudinal gene flow had been markedly enhanced (Table 2). To further verify the effects of the EAMS on gene flow among the groups of R. soongarica, we conducted Mantel tests on maternally inherited cpDNA and bi-parentally inherited nrITS along different monsoon directions. Our data strongly supported that the EAWM decisively accelerated maternal gene flow along the latitudinal direction (Table 3). Due to the hairy and light characteristics for the mature seeds of R. soongarica, it is reasonable that the seeds of R. soongarica could be easily dispersed by the EAWM during autumn, accelerating maternal gene flow across the QB "corridor. " Gene flow can minimize intra-and inter-specific genetic differentiation and induce the accumulation of heterozygosity, facilitating the colonization of new habitats and promoting adaptation to climate change 36,37 . Many well-adapted plant species such as forest trees would utilize wind systems for seeds or pollen dispersal to increase gene flow between diverged populations [38][39][40] . By dating the historical maternal gene flow occurred among the groups of R. soongarica, we concluded that the EAWM intensified at least after the origin of the chlorotype C9 in the western group (ca. 1.4 Mya, more reasonable after 1.3 Mya, Figs 1a and 6b). In fact, many stratigraphic investigations across the Loess Plateau have recorded that the EAWM intensified after 1.3 Mya 7,8,32 . Up to date, empirical data are lacking except our new evidence from the desert plant R. soongarica on the scale and direction of the gene flow induced by the EAWM.
More interestingly, phylogenetic data of chlorotypes suggested that R. soongarica has recently (since ca. 0.37 Mya, point "f " in Fig. 3b) colonized into the GuD region, which was concordant to that the Gurbantunggut Desert expanded due to the intensification of the EAWM ca. 0.4 Mya 41 . This phylogeographic pattern indicates that the aridification of ACA during this time had expanded to the Gurbantunggut Desert region, which had supplied new habitats for the xerophyte R. soongarica (Fig. 6c). The stepwise colonization of plant species has also been found in an alpine conifer species endemic to the QTP 42,43 . In this study, we provided the first evidence for the stepwise aridification processes of ACA from a molecular phylogeographic perspective and confirmed that the formation time of Gurbantunggut Desert was later than the other deserts, e.g., Badain Jaran and Tengger deserts and Taklimakan Desert in ACA 44. The onset of the desertification of ACA has been debated for several decades, with the time estimates ranging from ca. 5.3 Mya to 0.85 Mya 45-49 . The regional asynchronous diversification of the constructive shrub R. soongarica indicated that the desertification of the eastern Badain Jaran and Tengger deserts regions could have occurred at 1.72 Mya, the western Taklimakan Desert could have formed ca. 1.26 Mya, and the northern Gurbantunggut Desert could have formed ca. 0.4 Mya. Of course, as the corresponding 95% CIs of dating were still wide, more evidence from fossil and molecular records are needed to clarify the accurate formation time of the deserts across ACA.
In summary, our study supported a unique case focused on evaluating the roles of EAMS in impacting the evolutionary history of constructive shrub species endemic to ACA, by triggering the intraspecific divergence between the longitudinal populations and enhancing gene flow among the latitudinal Scientific RepoRts | 5:15823 | DOi: 10.1038/srep15823 populations. Furthermore, this work supplied new molecular evidence for understanding the spatial and temporal evolution of ACA from a phylogeographic view point.

Methods
Population sampling. A total of 272 individuals were sampled from 34 populations across the ACA geographical ranges (Supplementary Fig. S5; Table S5). Six of these populations were collected in Taklimakan Desert (TaD) region, three in Gurbantunggut Desert (GuD) region, five in Qaidam Basin (QB) region, four in Kumtag Desert and Gashun Gobi (KD-GG) region, and 16 in Hexi corridor and Badain Jaran -Tengger deserts (BJ-TD) regions. Additionally, three individuals of Tamarix amplexicaulis and Myricaria pulcherrima were collected as outgroups to Reaumuria 24 . In each population, samples of fresh leaves were collected from eight individuals at least 30 meters apart to reduce sampling bias and dried by silica gel. The information on the location of each population, including latitude, longitude and altitude, were recorded with an e-Trex Global Positioning System (Garmin, Taiwan; Supplementary Table  S5 Table S6). Due to the maximum polymorphic sites, five cpDNA fragments (petB-petD, trnS2-trnG2, atpH-atpI, ndhA and petN-psbM) and one nrITS fragment were selected for full phylogeographic surveys (Supplementary Table S7 Population genetic structure. Population gene diversity (H S , H T ) and between-population divergence (G ST , N ST ) were estimated for each region based on the cpDNA using the program PERMUT with 1,000 permutation tests (available at http://www.pierroton.inra. fr/genetics/labo/Software/PermutCpSSR) 54 . In addition, spatial genetic structures of chlorotypes and ribotypes were analyzed with SAMOVA Version 1.0 55 , using a simulated annealing procedure to identify the number of population groups (K) that were geographically homogeneous and maximally differentiated from each other. According to the various geographical locations and SAMOVA results, all populations of R. soongarica were divided into five groups as follows: QB, BJ-TD, TaD, GuD and KD-GG. Notably, the populations from the KD-GG group were considered one group due to the geological "junction" location of the GuD, TaD, QB and BJ-TD groups. The populations collected from the Hexi Corridor and from the Badain Jaran-Tengger deserts were combined as the BJ-TD group according to the analysis of the initial polymorphism survey and geographical locations (Supplementary Fig. S5). Further genetic differentiation among groups (F CT ), among populations within groups (F SC ), and within populations (F ST ) was calculated with a hierarchical analyses of molecular variance (AMOVA) using Arlequin version 3.11 56 , with significance tests of variance components and differentiation indices with 1,000 permutations. Pairwise differentiations (Φ ST ) between all pairs of groups were also estimated with 1,000 permutations.
We used the software BARRIER version 2.2 57 to identify the areas of genetic discontinuity or barriers to gene flow across the sampled species' ranges. One hundred bootstrapped datasets were generated by re-sampling individuals within populations with replacement using PopTools version 3.2.5 (available at http://www.poptools.org/). Under ideal circumstance, the pairwise genetic differentiation between the populations could not reject the Mantel test supporting isolation by distance 58 . To test the significance of the impact of wind systems on maternal gene flow limited by physical distance, Mantel tests along different East Asian monsoon directions were performed in IBDWS (available at http://ibdws.sdsu. edu/~ibdws/) based on the pairwise F ST between populations using 1,000 randomizations 59 .
Phylogenetic analyses and estimates of divergence time. The phylogenetic relationships among chlorotypes and ribotypes were reconstructed by PhyML version 3.0 with maximum likelihood (ML) 60 and by MrBayes Version 3.2.1 with mixed models using Bayesian inference (BI) and MCMC sampling 61,62 , respectively. The best-fit nucleotide substitution model for the codon position was evaluated by the Akaike Information Criterion 63 with the program jModeltest Version 0.1.1 under the ML criterion 64 and then by GTR+ G and TrN+ I+ G models for cpDNA and nrITS, respectively. The ML analyses were conducted with the heuristic search option over 100 random-taxon-addition replicates with tree bisection-reconnection (TBR) branch swapping. Trees were sampled for every 100 generations, resulting in 100,000 trees, and the first 25% were discarded as burn-in. Bayesian analyses were performed over 50 million generations with one cold and three incrementally heated Markov Chain Monte Carlo (MCMC) analyses, starting from a random tree and sampling every 1,000 generations. The first 10% of trees were treated as burn-in and discarded, and all remaining trees were used to construct a Bayesian consensus tree. All trees were then modified using FigTree version 1.3.1 65 .
To estimate the divergence time between the genealogical lineages, we carefully chose a mutation rate for cpDNA of 4.87 × 10 −10 s/s/y, which is suitable for shrub species 27 . Since neither any pollen nor fossil records of Reaumuria could be found from the deposit strata published, we chose two closely related genera (Tamarix and Myricaria) within Tamaraceae to calibrate the molecular clock. Based on this mutation rate and its assumed standard deviation (10%), the divergence time and coalescence time of lineages were estimated using a Bayesian approach implemented in BEAST Version 1.7.5 with the criteria set by BEAUti (available at http://BEAST.bio.ed.ac.uk) 65 . All simulations were based on an uncorrected log-normal relaxed clock distribution of branch lengths, started with a prior constant population size and a randomly generated starting tree. The posterior distributions of the parameters in the MCMC analyses were approximated with 80 million steps in each analysis and sampled every 8,000 generations. Convergence of the parameters sampled was checked with the program TRACER version 1.5 to examine the highest effective sampling size values (ESSs > 200) for all parameters 66 , and the burn-in steps were discarded to obtain an estimation of the posterior probability distribution of divergence time at the relevant node. The maximum lineage credibility tree was then summarized by TreeAnnotator version 1.7.5 (available at http://BEAST.bio.ed.ac.uk/TreeAnnotator/).

Inference of demographic history.
Possible population expansions were assessed by neutrality tests with Tajima's D and Fu's Fs and mismatch distribution analysis (Supplementary Method S2). MIGRATE-N version 3.6.4 25 was used to estimate the variation in the historical effective population size. The default settings were conducted in the MIGRATE-N program with some modified options as follows: the Bayesian inference module; one single long run using heating with temperatures of 1.0, 1.5, 3.0 and 10,000, totaling 50,100,000 visited parameter and genealogy changes in the cold (1.0) chain; sampling 10 replicates of each 50,000 in intervals of 100; after discarding the first 100,000 visits. The mutation-scaled population size theta (θ ) was displayed by the uniform prior distributions with ranges from 0 to 0.1, namely the effective population size Ne times the mutation rate μ per site and generation. Moreover, Bayesian skyline plots (BSPs) for population dynamics were also conducted on cpDNA (Supplementary Method S3).

Ecological niche modeling.
To predict the distribution of R. soongarica based on habitat suitability, ecological niche modeling was conducted with Maxent 26 with the WorldClim database from 1955 to 2000 (http://www.worldclim.org/). The pairwise correlation for all 19 biological factors were tested, and 12 independent biological factors were used for further predictions. The distribution sites consisted of those examined in the field investigation together with the specimen records from the Chinese Virtual Herbarium and with coordination from the Global Positioning System (http://www.cvh.org.cn/); 75% of the sites were randomly selected to design a prediction scheme on consideration of the commission error, and the remaining 25% of the sites were used for evaluations of the goodness of the prediction.
The contemporary species-climate relationships were then projected to past climate layers to predict the potential range of the species during LGM (ca. 21,000 years BP) and LIG (ca. 120,000 -140,000 years BP) periods. The project bioclimatic layers for LIG (resolution of 30 arc-second) and LGM (resolution of 2.5 arc-minute) were obtained from the WorldClim database version 1.4 67 . We employed a "maximum training sensitivity plus specificity" threshold to determine suitable/unsuitable habitat because sensitivity-specificity combined approaches have been suggested as good thresholds 68 . An AUC value estimation greater than 0.9 for each estimate was used as the criterion for a good model. Additionally, 100 permutations were conducted in the modeling process for LGM and LIG periods based on the MIROC model. The relative contribution of each variable was evaluated with the covariance in terms of increased predictive performance resulting from a given environmental variables. Only the ecological variables that explained greater than 10% of the prediction variance were viewed as the principal habitat suitability factors limiting the distribution range of R. soongarica.