Molecular Phylogeography Analysis Reveals Population Dynamics and Genetic Divergence of a Widespread Tree Pterocarya stenoptera in China

The geological events, past climatic fluctuations, and river systems played key roles in the spatial distribution, population dynamics, and genetic differentiation of species. In this work, we selected Pterocarya stenoptera, a widespread tree species in China, to test the roles of these factors. Four noncoding spacers, eight microsatellite (simple sequence repeat) markers, and species distribution modeling were used to examine the phylogeographical pattern of P. stenoptera. Based on chloroplast DNA data, populations of P. stenoptera were clearly clustered into three groups. The divergence time of these groups fell into the stage of the Qinghai–Tibet Movement, 1.7–2.6 Ma. For simple sequence repeat data, only one western marginal population YNYB could be separated from other populations, whereas other populations were mixed together. Our results indicated that the environmental heterogeneity resulting from the Qinghai–Tibet movement might be response for this genetic divergence. The climatic fluctuations in the Pleistocene did not cause the substantial range shift of P. stenoptera, while the fluctuations affected its population size. Moreover, we also confirmed the river systems did not act as channels or barrier of dispersal for P. stenoptera.


INTRODUCTION
The uplift of the Qinghai-Tibet Plateau (QTP) has not only changed the topography of the area and its surroundings but also played a remarkable role in the climate change of Asia Deng and Ding, 2015). The complex terrain and climate caused by this uplift exacerbated the rapid differentiation of species on the plateau, thus transforming the QTP and its adjacent areas into one of the world's richest regions in China (Wu, 1987). Furthermore, the uplift of the QTP caused a profound impact on species differentiation in other regions of Asia. The fast uplift of the QTP, i.e., the Qinghai-Tibet movement (since ca. 3.6 Ma), caused a remarkable increase in the eolian flux in the north Pacific, and this increase in eolian flux has remained steady until today (Li et al., 1979;Rea et al., 1998). As a consequence of this change, the Asian inland began to experience aridification, while southeast China became increasingly humid . To date, little attention has been paid to whether climate change in Asia caused by the Qinghai-Tibet movement has affected species differentiation (Liu et al., 2013). Most of the research attention has focused on the impact of the recent uplift of the QTP, i.e., Kunlun-Yellow River movement (Liu et al., 2012). This movement, which occurred approximately 1.1-0.7 Ma, affected the whole QTP and lifted the plateau surface to 3,000-3,500 m (Cui et al., 1998;Shi, 1998). Even some areas of the mountain ranges were uplifted to a height of 4,500-5,000 m, thus allowing the development of mountain glaciers in response to global glacial-interglacial cycles (Shi, 1998;Zhou et al., 2006). The movement also influenced the pattern of atmospheric circulation and formed the modern pattern of the East Asian monsoon (Jiang et al., 2005). Most of East Asia was not covered by glaciers, but the mountain glacial cycles after the Kunlun-Yellow River movement notably caused the climatic fluctuation in this area (Shi et al., 2006). The climatic oscillations during the Pleistocene substantially shaped the spatial distribution, population dynamics, and genetic differentiation of species (Comes and Kadereit, 1998;Hewitt, 2000;Hewitt, 2004).
Molecular phylogeographical methods seek to unravel historic biogeographical processes in response to climatic oscillations during the Pleistocene (Petit et al., 2003;Avise, 2009;Hickerson et al., 2010). The rise in phylogeographical studies for an increasing number of species and the resultant accumulated evidence suggest that the species have responded to climate fluctuations in two ways (Fu et al., 2014;Miao et al., 2017a). The first type of response is the adjustment in distribution range of these species when they migrated southward or to main refugia during glaciation and their recolonization after glaciation; the other type is the persistence in situ in the form of local adaptation during multiple glacial-interglacial cycles Miao et al., 2017b). Response modes are determined by the tolerance and migration capacity of species. However, regardless of the response patterns, species population dynamics are more likely affected by the Pleistocene climatic fluctuations. The glacialinterglacial cycle (0.7 Ma to the present), which was strengthened by the Kunlun-Yellow River movement, occurred many times in East Asia (Cui et al., 2011). Which cycle has the greatest impact on species population dynamics in this region during the Pleistocene?
Pterocarya stenoptera C. DC (Juglandaceae) is a dominant deciduous broad-leaved tree growing in forests below 1,500 m above sea level. The species is widely distributed in the warmtemperate and subtropical regions in China. P. stenoptera is intolerant of cold and drought, and it is likely to grow in warm and humid environments (Xu et al., 2002). Flowers bloom in April and May, and fruit ripening extends to August and September. The samara of P. stenoptera becomes dry after ripening and can be dispersed by wind and water flow. According to our field observations, P. stenoptera is mainly distributed along stream banks. Thus, knowing whether rivers have an impact on genetic patterns of P. stenoptera is worth discussing. The three relatively large rivers in the distribution area of P. stenoptera are the Yellow River, the Huaihe River, and the Yangtze River. In previous phylogeographical studies, rivers serve as physical barriers or dispersal channels Wang et al., 2014;Geng et al., 2015;Cazé et al., 2016;Liang et al., 2018). P. stenoptera, as a dominant deciduous broad-leaved tree of China, is therefore an ideal candidate in the investigation of the influences of geological events, climate change, and river systems on the genetic differentiation and population dynamics of species in this region.
Here, we employed chloroplast DNA (cpDNA) noncoding spacers, microsatellite [simple sequence repeat (SSR)] markers, and species distribution modeling (SDM) to investigate the phylogeographical pattern of P. stenoptera. Our specific aims were to address two main questions: (i) What is the genetic structure of P. stenoptera populations in China as revealed by cpDNA data and SSR markers? (ii) How do geological events, climate change, and river systems affect the genetic differentiation and population dynamics of P. stenoptera?

Population sampling
Fresh leaves of 440 samples of P. stenoptera distributed in 22 localities were collected covering its entire distribution area in China during May to November 2017 (Figure 1 and Table 1). Each population collects 20 individuals, and they are at least 20 m apart. The leaves were dried immediately in silica gel and stored at room temperature. The geographical information of these populations and numbers of individuals used in the cpDNA and SSR analyses were presented in Table 1. No specific permits are required for the collection of P. stenoptera, and all samples were collected under current national and local regulations.
DNA Isolation, ssR Genotyping, and DNA sequencing Total genomic DNA was extracted from 30 mg of silicadried leaves using a DNA Extraction Kit DP305 (Tiangen, Beijing, China). Four noncoding cpDNA spacers of ndhAx1-ndhAx2, rpL34-rpL36, trnL(UAG)-rpL32-F (Shaw et al., 2007), and trnV(UAC)x2-ndhC (5′-TGCATGTTGGGTCTTTGA, 5′-CATCGCCCATTGGTTCTA) and eight microsatellite markers (Ps23, Ps28, Ps59, Ps69, Ps70, Ps75, Ps82, and Ps83) (Wang et al., 2018) were amplified for 438 and 440 samples respectively (Table 1). For cpDNA fragments, all the polymerase chain reactions (PCRs) were performed in a total volume of 40 µl, containing 30 ng of template DNA, 10 mM Taq buffer, 0.2 mM dNTPs, 0.3 mM primers, and 2 units of Taq polymerase (Tiangen, Beijing, China). PCRs were performed in a Mastercycler nexus thermocycler (Eppendorf, Hamburg, Germany). The PCR program was set to: an initial denaturation at 95°C for 5 min, 35 cycles of denaturation at 95°C for 40 s, 40 s of annealing at 52-60°C and 1 min at 72°C, and a final extension of 72°C for 5 min. The PCR products were sequenced using an ABI 3730 DNA Sequencer. The quality of original chromatogram of sequence was checked and aligned using CLUSTAL_X version 2.1 (Larkin et al., 2007). For SSR markers, PCRs were performed in a 20-µl reaction mixture consisting of 20 ng genomic DNA, 10 mM Taq buffer, 0.2 mM dNTPs, 0.3 mM primers, and 1 unit of Taq polymerase. The forward primers of SSRs were 5′ labeled with FAM, HEX, or TAMARA according to Wang et al. (2018). The PCR program was set to: an initial denaturation at 95°C for 5 min, 35 cycles of denaturation at 95°C for 40 s, 40 s of annealing at 48-62°C and 40 s at 72°C, and a final extension of 72°C for 5 min. PCR products were mixed with 10 μl of HiDi formamide and 0.1 μl of internal standard ROX500 (Applied Biosystems, Foster City, USA). These products were genotyped on an ABI 3730 DNA Analyzer.

DATA ANAlYsIs cpDNA Data
The cpDNA sequences of each noncoding spacer were aligned by CLUSTAL_X version 2.1 (Larkin et al., 2007). The indels in the cpDNA sequences were encoded as described by Caicedo and Schaal (2004). The genetic parameters of haplotype diversity (h) and nucleotide diversity (π) were estimated using dnasp version 6.12 (Rozas et al., 2017). The number of private haplotypes (N Ph ) was counted by manual statistics. The networks of cpDNA haplotypes were constructed using TCS version 1.21 (Clement et al., 2000) with >20-step limitation. Divergence time, which was estimated by BEAST version 2.5.0 (Bouckaert et al., 2014), was used to relate the genetic differentiation of cpDNA lineages to the historic events of this species. The statistical parsimony network of 28 haplotypes were constructed. According to Zhang et al. (2013), Juglans mandshurica is a species close to P. stenoptera, while Platycarya strobilacea is a relatively distant species. Thus, we selected these two species as the outgroups before running BEAST module, and the best fit nucleotide substitution model of Hasegawa, Kishino, and Yano (HKY) was selected in MEGA version 10.0 (Kumar et al., 2018) in accordance with the Bayesian information criterion. Dating analysis was performed with a lognormal relaxed molecular clock model, four gamma categories, the Yule process of tree prior, HKY as site model, and a normal prior distribution for age constraints. Normal distribution prior root sets were used for the node ages to unite J. mandshurica = 39.72 Ma and P. strobilacea = 63.05 Ma in secondary calibration based on the results of Zhang et al. (2013). The sequences of J. mandshurica were deposited in GenBank with accession numbers MK002388, MK002394, MK002404, and MK002416. The sequences of P. strobilacea were downloaded from GenBank with accession number KX868670. We set the Markov chain Monte Carlo (MCMC) chain to 10 million generations at the sampling frequency of 1,000 generations each time. Tracer version 1.7.1 was used to visualize and check for convergence to a stationary distribution and ensure effective sampling size values (ESSs > 200). We discarded the first 10% as burn-in and constructed the maximum clade credibility tree by using TreeAnnotator version 2.5.1 (Bouckaert et al., 2014). The resulting tree, with ages appropriated to each node and 95% credibility intervals set as the divergence times, was drawn in FigTree version 1.4.3 (Rambaut, 2016). The genetic relationships of 22 populations of P. stenoptera were constructed the basis of their pairwise genetic distances in cpDNA and determined from the net average between populations of sequences (D A ) by using the unweighted pair-group method with arithmetic module of MEGA version 10.0 (Kumar et al., 2018). The D A value was estimated using MEGA version 10.0 (Kumar et al., 2018) under the maximum composite likelihood model (Tamura et al., 2004).
We retraced the demographic history of P. stenoptera by using the extended Bayesian skyline plot (EBSP) in BEAST version 2.5.0 to further infer potential relatively complex changes given an effective population size (Bouckaert et al., 2014). EBSP analyses were performed with a lognormal relaxed molecular clock model, HKY as the site model, and a normal prior distribution for age constraints. A normal distribution prior root set was used for the node-age results of the above calculations. The scale factor of 0.5 for cpDNA was adopted because only the female cpDNA contributes to effective population size, and then the MCMC chain was run in 10 million generations at a frequency of 1,000 generations each time.
Inferences for the historical demographic processes were generated by calculating Tajima's D values (Tajima, 1989) and Fu's Fs values (Fu, 1997) for the total level and for each group in arlequin version 3.5.1 (Excoffier and Lischer, 2010), respectively. f a demographic expansion hypothesis is true, then significant  Table 1). The dotted line represents the dividing line of three genetic groups that identified by cpDNA data. (B) The statistical parsimony network of haplotypes H1-H28. The size of circles corresponds to the frequency of each haplotype. Each solid line represents one mutational step interconnecting two haplotypes for which parsimony is supported at the 95% level, the black dotted line represents multiple mutational steps that connected with outgroups. The colorful dotted line represents the dividing line of three genetic groups that identified by cpDNA data.  Nei, 1987), allele frequencies, and allelic richness (A R ) were calculated by FSTAT version 2.9.3 (Goudet, 1995). The number of private alleles (N PA ) was manually counted according to the allele frequencies in FSTAT. The role of rivers in the transmission of P. stenoptera was explored by performing correlation analyses between the level of genetic diversity based on SSR and cpDNA and longitude using the Mantel test in vegan package (Oksanen et al., 2017). Similarly, correlation analyses between the level of genetic diversity based on SSR and cpDNA and latitude were also performed using the Mantel test in vegan package. Bayesian analysis was performed in structure version 2.3.4 (Pritchard et al., 2000) to infer the most likely number of genetic subgroups (K) in the SSR dataset. The program was run with K values from 1 to 10 in 10 replicates for each K, 20,000 burn-in periods and 100,000 MCMC iterations. An admixture model with correlated allele frequencies was chosen for this analysis. The optimal K value was determined by computing the ∆K values according to the method of Evanno et al. (2005); this method was executed in Structure Harvester (Earl and vonHoldt, 2012). The best K value was further estimated according to the value of the TI posterior probability by using RmavericK (Verity and Nichols, 2016). The average value of the admixture coefficients over 10 runs was calculated using CLUMPP version 1.1 (Jakobsson and Rosenberg, 2007). The bar plots of STRUCTURE were produced by DISTRUCT 1.1 (Rosenberg, 2004). A principal component analysis of the SSR data was performed using ggord package (Beck, 2017). For both SSR and cpDNA data, tests of isolation by distance (IBD) were performed by regressing the standardized values of F ST against geographical distance (m) and the ecological distance via the Mantel permutation procedure in vegan package (Oksanen et al., 2017). The values of F ST based on the SSR and cpDNA data were estimated by FSTAT and dnasp. Geographical distance was calculated by the geosphere package (Hijmans, 2015), and environmental distances were calculated as Euclidian distances along the principal components (PCs) according to the method of Pluess et al. (2016). The first four PCs can explain 91.11% of the variation for the present extracted 19 environmental variables. Nineteen environmental variables at 2.5 arcmin resolution for the present time were downloaded from the WorldClim website. Genetic differentiations of cpDNA data and SSR data among populations and for populations within each group (F ST ) were calculated by nonhierarchical analyses of molecular variance (AMOVAs) in ARLEQUIN version 3.5.1 (Excoffier and Lischer, 2010). The amount of genetic variations among groups, among populations within groups, and within populations was estimated with hierarchical AMOVAs. The significance of AMOVA was tested with 1,000 permutations. Gene flow (Nm) among populations for the cpDNA and SSR data was estimated as 1/4(1/F ST − 1) (Wright, 1951). Pollen/seed migration ratio (r) was estimated using a modified equation of Ennos (1994), according to Petit et al. (2005)

species Distribution Models
The potential range shift of P. stenoptera was determined by SDM in Maxent version 3.4.1 (Phillips et al., 2006). The geographic coordinates used in SDM analysis were based on a set of 40 presence points (Table S1), in which 18 points were adopted from the Chinese Virtual Herbarium, while the other 22 points represent our sampling sites. Nineteen environmental variables at 2.5 arcmin resolution at the present and the last glacial maximum (LGM) were downloaded from the WorldClim website. The calculation method used to determine potential climatically suitable areas was based on the study of Fu et al. (2014). However, logistic probabilities were used for the output. Logistic probabilities higher than 0.5 were taken as high climatically suitable areas in accordance to the classification of Bai et al. (2016). The significant difference of environmental variables among groups was determined by T-test in SPSS version 19 (Inc., Chicago, IL, USA). According to Li et al. (2015), the result of the Qinghai-Tibet movement can be used to represent the difference of warmth and drought between the northern and southern regions of China. Thus, five environmental variables (Table S2) related to warmth and drought, including annual mean temperature (Bio1), mean temperature of coldest quarter (Bio11), annual precipitation (Bio12), precipitation of coldest quarter (Bio19), and water vapor pressure in January (wvp1), were selected. The corresponding values, which were downloaded from the WorldClim website, covered the period of 1970-2000 at 2.5 arcmin resolution, then they were further extracted using DIVA-GIS 7.5 (Hijmans et al., 2001).
The populations of HBSN and FJWY had the highest nucleotide and haplotype diversities, while the populations of HBSN and HNHM had highest number of private haplotypes ( Table 1). As for SSR data, 440 individuals of P. stenoptera were genotyped. The results of genetic diversity analyses of each population are presented in Table 1. The level of genetic diversity (H E ) of each population ranged from 0.194 (JXLH) to 0.397 (HBJG). The level of allelic richness (A R ) of each population ranged from 1.931 (SDTM) to 2.842 (YNYB). Among these populations, six populations were observed with private alleles, and the population of SDMM showed the highest number of private alleles. Mantel test showed that the genetic diversity of SSR (r = 0.004, P = 0.390; Figure 2) and cpDNA (r = 0.016, P = 0.287; Figure 2) did not decrease significantly along the longitude direction from west to east. No significant decrease was observed in genetic diversity along the latitude direction (SSR: r = 0.063, P = 0.240; cpDNA: r = −0.020, P = 0.545; Figures 2C, D) from south to north.

Genetic structure
Populations of P. stenoptera were significantly structured as revealed by the population-based unweighted pair-group method with arithmetic tree based on cpDNA data, and the 22 populations were subdivided into three genetic groups ( Figure  S1). The three identified genetic groups were the North group, the South group, and the Southwest group. The network of haplotypes also supported this subdivision (Figures 1 and 3). These groupings were largely consistent with their geographic distributions. Three haplotypes were shared among groups. In particular, H7 was shared between FJWY from the South group and SDTM, SXWZ, JSBH, JSLM, HNJG, and HBSN from the North group; H17 was shared between ZJTM and HNTM from the South group and the North group; and H23 was shared between SCWD from the Southwest group and SCEM from the South group. According to the Bayesian analysis of population structure, the highest likelihood based on SSR data was obtained when K = 6 ( Figure S2). The best K = 1 was estimated according to the value of TI posterior probability ( Figure S3). The cluster results on TI posterior probability were also supported by principal component analysis. No significant genetic differentiation was found among 22 populations (Figure 4). However, the actual results when K = 6 in the structure analysis showed that only the population of YNYB can be separated from other populations, whereas the other 21 populations were mixed together and cannot be separated from one another (Figure 5). Taken together, the cluster results based on SSR data did not correspond to the separate geographical regions supported by cpDNA data (Figure S1). As for cpDNA data, a strong population genetic differentiation (F ST = 0.825; P < 0.001) was detected by the analyses of molecular variance (AMOVA) ( Table 2). The relatively low gene flow (Nm = 0.053) among populations resulted in significant genetic differentiation. No significant isolations by geographical distance (r = −0.019, P = 0.695) and ecological distance (PC1: r = 0.036, P = 0.128; PC2: r = −0.055, P = 0.946; PC3: r = −0.067, P = 0.988; PC4: r = −0.026, P = 0.818) for cpDNA were detected at the species-range scale. Nonsignificant isolation by ecological distance for SSR was also detected along PC1, PC3, and PC4 (r = 0.033, P = 0.150; r = −0.029, P = 0.699; r = −0.027, P = 0.791, respectively), but a significant IBE existed along PC2 (r = 0.703, P = 0.001). However, significant isolation by geographical distance (r = 0.488, P = 0.001) for SSR was detected. AMOVA showed that 54.93% of the total cpDNA genetic differentiation existed among the three groups (F CT = 0.549), and 31.50% among populations within groups (F SC = 0.699), and 13.58% within populations (F ST = 0.864) ( Table 2). As for SSR data, 2.89% of the total genetic differentiation was found among the three groups (F CT = 0.029), 6.34% among populations within groups (F SC = 0.065), and 90.78% within populations (F ST = 0.092) ( Table 2). The high pollen/seed migration ratio (r) = 52.0 suggests that the high gene flow was mostly transmitted by pollen.

Divergence Time
Divergence time in the time-calibrated tree (Figure 3) ranged from 3.16 Ma (95% HPD = 1.79-5.02) to 0.37 Ma (95% HPD = 0.01-0.88). The divergence time of the North group and the South-Southwest groups was 3.16 Ma (95% HPD = 1.79-5.02), which falls within the A to B stages of the Qinghai-Tibet movement (2.6-3.6 Ma). The divergence time of the South group and the Southwest group was 2.54 Ma (95% HPD = 1.12-3.04), which falls within the B to C stages of the Qinghai-Tibet movement (1.7-2.6 Ma).

Demographic history
The neutrality test with Tajima's D and Fu's Fs statistics showed negative values at the total level (Tajima's D = −0.438, P = 0.395; Fu's Fs = −3.790, P = 0.201; Table 3), suggesting that the demographic expansion has occurred in the recent past; however, the obtained  Table 1.
FIGURe 5 | Estimated genetic structure for K = 6 obtained with the program STRUCTURE for 22 populations of Pterocarya stenoptera based on simple sequence repeat (SSR) data. Each vertical bar represents an individual and its assignment proportion into one of five (colored) population clusters (K).
values were statistically nonsignificant. This demographic expansion was obviously supported by the nonsignificant sum of squared deviations (SSD) and raggedness (RAG) index values of mismatch distribution analysis at the total level (Table 3). Tajima's D and Fu's Fs statistics and mismatch distribution analysis were examined for each of the three groups. The demographic expansion of the South group was supported consistently by the above three methods, whereas the demographic results of the other two groups were inconsistent. EBSP analysis was performed to infer the time of demographic expansion in this study. The results of EBSP indicate a continuous demographic expansion that began approximately 0.40 Ma (Figure 6), and this result is highly consistent with the time of Marine Isotope Stage (MIS) 11 (0.33-0.46 Ma).

DIsTRIBUTION INFeReNCe BY sPeCIes DIsTRIBUTION MODels
The values of the areas under the receiver operating characteristic curve based on both training and test presence data at the present and lgm periods were 0.994 And 0.990 (Present) 0.994 And 0.992 (Lgm), respectively. The results of areas under the receiver operating characteristic curve values demonstrated good model performance. The results of the predicted distribution of p. Stenoptera indicate that the climatically suitable areas (logistic probabilities > 0.25) At present and lgm were both continuous in china, and no habitat fragmentations were isolated by intervening unsuitable habitats (Figure 7). The area of climatically suitability (logistic probabilities > 0.25) at the lgm (2,081,410 km 2 ) was not compressed compared with that at the present time (1,950,755 km 2 ). Similarly, no obvious difference existed between the area of high climatically suitability (logistic probabilities > 0.50) At the lgm (1,221,637 km 2 ) and that at the present time (1,172,061 km 2 ). The significant difference of the environmental variables (bio1, bio11, bio12, bio19, and wvp1; Table S2) among the groups was further determined by t-test in spss. The results showed that the environmental variables, bio1, bio11, bio12, and wvp1 between the north group and the South-southwest groups were significantly different (p < 0.05). Moreover, bio12 and bio19 between the South group and the Southwest group were significantly different (p < 0.05).

DIsCUssION
spatial Population Genetic structure of P. stenoptera The spatial genetic structure of species is widely known to be linked to multiple factors (Ohsawa and Ide, 2008;Hickerson et al., 2010;Yang et al., 2017). In this study, we investigated the spatial genetic structure of P. stenoptera across its distribution area in China to provide insights into the complex interplay among multiple factors. Our survey showed that the 22 populations of P. stenoptera were divided into three groups, i.e., North group, South group, and Southwest group.  Previous phylogeographical studies reported large-scale intraspecific disjunctions in many species that can alternatively be explained by geographical barrier, ecological barrier, and recolonization from multiple refugia (Li et al., 2011;Li et al., 2012;Bai et al., 2016;Luo et al., 2017;Ye et al., 2017). If the subdivision of P. stenoptera was caused by geographical barriers among groups, then the premise of this hypothesis was that geographical barriers existed among these groups. For the distribution area of P. stenoptera, no large mountains isolated the groups, but three major rivers, i.e., the Yellow River, Huaihe River, and Yangtze River, may have divided them. According to the literature, the Yellow River and Yangtze River, as phylogeographical boundaries, have promoted the genetic divergence of plant species that distributed across these rivers (Wang et al., 2014;Geng et al., 2015). However, the Yangtze River and Yellow River formed at 1.9-1.7 Ma (Wang et al., 1998;Li et al., 2001), while the Huaihe River formed at the end of the Late Pleistocene (ca. 0.15 Ma; Shao et al., 1989), which are clearly later than the divergence time of the North group and the South-Southwest groups (3.16 Ma), and the divergence time of the South group and the Southwest group (2.54 Ma). The findings indicate that these river system formations were not the vicariant factor that have caused the divergence of the three groups. Furthermore, the topologies of these rivers were inconsistent with the geographical boundaries of the three groups (Figure 1), indicating that these rivers also did not play a role in enhanced isolation after genetic differentiation. Thus, the divergence of the three groups is unlikely caused by the geographical isolation of P. stenoptera. Moreover, the possibility that these larger rivers served as dispersal channels for P. stenoptera still needs to be explored. The flow direction of the three large rivers are from west to east, but our study did not find a gradually decline of genetic diversity in this direction ( Figure  2). Therefore, these rivers did not act as channels of dispersal of P. stenoptera.
We used SDM to simulate the distribution of P. stenoptera at the present and at LGM to further test whether this divergence was caused by ecological isolation. The results did not show habitat fragmentations isolated by intervening unsuitable habitats at the present and at LGM. Our results of nonsignificant IBE based on cpDNA and SSR data (except PC2) also do not support the assumption that divergence is caused by ecological isolation. Moreover, the Pleistocene glacial epoch occurred after the Kunlun-Yellow River movement (since ca. 1.1 Ma; Li et al., 2015), suggesting a much later than the time of species group differentiation. These findings indicate that the much even earlier glaciations cannot be the reason for the differentiation of these groups. In the assumption that the genetic divergence of P. stenoptera is caused by recolonization from multiple refugia, the movement of populations away from a refugia will likely cause a gradual decline in genetic diversity. However, neither the North group nor the South group underwent a gradient descent from the presumed refugia, i.e., these populations have the highest genetic diversity. Here, the Southwest group was not surveyed considering the relatively small number of populations in this area. We also did not find a gradual decline in genetic diversity from south to north. The results of nonsignificant IBD based on cpDNA data also do not support this particular hypothesis. The IBD for SSR data was significant, but they were likely influenced by the strong gene flow via pollination. Thus, the results cannot clearly reflect a migration route that occurred in the past. Obviously, our results do not support this interpretation of recolonization from multiple refugia. Taken together, the reasons mentioned above do not seem to explain the genetic differentiation among the three groups.
We then focused on the divergence time between groups. The divergence time of the North group and the South-Southwest group was 3.16 Ma. This divergence time falls within the time of the Qinghai-Tibet movement (1.7-3.6 Ma; Li and Fang, 1998). Consequently, the Qinghai-Tibet movement caused the southeast China, i.e. the South and Southwest groups for P. stenoptera, to become more humid until today .
Our T-test results also showed that the environmental variables related to temperature, precipitation, and air humidity (Bio1, Bio11, Bio12, and wvp1; Table S2) between the North group and the South-Southwest groups were significantly different. The divergence time of the South group and the Southwest group was 2.54 Ma, which falls within the time of the Qinghai-Tibet movement. The Southwest group was closer to Qinghai-Tibet than the South group, and the environmental variables related to precipitation (Bio12 and Bio19) between the South group and the Southwest group were significantly different. The environmental heterogeneity caused by the Qinghai-Tibet movement may be related to the genetic differentiation of the three groups of P. stenoptera. Nonetheless, we cannot completely rule out genetic drift and ancient hybridization. Random genetic drift or ancient hybridization to capture the chloroplast genes from other species may also give rise to the result of our study.
For the structure analysis of the SSR data, the results showed that only one population of YNYB can be separated from the other populations, while the other 21 populations were mixed together. This phenomenon is due to the high gene flow among populations (Nm = 2.836), which is mainly derived from pollen transmission with the high pollen/seed migration ratio (r = 52.0). Efficient gene flow (Nm > 1) will break the isolation between populations (Wright, 1931). Although the pollen-mediated gene flow obscured previous genetic structures, the YNYB population at the edge of the distribution area is not mixed with other populations. Efficient gene flow can also explain the access to private haplotypes of only a handful population. However, the eastern marginal population SDMM has an unusual number of private genes, which may be the result of genetic drift.
Population Dynamics for P. stenoptera in Response to Climatic Oscillations During the Pleistocene Previous phylogeographical studies suggested that the species distributed in China during the Pleistocene glacial period migrated southward or retreated to main refugia as a means to cope with the climate fluctuation Liu et al., 2012). Recent phylogeographical studies offer another alternative hypothesis of microrefugia, i.e., species persisted in situ by local adaptation during the multiple glacial-interglacial cycles during the Pleistocene (Fu et al., 2014;Wang et al., 2014;Zhang et al., 2015;. In the case of P. stenoptera, the genetic differentiation of the three groups with obvious boundaries and the simulation of species distribution for the LGM and present periods obviously did not support the assumption of range shift by southward migration. The results of SDM indicate that the area of climatically suitability at LGM was not compressed compared with that at the present time. Moreover, a gradient decline in genetic diversity was lacking due to the post-glacial colonization of the groups of P. stenoptera. Thus, the range in the form of retreat to the main refugia for P. stenoptera during the Pleistocene glacial period was not supported. On the basis of the cpDNA data, most of the populations (13 of 22 populations) of P. stenoptera have their own private genes. The results indicate the microrefugia hypothesis are likely suitable for P. stenoptera. Although our results did not support the range shift of P. stenoptera during the climatic fluctuations, its population size was affected. Tajima's D and Fu's Fs statistics and the mismatch distribution analysis all supported P. stenoptera experienced a population expansion. The results of EBSP indicate that this expansion occurred at approximately 0.40 Ma whose timeline is highly, consistent with the interglacial of MIS 11 (0.33-0.46 Ma; Cui et al., 2011). In other words, P. stenoptera populations experienced continuous population expansion after MIS 12, i.e., the Zhonglianggan glaciation. Overall, the climatic fluctuation in the Pleistocene did not cause the substantial range shift of P. stenoptera, but its population size was affected. Therefore, our results support the assumption that P. stenoptera has not been compressed into several main refugia during the Pleistocene. The previously proposed microrefugia hypothesis (Miao et al., 2017a) seems to be more suitable for P. stenoptera, which underwent local adaptation in response to the climate fluctuations during the Pleistocene.
The main goal of this study was to understand how geological events, climate change, and river systems affect the genetic differentiation and population dynamics of P. stenoptera. The phylogeographic study of P. stenoptera showed that the species is genetically highly diverse based on the cpDNA data. Our results support the assumption that the environmental heterogeneity, which was caused by the climate change resulting from the Qinghai-Tibet movement may be linked to the genetic divergence. The climatic fluctuations during the Pleistocene did not cause the substantial range shift of P. stenoptera, but the fluctuations affected its population size. We also confirm in this study that the river systems in the area did not act as channels or barriers of dispersal for P. stenoptera.

DATA AVAIlABIlITY sTATeMeNT
Publicly available datasets were analyzed in this study. This data can be found here: accession numbers MK002380 to MK002387, MK002389 to MK002393, MK002395 to MK002403, and MK002405 to MK002415.

AUThOR CONTRIBUTIONs
YL conceived the research project and wrote the paper. Z-HQ, M-WL, and J-XL collected the data. Z-HQ, Y-XH, M-WL, and X-FY analyzed the data. All authors are in agreement with the content of the manuscript.

ACKNOWleDGMeNTs
We are grateful to Jie Yang and Cai-Yun Miao for help during the field work for collecting samples.