Late Pleistocene speciation of three closely related tree peonies endemic to the Qinling–Daba Mountains, a major glacial refugium in Central China

Abstract Determining the factors promoting speciation is a major task in ecological and evolutionary research and can be aided by phylogeographic analysis. The Qinling–Daba Mountains (QDM) located in central China form an important geographic barrier between southern subtropical and northern temperate regions, and exhibit complex topography, climatic, and ecological diversity. Surprisingly, few phylogeographic analyses and studies of plant speciation in this region have been conducted. To address this issue, we investigated the genetic divergence and evolutionary histories of three closely related tree peony species (Paeonia qiui, P. jishanensis, and P. rockii) endemic to the QDM. Forty populations of the three tree peony species were genotyped using 22 nuclear simple sequence repeat markers (nSSRs) and three chloroplast DNA sequences to assess genetic structure and phylogenetic relationships, supplemented by morphological characterization and ecological niche modeling (ENM). Morphological and molecular genetic analyses showed the three species to be clearly differentiated from each other. In addition, coalescent analyses using DIYABC conducted on nSSR variation indicated that the species diverged from each other in the late Pleistocene, while ecological niche modeling (ENM) suggested they occupied a larger area during the Last Glacial Maximum (LGM) than at present. The combined genetic evidence from nuclear and chloroplast DNA and the results of ENM indicate that each species persisted through the late Pleistocene in multiple refugia in the Qinling, Daba, and Taihang Mountains with divergence favored by restricted gene flow caused by geographic isolation, ecological divergence, and limited pollen and seed dispersal. Our study contributes to a growing understanding of the origin and population structure of tree peonies and provides insights into the high level of plant endemism present in the Qinling–Daba Mountains of Central China.

quences to assess genetic structure and phylogenetic relationships, supplemented by morphological characterization and ecological niche modeling (ENM). Morphological and molecular genetic analyses showed the three species to be clearly differentiated from each other. In addition, coalescent analyses using DIYABC conducted on nSSR variation indicated that the species diverged from each other in the late Pleistocene, while ecological niche modeling (ENM) suggested they occupied a larger area during the Last Glacial Maximum (LGM) than at present. The combined genetic evidence from nuclear and chloroplast DNA and the results of ENM indicate that each species persisted through the late Pleistocene in multiple refugia in the Qinling, Daba, and Taihang Mountains with divergence favored by restricted gene flow caused by geographic isolation, ecological divergence, and limited pollen and seed dispersal. Our

| INTRODUC TI ON
It is widely known that climatic fluctuations during the Pleistocene strongly impacted the distribution and genetic structure of plant and animal species in both Northern and Southern Hemispheres (Hewitt, 2004;Lu, Heckel, Liang, & Zhang, 2013;Song et al., 2009). Although China experienced less severe glaciations during this period relative to Europe and North America (López-Pujol, Zhang, Sun, Ying, & Ge, 2011;Qian & Ricklefs, 1999Rost, 2000), climatic oscillations nonetheless affected the genetic diversity and evolution of many extant species in the region (Feng, Zheng, & Gong, 2016;Song et al., 2009;Zhou et al., 2010). Recent phylogeographic surveys in China indicate that during Pleistocene glaciations montane species might frequently have persisted in topographically complex regions that served as refugia, rather than being restricted only to lower latitudes Liu et al., 2013;Qiu, Fu, & Comes, 2011;Song et al., 2009;Zhou et al., 2010). As one of the major areas of endemism in China (López-Pujol et al., 2011), central China is characterized by a complex topography with altitudes ranging from 100 to 3,767 m. The relatively stable habitats within this region have enabled endemics to persist and also provided conditions for new species to form throughout Quaternary oscillations (López-Pujol et al., 2011;Ying, 1994;Ying & Zhang, 1993). However, in contrast to numerous phylogeographic surveys of plants from the Qinghai-Tibet Plateau (Liu et al., 2013;Luo et al., 2016;Zhang et al., 2018) and subtropical China (Fan et al., 2016;Wang et al., 2015) where the effects of climatic oscillations on genetic diversity have been well studied, relatively few phylogeographic analyses of plants in central China (Shao & Xiang, 2015;Zhou et al., 2010) have been conducted, and consequently, little is currently known about mechanisms responsible for the high level of plant endemism found there.
Tree peonies, which belong to section Moutan of the genus Paeonia (Hong & Pan, 2005;Stern, 1946), are distributed in central and northwestern China and include eight wild species (Hong, Pan, & Turland, 2001  Mountains (Hong, 2010;Li, 2011), although P. ostii is listed as "critically Endangered" in the Threatened Species List of China's Higher Plants (Qin et al., 2017) with extremely small populations owing to human overexploitation and domestication (Hong & Pan, 1999). Of the remaining species, Paeonia rockii is widely distributed throughout a large part of the Qinling-Daba Mountains and extends northwards in China at altitudes between 1,100 and 2,800 m, while Paeonia jishanensis and P. qiui have more restricted distributions, being confined to the northern and southern Qinling-Daba Mountains (Figure 1a), respectively, at altitudes between 750 and 1,700 m (Hong et al., 2001;Xu, Cheng, Xian, & Peng, 2016).
The three tree peony species are closely related showing minor morphological differences (Hong et al., 2001). Paeonia qiui has leaves with nine leaflets that usually lack obvious purple-red color above and is very similar to P. jishanensis, making their separate identification extremely difficult . Recent molecular phylogenetic analysis has confirmed a close relationship among the three species which are assigned to subsection Vaginatae (Zhang, Wang, Xia, & Zhou, 2009;Zhao, Zhou, Lin, Pan, & Li, 2008) with P. qiui and P. jishanensis considered as the most recently diverged sister species pair (Zhang et al., 2009;Zhao et al., 2008;Zhou et al., 2014). In phylogenetic trees produced from chloroplast DNA sequences, placement of P. qiui/P. jishanensis and P. rockii indicates that introgression may have occurred between them . A recent coalescent analysis using approximate Bayesian computation (ABC) of nuclear simple sequence repeat (nSSR) variation indicated that P. jishanensis and P. rockii diverged from their common ancestor ca. 140K-220K years ago (Yuan, Cornille, Giraud, Cheng, & Hu, 2014). Paeonia qiui was not included in this analysis and information on its phylogeography is study contributes to a growing understanding of the origin and population structure of tree peonies and provides insights into the high level of plant endemism present in the Qinling-Daba Mountains of Central China.

K E Y W O R D S
ecological niche modeling, genetic divergence, multiple refugia, niche divergence, phylogeography, speciation, tree peony therefore lacking, while that for P. jishanensis and P. rockii is limited.
The distributions of all three tree peony species have recently been decreasing based on our own field observation (Cheng, 2007;Xu et al., 2017;Yuan et al., 2012). Consequently, they have been listed as rare and endangered species in the Red Book of Chinese Plant Species (Fu, 1992). In addition, P. qiui and P. rockii have been classified as "endangered" in the Threatened Species List of China's Higher Plants, and P. jishanensis is listed as "vulnerable" (Qin et al., 2017). Hence, developing conservation programs for them has become urgent.
Here, we report an investigation of morphological divergence and evolutionary history of the three peony species, P. qiui, P. jishanensis, and P. rockii, based on nSSR and chloroplast DNA sequence variation. In addition, we examine niche divergence between the species using ecological niche modeling (ENM) and compare spatial distributions at present with those at the Last Glacial Maximum (LGM). Our research aimed to (a) investigate patterns of interspecific and intraspecific genetic differentiation among the three species; (b) estimate times of divergence; (c) test modes of origin and infer the demographic history of each species; and (d) consider effective conservation and management strategies separately for each of the three peony species in Central China.  (Table S1) distributed over five linkage groups according to the tree peony high-density genetic linkage map (Cai, 2015;. Details of nuclear microsatellite genotyping are described in Wu, Cai, Cheng, Cui, and Zhou (2014). In addition, three chloroplast (cp) DNA intergenic spacers (petB-petD, accD-psaI, and psbE-petL) were sequenced across a subsample of 296 individuals (Table 1) using the methods described in Grivet, Heinze, Vendramin, and Petit (2001) and Suo et al. (2012). All cpDNA sequences have been deposited in GenBank under accession numbers of KY200296-KY200338.

| Morphological differences
Six floral traits (color of petal, carpel, filament, stigma and flare on the petal base, and number of carpels) and seven vegetative characters (leaf color, type of compound leaves, stolon, lobed leaflet, shape of tip leaflet, number of leaflets, and plant height) were selected for recording (Table S2A). Each trait was recorded on 39 individuals across all three species (Table S2B). Differences for each trait among species were evaluated using the Kruskal-Wallis multiple-range test in Agricolae (De, 2014) and Vegan (Oksanen et al., 2013). In addition, principal component analysis (PCA) was performed to examine multivariate differences between species. All statistical analyses were performed in R (R Core Team, 2015).
Phylogeographical structure was estimated by testing whether N ST was significantly larger than G ST using PERMUT 1.0 with 10,000 permutations (Pons & Petit, 1996).
TA B L E 1 Geographic locations and sample sizes of 40 P. jishanensis, P. qiui, and P. rockii populations used in this study for seed plants (Graur & Li, 2000) was applied to calculate divergence times. Parameters were sampled every 1,000 steps for 10 7 MCMC steps, with the first 10% of samples discarded as burn-in.
Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) tests implemented in ARLEQUIN 3.5 were used to determine that all mutations were selectively neutral and to evaluate the hypothesis of demographic expansion. Demographic processes were examined using mismatch distribution analysis in ARLEQUIN 3.5. Multimodal mismatch distributions of pairwise differences between haplotypes indicate that a population size remained relatively stable over time, whereas a unimodal distribution implies a recent demographic expansion (Fu, 1997;Tajima, 1989). The validity of the estimated demographic model was tested by determining the distribution of the sum of squared devi-
First, population structure was examined using STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). For all populations of three peony species, we computed the likelihood for K clusters from 1 to 20 with ten replicates per K. For each species, ten replications were conducted for each K in the range K = 1-10.
Each run had a burn-in of 100,000 steps, followed by 200,000 MCMC iterations using the admixture model and assuming correlated allele frequencies. The optimal K was evaluated according to the methods of Pritchard et al. (2000) and Evanno, Regnaut, and Goudet (2005). The estimated admixture coefficients across replicated runs were permuted using CLUMPP v.1.1.2 (Jakobsson & Rosenberg, 2007). Graphics were produced using DISTRUCT v.1.1 (Rosenberg, 2004). To detect genetic grouping further, principal component analysis (PCA) was also conducted for all and each species implemented with GenALEx 6.5. An unrooted neighborjoining tree was also generated based on the shared-allele distance using PowerMarker 3.25 (Liu & Muse, 2005) to construct genetic relationships among individuals. To examine the partition of nSSR variation within and between species, hierarchical analysis of molecular variance (AMOVA) with 10,000 permutations was conducted using ARLEQUIN 3.5. We estimated pairwise F ST values between populations for all peony species in FSTAT 2.9.3 software (Goudet, 2001). This was followed by an evaluation of isolation by distance (IBD) between populations by testing the correlation between the matrix of pairwise F ST /(1-F ST ) and the matrix of geographic distances using GENALEX 6.5 with 999 permutations.
The directional relative migration network implemented in divMigrate (Sundqvist, Keenan, Zackrisson, Prodöhl, & Kleinhans, 2016) was used to estimate potential past gene flow among three peony species.

| Demographic model test with ABC
To compare plausible scenarios of divergence and to infer historical parameters of genetic divergence among the three species, we used DIYABC 2.0 (Cornuet et al., 2014) on the microsatellite dataset. Different evolutionary scenarios were constructed, taking account of the results obtained from the morphological divergence and STRUCTURE analyses. In step 1, six possible scenarios among the three species were tested; these differed in the relationships among taxa ( Figure 2a). Subsequently, in step 2, two competing scenarios based on the previous best-fit scenario that P. qiui diverged from a common ancestor with P. jishanensis were compared ( Figure 2b). In each scenario, current population sizes of the three taxa P. jishanensis, P. qiui, and P. rockii were denoted as NJ, NQ, and NR, respectively, while A1 represented the population size of the ancestral population at t2, and A2 represented the population size of ancestral population at t1 ( Figure 2). Because individuals of tree peony start flowering 5 to 8 years after germination and vegetative reproduction by rhizomes is very common (Sang, Crawford, & Stuessy, 1997), generation time was set at 10 years for all three species in our models.
We assumed a uniform prior distribution for mean mutation rate (μ) from 10 −4 to 10 −3 , which is consistent with rates for woody plants (Petit & Hampe, 2006). Summary statistics used for analysis were mean number of alleles, mean genetic diversity, F ST , and (dμ) 2 distance. The posterior probability of the competing scenarios was compared using a polychotomous logistic regression (Fagundes et al., 2007).

| Ecological niche modeling
To predict the geographic distribution of suitable habitat for each species, ecological niche modeling was conducted with MAXENT v.
F I G U R E 2 Evolutionary scenarios tested for the origin of P. jishanensis (PJ), P. qiui (PQ), and P. rockii (PR) using approximate Bayesian computation (ABC). In the analysis, current population sizes of the three taxa PJ, PQ, and PR were denoted as NJ, NQ, and NR, respectively, A1 represented the population size of the ancestral population at t2, and A2 represented the population size of ancestral population at t1. The size of an ancestral population was varied in all scenarios. (a) Six scenarios for the origins and relationships among three peony species. (b) Two competing scenarios based on the best-fit scenario of PQ from PJ. Scenario 1 predicted that PQ split from PJ at time t1, while PR split from PJ at time t2. P. jishanensis is the ancestor of P. rockii. Scenario 2 predicted that PQ split from PJ at time t1, while PJ split from PR at time t2. In scenario 2, P. rockii is the ancestor of P. jishanensis The same eight bioclimatic variables were used for the Community Climate System Model (Collins et al., 2010). We ran 20 replicates with 25% of the occurrence points used for model testing and determined logistic probabilities for the output. A presence-absence map was produced using the "maximum training presence threshold." The area under the curve (AUC) was used to evaluate model performance.

| Morphological divergence
ANOVAs and multiple-range tests ( Table 2; Table S5) showed that there were significant differences between species for all floral and vegetative traits recorded except plant height. Greatest divergence in terms of the proportion of total variance due to differences between species was evident for flare of the petal base (95%), carpel color (90%), filament color (90%), stigma color (90%), leaf color (80%), and numbers of carpels (61%;

| Divergence time and demographic history based on cpDNA data
Analysis of cpDNA haplotype variation using BEAST detected two clades with relatively high support (posterior probability, PP = 0.84; Tajima's D and Fu's Fs estimated for each species were non-significantly positive (Table 5), and the observed mismatch distributions of pairwise nucleotide differences were multimodal ( Figure S2). The sum of squares deviations (SSD) between the observed and expected distributions and the Harpending's raggedness index (H Rag ) were mainly significant (Table 5), indicating each species is in approximate demographic equilibrium.

| Nuclear microsatellite diversity and population structure
The frequency of null alleles at all loci was very low for P. jishanensis (0.051), P. qiui (0.087), and P. rockii (0.089), respectively. qiui (with a mean of 1.7) than for P. jishanensis (with a mean of 0.2) and P. rockii (with a mean of 0.7; Table 3).
STRUCTURE analysis of the complete dataset indicated that the highest peak for ΔK was at K = 2, followed by K = 3 ( Figure S4).
When K = 2, P. jishanensis and P. rockii individuals clustered into two separate group, while some individuals of P. qiui were admixed ( Figure 5a). Approximately 85% of the mean genetic composition of P. qiui came from P. jishanensis, while the remaining 15% came from P. rockii. However, when K = 3, P. jishanensis, P. qiui, and P. rockii were clearly assigned to different genetic groups, with each group representing a different species. Furthermore, the PCA of all individuals identified three groups (Figure 5b) as did the neighbor-joining tree constructed for all samples (Figure 5c). When a separate STRUCTURE analysis for each species was conducted, two genetic groups were identified for P. jishanensis and for P. qiui, whereas three genetic groups were found for P. rockii ( Figure S5).
The directional relative migration networks as estimated by divMigrate ( Figure 7) revealed a high degree of potential bidirectional gene flow between P. jishanensis and P. qiui, and between P.

| Mode of divergence and speciation
Results of our ABC analysis of alternative model scenarios for divergence and speciation based on nSSRs data are summarized in Table 6. In step 1, scenario 1 (Figure 2a) had the highest posterior probability with a 95% confidence interval not overlapping those of the other scenarios. In step 2, the same scenario was the best supported (PP = 0.6104), indicating that P. qiui split from P. jishanensis at time t1, while P. rockii split from P. jishanensis at time t2 (Figure 2b).
Our evaluation of confidence in scenario choice revealed low type I error and type II error, suggesting that our model choice analyses were reliable.
Results from the ABC analysis also indicated that the current population sizes of P. jishanensis (NJ) and P. qiui (NQ) decreased to 9.21 × 10 2 and 3.29 × 10 3 from a larger ancestral population size of 9.52 × 10 3 (A2), respectively ( Figure S6). In contrast, the current population size of P. rockii (NR) is 5.42 × 10 3 , which is slightly larger than the ancestral population size of A1 ( Figure   S6). Importantly, the median divergence time between P. jishanensis and P. rockii was estimated to be 50.  Figure S6).

| Ecological niche modeling
Ecological niche modeling revealed that predictions for the present distributions of P. jishanensis and P. rockii were mostly congruent with observed distributions (Figure 8). However, for P. qiui it was evident that in some predicted areas the species is not currently present (e.g., northeast of Sichuan

| Species identification and cytoplasmic-nuclear discordance in Paeonia
Our analyses based on morphological and nSSR data support the view that P. jishanensis, P. qiui, and P. rockii are distinct species. In agreement with previous morphological analyses (Zhao et al., 2008;Zhou, Pan, & Hong, 2003), our morphological analysis showed that the three peony species can be distinguished according to floral and leaf characters. In addition, they are highly divergent based on an analysis of nSSR variation. Thus, a hierarchical AMOVA of the nSSR dataset revealed significant genetic differentiation between the three species, while PCA and a neighbor-joining tree placed P. jishanensis, P. qiui, and P. rockii individuals into different clusters representing the three species.
Finally, ENM and multivariate niche space analysis identified significant ecological differentiation between the three species. STRUCTURE analysis of the complete dataset suggested the presence of only two genetic groups, with some P. qiui individuals indicated to be admixed, implying potential gene flow between P. jishanensis and P. rockii, which was also supported by the directional relative migration network as estimated by divMigrate. The latter analysis indicated high potential rates of gene flow between all three peony species. However, when K = 3, P. qiui individuals were assigned to a separate genetic group.
A greater genetic similarity of P. qiui to P. jishanensis than to P.
rockii was made further evident by the pairwise F ST comparisons between species. The possibility that P. qiui diverged from P. jishanensis was supported by ABC analysis of nSSR data. Despite P.
qiui having a more restricted distribution than P. jishanensis and P.
rockii, it exhibited the highest levels of genetic diversity (at both cpDNA and nuclear microsatellite markers). Furthermore, P. qiui contained the highest number of private alleles in most populations based on the nSSR analysis.
F I G U R E 4 BEAST-derived trees among the three peony species based on cpDNA haplotypes TA B L E 5 Results of neutrality tests for regional and mismatch distribution analysis for three peony species In contrast to the clear separation of species according to morphological and nSSR variation, the cpDNA haplotype network showed that the 18 haplotypes resolved did not cluster by species and that P. qiui shared two haplotypes (H2 and H12) with P. rockii. Such cytoplasmic-nuclear discordance in species delimitation commonly occurs in plants and has also been reported in other tree peonies of the group P. delavayi/ludlowii (Zhang et al., 2018) and may reflect the effects of incomplete lineage sorting or introgressive hybridization (Comes & Abbott, 2001;Wu & Campbell, 2005). Because the effective population size of nuclear DNA is four times that of cpDNA (Palumbi, Frank, & Hare, 2001;Schaal, Hayworth, Olsen, Rauscher, & Smith, 2010), lineage sorting should be faster for cpDNA (Guillemin et al., 2016;Qu et al., 2012). Thus, the most plausible explanation for the discordance in the present study is that hybridization and introgression occurred between the species, particularly between P. qiui and P. rockii, at times of contact in the Daba and central eastern Qinling Mountains during and/or after the LGM. There is good evidence for reticulate evolution having occurred frequently in Paeonia (Sang, Crawford, & Stuessy, 1995;Sang & Zhang, 1999;Yuan, 2010), and our results support this for P. qiui and P. rockii and further emphasize that biparentally inherited SSR markers are more effective than maternally inherited cpDNA markers in delimiting these species.

| Genetic differentiation among the three tree peony species
High levels of population differentiation were observed within the three species for nSSR data and especially for cpDNA (Table 4), in comparison with the average values for plant species (Nybom, 2004 andPetit et al., 2005). Similar levels of high genetic differentiation have been reported in P. rockii before (cpDNA: Yuan, Cheng, & Zhou, 2011;nSSRs: Yuan et al., 2012), P. jishanensis (nSSRs:  and in another tree peony, P. delavayi (nSSRs and plastid DNA: Zhang et al., 2018). In contrast, levels of genetic differentiation recorded among the three closely related peony species for both markers (F STnSSR = 0.498 and F STcpDNA = 0.929) were slightly higher than those detected between two other peony species, P. delavayi and P. ludlowii (F STnSSR = 0.403 and F STplastid DNA = 0.420), endemic to the Himalayan-Hengduan Mountains (Zhang et al., 2018). Although explanations for high levels of genetic differentiation have been sufficiently addressed in previous studies for P. rockii (Yuan et al., 2011(Yuan et al., , 2012 and P. jishanensis , relatively little is known about the factors affecting genetic divergence among the three species (Sang et al., 1995;Yuan et al., 2014).
In this study, such significant genetic differentiation indicates that despite evidence of introgressive hybridization having occurred in the past, barriers to gene flow are very strong  (Luo, Pei, Pan, & Hong, 1998;Yuan et al., 2011). Similarly, gene flow via seed dispersal will be limited as seeds are normally dispersed by gravity or by rats (Hong, 2010) and when China experienced two abrupt climatic warming phases during the last glacial cycle (Shi, Zhao, & Wang, 2011;Zhao, Shi, & Jie, 2011).
This estimate of time of divergence between these two species is more recent than that reported previously (140-220 Ka; Yuan et al., 2014). The difference in estimated age likely results from differences in mutation rates used in the ABC computation. In the present study, we set mutation rates to 10 −4 and 10 −3 , which is consistent with rates for EST-SSRs in woody plants (Petit & Hampe, 2006), producing an estimated median mutation rate of 1.74 × 10 −4 . In contrast, Yuan et al. (2014) set mutation rates to 10 −3 and 10 −2 , yielding a much higher and what we regard as less accurate median mutation rate of 1 × 10 −3 .
The most likely model of divergence and speciation also suggested that P. qiui diverged from P. jishanensis ~24.3 Ka, which is at the time of the LGM (18-25 Ka), and that its origin was not a consequence of hybridization between P. jishanensis and P. rockii. Our results indicate that Quaternary climate oscillations had an important effect in driving speciation of these sister species. Although Central China has never been directly impacted by extensive and unified ice-sheets (Qiu et al., 2011), it nonetheless experienced climatic oscillations throughout the Quaternary with dramatic F I G U R E 6 Plot of geographical distance against genetic distance for population comparisons of (a) all three species, (b) P. jishanensis, (c) P. rockii, and (d) P. qiui effects on the evolution and distribution of species (Harrison, Yu, Takahara, & Prentice, 2001;Qiu et al., 2011). Quaternary climate oscillations and associated environmental changes promoted range fragmentation and population isolation, thus providing opportunities for allopatric speciation (Comes, Tribsch, & Bittkau, 2008;Qiu et al., 2011). According to our field observation, P. jishanensis and P. qiui have different distributions endemic to the north and south of the Qinling mountain chain, respectively. Most species in temperate areas in Central China displayed a fragmented distribution pattern during the glacial cycles of the Quaternary (Qian & Ricklefs, 2000;Hewitt, 2000;. By comparing the relative importance of climatic variables on the distributions of the two species, we found that P. jishanensis and P. qiui show preferences for different climatic conditions. As a consequence, the geographic distributions of P. jishanensis and P. qiui likely became fragmented, thus promoting conditions for allopatric divergence among isolated populations during Quaternary climate oscillations. Our estimates of divergence times should be treated with caution due to the wide 95% confidence intervals attached to them. In addition, vegetative reproduction by rhizomes among tree peonies may prolong generation time, causing us to underestimate divergence time (Sang et al., 1995). Nonetheless, the estimated time for divergence of all three tree peony species approximates to the most recent Pleistocene glaciation, which if correct suggests that their origin was not triggered by the rapid and massive uplift of the Qinling Mountains estimated to have occurred between 1.2 and 2.4 Ma (Yuan et al., 2012), but more likely was triggered by climate change during the late Pleistocene and the effects this had on the distribution of ancestral forms at that time.
Climate change during late Pleistocene glacial-interglacial cycles is well known to have dramatically influenced species distributions (Comes & Kadereit, 1998;Hewitt, 1996;Klicka & Zink, 1997 TA B L E 7 Posterior median estimate and 95% highest posterior density interval (HPDI) for species effective population sizes, t1, t2, and μ in the best supported scenario 2 in step 2 of the ABC analysis, based on nSSRs data for three tree peony species to higher latitudes and altitudes during interglacials (Davis & Shaw, 2001;Hewitt, 1996). Reductions in genetic diversity have often been detected in higher latitude populations due to loss of such diversity during the recolonization process (Hewitt, 1996;Nason, Hamrick, & Fleming, 2002). However, in the present study, for both nSSRs and cpDNA, we did not detect significant declines in genetic diversity among populations of the three peony species with increasing latitude.
Based on cpDNA variation, none of the three species exhibited a strongly unimodal mismatch distribution or significant negative values of Tajima's D or Fu's F S , which suggests that they maintained relatively stable population sizes throughout the last glacial period and into the current interglacial. However, a different picture with regard to the stability of population sizes of P. jishanensis and P. qiui emerged from the ABC analyses of nSSR variation. This indicated that the population sizes of these two species were much larger during the LGM than those currently, indicating that the two species experienced historical expansions of population size throughout the LGM rather than contraction.
This hypothesis is further supported by the results of ecological niche modeling, which indicated larger distribution ranges for both P. jishanensis and P. qiui at the LGM compared to currently. In contrast, the population size of P. rockii (NR) did not greatly alter, indicating that it may have persisted in situ and maintained a stable population size throughout the last glacial period, as suggested by the results of analyses of neutral tests and mismatch distribution of its cpDNA variation. This mode of in situ persistence during the Pleistocene glacial periods is also supported by the phylogeographic studies of other peony species, P. ludlowii and P. delavayi, in the Hengduan Mountains, although P. delavayi is believed to have shown both in situ persistence and retreat to refugia during these periods, Zhang et al., (2018). In contrast with the general "contraction-expansion" scenario (Hewitt, 1996), our analyses reveal population expansion or local persistence throughout the LGM, rather than contraction. Such a phylogeographic pattern of expansion or stability during the last glacial period has also been highlighted in a few recent studies (Fan et al., 2016;Liu et al., 2013;Xie et al., 2012).
The unusual demographic history of the three tree peony species might be partly explained by the climatic characteristics of the Qinling-Daba Mountains. The east-west orientation of this long mountain range in central China forms a natural division between the current northern temperate and subtropical regions (Qian & Ricklefs, 1999. Topographically diverse landscapes in the Qinling-Daba Mountains can buffer regional climate variability and therefore create stable climatic conditions (Hewitt, 2004;Tang et al., 2018;Tzedakis, Lawson, Frogley, & Hewitt, 2002) for these peony species to maintain large population sizes in this region during the last glacial period and current interglacial, respectively.
In general, refugia are characterized by climatic stability allowing species to survive through glacial and/or interglacial periods (Tang et al., 2018;Tzedakis et al., 2002). In addition, populations in refugia usually harbor not only relatively higher genetic diversity but also haplotypic uniqueness (Petit et al., 2003;Zhou et al., 2010). In our study, as revealed by our analysis of nSSR and cpDNA variation, potential and into the Holocene. In addition, Shennongjia probably was the main refugium, because in this region population SBL harbors the highest number of private alleles. Indeed, the stable refugia identified for P.
rockii are consistent with previous finding that populations of P. rockii located in the western Qinling Mountains harbored relatively higher genetic diversity as revealed by both cpDNA and nSSR markers (Yuan et al., 2011(Yuan et al., , 2012. Our results therefore support the hypothesis that the Qinling-Daba Mountains with a variable and complex geography could provide stable conditions for species to survive in this region, therefore fostering the accumulation of genetic diversity during the LGM (Tzedakis et al., 2002). The Qinling-Daba Mountains have long been regarded as refugia for fauna and flora during the Quaternary (Dong et al., 2011;Ying, 1994), as indicated by other studies on plants Cathaya argyrophylla (Wang & Ge, 2006), Saruma henryi , and members of the Chrysanthemum indicum complex (Li et al., 2014), and also in some animals, for example, Alcippe morrisonia (Song et al., 2009) and Feirana taihangnica (Wang, Jiang, Xie, & Li, 2013). In the case of Paeonia jishanensis, the possible refugial population (JYA) is located in the Taihang Mountains of the northern Qinling range, which is also believed to have served as a refugium for the plants Opisthopappus longilobus and O. taihangensis at the time of the LGM (Wang & Yan, 2014). In summary, our study suggests that multiple refugia occurred in Qinling-Daba Mountains, which allowed the three peony species to persist in situ during the Pleistocene glaciations. The current distribution pattern of refugia is similar to that of other regional endemics including Paeonia delavayi distributed in Hengduan Mountains (Zhang et al., 2018) and temperate plants in Central and South China (Qiu et al., 2011) such as Castanopsis tibetana and Schima superba (Fan et al., 2016), with multiple refugia sustained across their distribution ranges. Stable climate refugia located in Qinling-Daba Mountains have long been recognized as likely centers of plant endemism in Central China (López-Pujol et al., 2011;Ying, 1994Ying, , 2001 and foci of speciation (Tang et al., 2018). Consequently, they should be protected against deforestation and habitat destruction in order to conserve diversity at the present time and into the future.

| Conservation implications
Knowledge of the level and structure of genetic diversity in endangered species is an important element in designing conservation programs (Ouborg, 2010  The eight bioclimatic variables (bio2, bio3, bio4, bio5, bio9, bio13, bio15, and bio19). **p < 0.001 for correlations between PC axes and geographical variables.
However, the distribution of P. rockii has recently been decreasing due to habitat destruction and genetic fragmentation Yuan et al., 2012). Therefore, exploitation for medicinal use and anthropogenic habitat destruction should be strictly prohibited to preserve the remnant populations of these species.
In addition, because ex situ conservation can be critical to protecting endangered species (Heywood & Iriondo, 2003), we advise the establishment of a germplasm resource nursery from seeds and rhizomes collected from populations of P. jishanensis and P. qiui that harbor high levels of genetic diversity.

ACK N OWLED G M ENTS
The authors are sincerely grateful to the editor, several anonymous reviewers, and to Kermit Ritland for helpful comments on previous versions of this manuscript. We thank Wei Cheng, Ye Yuan, and Ling Yu for their generous assistance in collecting plant materials.
We also thank our friends from the Forestry Department of Shaanxi

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

AUTH O R CO NTR I B UTI O N S
X-XXu, J-FM, and F-YC conceived the overall study. X-XXu performed the experiments. X-XXu, Y-QS, X-GH, K-HJ, and L-PP analyzed the data. X-XXu, J-FM, F-YC, and RA wrote the manuscript. All authors read and approved the final manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The three cpDNA sequences have been deposited in GenBank (accessions numbers KY200296 and KY200338).