Phylogeography and genetic effects of habitat fragmentation on endemic Urophysa (Ranunculaceae) in Yungui Plateau and adjacent regions

Urophysa is a Chinese endemic genus with only two species (U. rockii and U. henryi) distributed in Yungui Plateau (Guizhou Province) and adjacent regions (i.e., Provinces of Hunan, Hubei, Chongqing and Sichuan). The aim of this study was to determine the genetic diversity and population differentiation within Urophysa and investigate the effect of the Yungui Plateau uplift and climate oscillations on evolution of Urophysa. In this study, micro-morphological characteristics, nine microsatellite loci (SSR), two nuclear loci (ITS and ETS) and two chloroplast fragments (psbA-trnH and trnL-trnF) were used to analyze the phylogenetic relationships and assess genetic and phylogeographical structure of Urophysa. Isolation by distance (IBD) was performed to research the effects of geographical isolation. We detected high genetic diversity at the species level but low genetic diversity within populations. Striking genetic differentiation (AMOVA) among populations and a significant phylogeographical structure (NST > GST, p < 0.01) were detected among U. henryi populations, along with significant effects of isolation by distance (IBD). Molecular clock estimation using calibration strategy and cpDNA substitution rate indicated that the divergence of U. henryi occurred during late Miocene to early Quaternary, when the orogeny of Yungui Plateau was violent. U. rockii originated at the early Quaternary and further differentiated at early Pleistocene. Our results suggested that habitat fragmentation played an important role in the genetic diversity and population differentiation of U. rockii and U. henryi. Heterogenous geomorphological configuration and complicated environment resulted from rapid uplift of the Yungui Plateau were inferred as important incentives for the modern phylogeograhpical pattern and species divergence of Urophysa. The geographical isolation, limited gene flow, specialized morphologies and the Pleistocene climatic oscillation greatly contributed to the allopatric divergence of U. rockii. Significant genetic drift and inbreeding were detected in these two species, in situ measures should be implemented to protect them.


Introduction
Geological history and climate oscillations are important drivers in the evolution and genetic structure of plant species [1,2]. Previous researches have indicated that landscape heterogeneity caused by the uplift of the Qinghai-Tibetan Plateau (QTP) greatly contributed to the evolution of many plants (i.e., Aconitum gymnandrum, Taxus wallichiana and Rhodiola kirilowii) in the QTP [3][4][5][6]. However, it is uncertain how landscape heterogeneity from uplift in adjacent areas (to the QTP) affected plant evolution [7,8], such as the Yungui Plateau. The Indian plate has kept moving northward since the Cenozoic, colliding with the Eurasian Plate between 40 and 50 million years ago (Mya). As a result, the orogeny of the QTP is violent, which induced the continuing uplift of the Yungui Plateau, contributing to the formation of its unique geomorphological characteristics and heterogeneous landscape environment [9]. The heterogeneous landscape was typically characterized by low soil water content, periodic water deficiency, and poor nutrient availability, which exert strong selective forces on plant evolution, resulting in remarkably high species richness and endemism in the Yungui plateau [10]. Underlying this species differentiation of the plateau, there was a massive divergence in population genetics and a promotion of ecological diversity. For example, the Yungui Plateau and its adjacent regions have been regarded as an important center of origin for the East Asiatic flora [11]. In addition, the unique characteristics of the plateau have given rise and refuge to a variety of endemic plants [12], and now is a refuge for plants that are threatened elsewhere.
Orogeny always resulted in regions of most plants shifted and divided populations into different spatial-temporal scales, along with breaking up of one patch of habitat into several smaller patches [13,14]. This phenomenon was called habitat fragmentation, which plays a major role in threatening plant species survival. Habitat fragmentation is often considered to be the main driving force for local species extinction [15,16] as large populations split into smaller populations with increasing geographic isolation [17]. The spatial isolation of populations may restrict connectivity, resulting in low levels of gene flow between fragments, with subsequently lower genetic diversity and higher genetic differentiation in/among remnant populations [1,18,19]. Habitat fragmentation can reduce the fitness of remnant plants by affecting population genetics and dynamics in several ways [19], including genetic erosion [20] population divergence and random genetic drift [1,21]. Smaller isolated populations often experience greater inbreeding, reduced reproduction rate and offspring survival [22]. Conversely, genetic variation could be maintained or even increased in fragmented populations [23,24], and habitat fragmentation may play an important role in plant divergence and allopatric speciation [1].
Urophysa Ulbr. (Ranunculaceae) is an Chinese endemic genus with two species, U. rockii and U. henryi. These two species are morphologically distinct: the former has sacs near the base of its petals while the latter has no sac (Fig 1). Additionally, U. rockii has a more restricted distribution (only located in Jiangyou county, Sichuan Province) than U. henryi (distributed in Yungui Plateau, south Chongqing, north Hunan and west Hubei). The two species' natural populations are restricted to small and isolated areas separated by high mountains and deep valleys, and grow in steep and karstic cliffs. The Yungui Plateau is bordered to the east by the Mountains (Mts.) Daba, Wuling and Xuefeng, separating U. henryi Yungui Plateau populations (populations in Guizhou Province) from neighboring populations (populations in Hunan, Hubei and Chongqing Provinces). Moreover, these two species have strict habitat requirements and are sensitive to environmental change and human activities (scenic spots, power station). By field observations and laboratory experiments, we found that U. rockii and U. henryi can not survive once leaving the karst limestone, and populations that located in hydroelectric dams and tourist attractions possess less individuals than that have not been effected by human activities. The plants are collected for Chinese traditional medicine for the treatment of contusions and bruises. However, the population of the two species are in decline [25,26]. There is ongoing research on the endangered U. rockii, including its growing environment and conservation strategies [25], biological and ecological characteristics [27,28] and reproductive biology [29]. However, there have been no studies on the evolutionary history and distribution of these two species. Similarly, species differentiation and the effect of habitat fragmentation in the Yungui Plateau have been little studied. Few studies that have been conducted have found that the occurrence of physical barriers in the Yungui Plateau and its adjacent regions have affect the genetic structure of species such as Eurycorymbus cavaleriei, Saruma henryi and Dipentodon [30][31][32]. Therefore, we determined the phylogenetic relationship between the two species, their lineage and cause of speciation and why they exhibit such special distribution and distinct morphologies, particularly for the endangered U. rockii.
Phylogeographical approaches are useful for reconstructing evolutionary histories of species and their close relatives [33], and are reliable for examining interspecific divergence and speciation [1]. Here, we used nine nuclear microsatellites (nSSR), two nuclear gene fragments (ITS and ETS) and two cpDNA regions (psbA-trnH and trnL-trnF) to investigate the phylogeography structure of U. rockii and U. henryi, and determine their divergence and evolution. It is commonly believed that chloroplast genomes are susceptible to the influence of fragmentation because they are maternally inherited, which can only be dispersed via seeds in sexual reproduction. This is in contrast to nuclear regions, that are diploid and can be dispersed via seeds or pollen. Therefore, we considered that the cpDNA markers were good indicators of population historical dynamics and genetic effects [34]. Our specific aims were to: (1) Evaluate the genetic diversity and population differentiation of U. rockii and U. henry. (2) Investigate phylogeographical patterns of Urophysa and (3) Investigate the effect of the QTP uplift and climate oscillations on the evolution of Urophysa, especially U. rockii. We believe that this study will be useful for developing conservation strategies and restoration of these two endangered U. rockii and endemic U. henry.

Materials and data
Plant sampling and morphological observation A total of 190 individual plants from 14 populations were collected, covering almost the entire geographic range of both species (Table 1), Nine to 16 individual plants were sampled from the 14 populations in the Yungui Plateau (Guizhou Province) and its adjacent regions (i.e., Provinces of Hunan, Hubei, Chongqing and Sichuan Provinces). Healthy leaves were sampled and dried in silica gel until total DNA was extracted. Voucher specimens were deposited in the Herbarium of Sichuan University (SZ). The detailed population information of sampled individuals was measured using a handheld GPS unit. Morphological characteristics were identified using herbarium specimens and fresh materials. The mature pollen grains and leaves dehydrated by graded ethanol were directly mounted on aluminum stubs using conducting carbon adhesive tab, sputter-coated with gold, and then observed using the JSM-7500F scanning electron microscope (SeM, Japan).

DNA extraction, sequencing and microsatellite genotyping
Total genomic DNA was isolated from silica-gel dried leaves using plant genomic DNA kit (Tiangen Biotech, Beijing, China). Internal transcribed spacer (ITS) sequences were amplified using the primers ITS4 and ITS5 [35] and ETS fragments were amplified based on the primers ETS9bp and ETS18s [36]. The primers trnH and psbA [37]were used to amplify the psbA-trnH sequences, and the trnL-trnF sequences were amplified by trnL and trnF [38]. Polymerase chain reaction (PCR) applications were carried out in 30μL reaction volumes. For the ITS and ETS regions, reactions were conducted with the following program: an initial 4-min denaturation at 94˚C followed by 30 cycles of 45-sec denaturation at 94˚C, 45-sec annealing at 56˚C for ITS (53˚C for ETS) and 90-sec extension at 72˚C with a final 5-min extension at 72˚C [35]. For the psbA-trnH and trnL-trnF intergenic spacers, the PCR program began with 4-min initial denaturing at 94˚C followed by 30 cycles of 1-min denaturation at 94˚C, 1-min annealing at 54˚C for psbA-trnH (45 sec at 64˚C for trnL-trnF), and 2-min extension at 72˚C. A final extension was run for 7 min at 72˚C for psbA-trnH (5 min at 72˚C for trnL-trnF) [38]. The PCR products were separated in 1.5% (w/v) agarose TAE gel and purified using Wizard PCR preps DNA Purification System (Promega, Madison, WI, USA) following the manufacturer's instructions. The purified PCR products were sequenced in an ABI Genetic Analyzer (Applied Biosystems Inc., Foster City, CA, USA) in both directions using the PCR primers. The chloroplast DNA sequences and nuclear fragments were edited and assembled using seqMan (DNAstar) [39], then aligned employing CLUSTAL W in MEGA 5 [40] and adjusted manually. All haplotype sequences were deposited in the GenBank database under accession numbers KR820593 to KR820702 (S1 Table). According to the primers and amplification protocols developed for Li et al. [41], we selected nine pairs of SSR primers for our population genetic data analysis (S2 Table). A preexperiment was performed using 20 pairs of EST primers, of which, nine pairs of primers that showed significant polymorphism and high homology with relatives Aquilegia Li et al. [41] were employed in following analyses (S3 Table). PCR products were separated on 3.5% of agarose gel followed by staining with ethidium bromide. Alleles were sized using PeakScanner v. 1.0 software (Applied Biosystems).

Data analyses
Genetic diversity and divergence. The DnaSP version 5.0 [42] was used to identify different haplotypes and to calculate haplotype diversity (H d ) and nucleotide diversity (π). To assess the level of genetic variation and population differentiation, the average within populations gene diversity (H S ), total gene diversity (H T ) and two population differentiation parameters, G ST [43] and N ST [44] were estimated following the methods described by Pons and Petit [45] using PERMUT [45]. The analysis molecular variance (AMOVA) in Arlequin version 3.5 [46] was used to further quantify genetic differentiation between groups or subgroups, as well as between populations within groups and among individuals within populations. In addition, maximum parsimony median-joining method in NETWORK 4.2.0.1 [47] was applied to construct haplotype networks.
Fst is a measure of checking population differentiation and inferring the effects of gene flow and drift [48]. Through constructing the regional scatterplots of the genetic differentiation (quantified as Fst/(1 -Fst)) on geographical distances and calculating the correlation coefficients between them, we can evaluate the relative historical influences of gene flow and genetic drift on regional population structure [49]. Based on this method, the Fst values between pairwise populations were calculated using Arlequin version 3.5 [46]. Scatterplots of the Fst/(1 -Fst) on geographical distances were constructed, and correlation coefficients were calculated along with the significance of correlation in GenAlEx version 6.0 [50], The linearized Fst statistic [Fst/(1 -Fst)] was compared with the matrix of geographical distance by means of a simple Mantel test to detect isolation by distance and to evaluate the relative influences of gene flow and genetic drift on the regional population structure. We used 999 random permutations to test for the Mantel statistical significance.
For SSR analysis, the total number of detected alleles (N A ), allelic richness (A R ) and inbreeding coefficient (F IS ) were calculated using the software fstat 2.9.3.2 [51]. The MICRO-CHECKER v. 2.2.3 [52] was used to test for evidence of deviations from the Hardy-Weinberg equilibrium, genotyping errors and null alleles, and compute the expected heterozygosity (H E ) and observed heterozygosity (H O ). A cluster analysis was performed to characterize population structure in STRUCTURE version 2.2 [53], along with 100 000 Markov chain Monte Carlo (MCMC) cycles following 10 000 burn-in cycles, using the admixture model with independent allele frequencies. Eleven replications were performed for each K, in the range K = 2-10, and the optimal K value was determined according to Evanno et al. [54]. In addition, the analysis molecular variance (AMOVA) in Arlequin version 3.5 [46] was used to measure the partitioning of genetic variability within and among populations. In addition, the IBD (isolation by distance) analysis was performed in GenAlEx version 6.0 [50].
Phylogenetic analysis. Maximum parsimony (MP) and Bayesian inference (BI) analysis were conducted using the program PAUP Ã version 4.0 beta 10 [55] and MrBayes v3.1.2 [56], respectively. For parsimony analyses, heuristic searches were carried out with 1000 random sequence replicates, with the tree bisection-reconnection (TBR) branch swapping and the Mul-Trees options selected. All characters were weighted and unordered equally and gaps were treated as missing data. Branch supporting values were estimated with 1000 replicates by bootstrap analysis. For Bayesian inference analyses, the best-fit evolutionary model was selected by the Akaike information criterion (AIC) using MrModeltest 2.2 [57]. The Bayesian Markov chain Monte Carlo (MCMC) searches was performed for 2 × 10 −7 generations with four chains, sampling trees every 1000 generations and discarding the first 20% of sampled trees as burn-in sample. The remaining trees were used for constructing a 50% majority-rule consensus tree and calculating posterior probability (PP). In addition, maximum likelihood (ML) analyses were performed by RAXML v7.2.8 [58] with 1000 bootstraps under the GTR + G substitution model, Nodes with a bootstrap value of ! 95% were considered well-supported in this analysis. Additionally, we made a model test for each locus, and performed the phylogenetic analysis again both in RAxML and MrBayes.
Demographic history analysis. Mismatch distribution analysis was undertaken to assess whether genealogy experienced historical population expansions. It is assumed that if populations experienced a sudden demographic expansion should display a unimodal and smooth distribution, and the multimodal distributions are related to demographic equilibrium or decline [59]. Gene flow between populations (N M ) and Neutrality test including the Tajima's D test [60] and Fu and Li's D Ã and F Ã statistics [61] were calculated by Arlequin version 3.5 [46].
Taking advantage of fossils from Prototinomiscium vangerowii (Menispermaceae) [62] and calibration points used in previous studies of Aquilegia [63], we estimated the divergence times of Urophysa. Bayesian searches for tree topologies and node ages of cpDNA and nrDNA dataset were conducted in BEAST using a GTR + G substitution model selected by JMODELTEST [64] and an uncorrelated lognormal relaxed clock [65]. A Yule process was specified as tree prior and the calibration priors was modelled as normal distributions with a mean time. Unfortunately, no fossils are known for Urophysa or any closely related lineages. Therefore, we chose the age of the fossil Prototinomiscium vangerowii (Menispermaceae). The Menispermaceae are represented in the Cretaceous of Europe by endocarps assigned to the fossil genus Prototinomiscium. Based on the oldest record, from the Turonian of Central Europe, dated to 91.0 (Mya) [62], sets the minimum age of the split between Menispermaceae and Ranunculaceae [66]. In addition, we employed the age interval 51-66 Mya reported by Wikström et al. [67] for the split of Ranunculus (subfamily Ranunculoideae) and Xanthorhiza (subfamily Coptidoideae), the former genus being in the sister subfamily to Thalictroideae, which includes the genus Urophysa [68]. Based on this information, we set priors of 91.0 ± 7.5 and 58.0 ± 2.5 Mya respectively for each calibration point. A Yule speciation model was selected as the tree prior. This is a simple model of speciation that is more appropriate when considering sequences from different species. Considering Bastida et al. [63], we dated the origin of Thalictroideae as a mean age of 27.61 Mya with a normal distribution and with a specified 95% confidence interval (CI) of 26.59-28.56 Mya. The stem node of Aquilegia was dated to 10.18 Mya with a normal distribution and with the specified 95% confidence interval (CI) of 9.21-11.14 Mya [63]. The program BEAST version 1.8.0 [69] was employed to estimate the divergence of major lineages of Urophysa, MCMC runs were performed, each of 2 × 10 7 generations with sampling every 2000 generations, following a burn-in of the initial 10% cycles. MCMC samples were inspected in Tracer to confirm sampling adequacy and convergence of the chains to a stationary distribution. All of the outgroup sequences used for time calibration were listed in S4 Table. The substitution rates were also used to estimate time of divergence of Urophysa. Aquilegia incurvata and Semiaquilegia adoxoides were chosen as outgroups [63,70]. Considering the dominant geographical distribution of cpDNA haplotypes and poorly resolved relationships among nrDNA haplotypes, the cpDNA dataset was used to estimate the divergence time. The best-fit model was GTR + G for cpDNA data inferred from the Akaike information criterion (AIC) with MrModeltest 2.2 [57]. An uncorrelated lognormal model was selected to describe the relaxed clock. The constant cpDNA substitution rates for most angiosperm species have been estimated to be in the range 1-3 × 10 −9 s/s/y [71]. Given Urophysa are perennial plants, the general substitution rates of the plastid sequence (u = 1.52 × 10 −9 s/s/y) was more suitable. The Markov chain Monte Carlo (MCMC) analyses were run for 2 × 10 7 generations and trees were sampled every 2000 generations. The first 10% of sampled trees was discarded as the burn-in sample as checked with Tracer 1.5 [69]. The Tree Annotator version1.4.8 [69] was used to summarize the samples in the maximum clade credibility tree with the posterior probability limit set to 0. 5

cpDNA sequence analysis
The two cpDNA regions were aligned along a total length of 1567 bp (psbA-trnH, 572 bp; trnL-trnF, 995 bp) and 17 chlorotypes were generated, including 71 polymorphisms (S5 Table). Chlorotypes H1-H14 were found only in U. henryi and H15-H17 were found in U. rockii (Fig 2 and Table 2). Only three chlorotypes (H1, H3, H15) were shared by two or more populations and no chlorotype was shared by these two species. In addition, we detected no haplotype located in the center of the chlorotype network, while many private cpDNA haplotypes were detected in each population. At the species level, for U. rockii, H d = 0.742 and π = 0.00117, for U. henryi, H d = 0.917 and π = 0.01398. Population JY2 of U. rockii and SM of U.
henryi possessed the highest haplotype diversity and the maximum number of haplotypes. The total genetic diversity (H T ) was higher than the diversity within populations (H S ) for both species, and H T and H S were higher in U. henryi than in U. rockii. Additionally, N ST was significantly higher than G ST (N ST = 0.921, G ST = 0.716, P < 0.01) but only in U. henryi (Table 3), indicating a significant phylogeographical structure existing between populations of U. henryi.
AMOVA indicated that 41.32% of the total cpDNA variation was partitioned between species and 57.06% of the variation could be attributed to variation between populations within species (Table 4). For each species, statistically significant variation was detected between populations (88.24% for U. rockii and 97.30% for U. henryi). The observed value in the mismatch distribution analysis fitted multimodal curves (S2A-S2E Fig) and the results of Tajima's D and Fu and Li's D Ã and F Ã statistics were not significantly negative (S6 Table), which indicated that the populations did not undergo expansion. Mantel test results indicated a significant effect of isolation by distance (IBD) at the genus range scale (r = 0.322, p = 0.001), as well as the species range of U. henryi (r = 0.246, The topological structure of haplotypes derived from Parsimony analyses was similar to that from Bayesian tree and maximum likelihood tree analyses, therefore only the Bayesian tree with maximum parsimony and maximum likelihood bootstrap support values was shown in S4 Fig    haplotypes was consistent with the strict consensus tree produced by Bayesian analysis ( Fig  2C). In addition, after testing the substitution model of each locus and combining to perform phylogenetic analysis, we found that the topology of RAxML and MrBayes trees were similar, and was consistent with topology which was obtained from cpDNA haplotypes (S7 Table and  S4 Fig).
From two time-calculated ways (calibration strategy and substitution rate), we estimated the divergence time of Urophysa, and found that the calibrated time estimated by cpDNA data is more accurate and reliable compared with time obtained by substitution rate. The origin time of Urophysa was dated to approximately 10.29 Mya (95% HPD = 8.99-11.69 Mya) (Fig 3).

nrDNA sequence analysis
The length of aligned sequences was 665 bp for ITS, and 521 bp for ETS. Based on 69 nuclear polymorphic sites combined the two data sets, 34 haplotypes (H1-H34) were detected (Fig 4A  and S8 Table). All of the haplotypes were unique to a population except H33 (Table 2), which was shared by population JY3 and JY4. The haplotype diversity (Hd) ranged from 0 to 0.378 and the nucleotide diversity (π) from 0 to 0.00051 within U. rockii. The haplotype diversity (Hd) within U. henryi ranged from 0 to 0.978 and the nucleotide diversity (π) from 0 to 0.00765. Diversity at the species level for U. rockii was H d = 0.756 and π = 0.00287, and for U. henryi was H d = 0.943 and π = 0.00930. The highest haplotype diversity and the maximum number of haplotypes were observed in population SM of U. henryi and JY2 of U. rockii. AMOVA showed that 62.14% of the genetic variation occurred between species and 31.98% of the variation was partitioned among all populations within species. Within both species, variation between populations was significant and accounted for 82.06% and 53.11% of total variation in U. henryi and U. rockii, respectively (Table 4). Similar to the findings for cpDNA, mismatch distribution analysis was clearly multimodal (S2F- S2H Fig), as well as a positive and non-significant result of neutrality test (p > 0.10) (S6 Table), which indicated the populations did not undergo expansion.
A high total genetic diversity (H T = 0.989) and low average within-population genetic diversity (H S = 0.397) was detected in U. henryi and in U. rockii (H T = 0.900, H S = 0.116). Similar to cpDNA, the results of nrDNA (N ST > G ST , P < 0.01, Table 3) indicated a significant phylogeographical structure existed between populations of U. henryi. Mantel test results indicated a significant effect of isolation by distance (IBD) at the genus range scale (r = 0.301, p = 0.004), whereas no significant IBD pattern presented in the two species (S3D-S3F Fig).
Only the Bayesian tree with parsimony bootstrap and maximum likelihood support values was shown in S6 Fig because it has the same topology as the maximum parsimony tree and maximum likelihood tree. U. henryi and U. rockii were clustered into a single clade with high bootstrap support and posterior probability values (S6 Fig). The network of nrDNA haplotypes was roughly consistent with the phylogenetic tree, but no ancestral haplotypes were detected (Fig 4B).
Divergence time estimation by nrDNA data (ITS) indicated that Urophysa began to diverge at approximately 9.65 Mya (95% HPD = 7.48-11.48 Mya) (S7 Fig), which was slightly older than that by cpDNA data. Estimations of divergence times of species within Urophysa are not reliable because of low branches posterior probabilities in the nrDNA tree, and it is better to use the time estimated by cpDNA data.  (Table 5). To test the impact of these null alleles, we removed the microsatellite markers (A41, EST2, EST9) with null alleles in both species and used the remaining six markers to perform the population structure analysis. We found that excluding microsatellite markers with null alleles had little impact on our results (S8 Fig Table 5). All pairwise population genetic differentiations were significant (p < 0.001), except between adjacent U. rockii populations JY3 and JY4 (pairwise F ST = 0.054, p > 0.05). Pairwise F ST estimates are listed in S9 Table. F IS values ranged from 0.103 to 0.930 across the U. rockii with an average value of 0.476. F IS values ranged from -0.148 to 1.000 for U. henryi, with an average value of 0.384.
The STUCTURE analysis, using the ΔK method, showed that the optimal K value was K = 2 (Fig 5A), which strongly supported two genetic clusters among our samples and generally corresponded to the two respective species. When K = 5, we detected a small peak, which further divided populations of U. henryi into four groups (Fig 5). AMOVA indicated that 31.7% of the genetic variation occurred between species (F CT = 0.317, p < 0.001). Within each species, most of the genetic variation was partitioned among populations (Table 4). A significant effect of isolation by distance (IBD) was detected between the 14 populations of Urophysa (r = 0.646,

Genetic diversity and significant population differentiation
Based on cpDNA and nrDNA data sets, high haplotype diversity (H d ) and nucleotide diversity (π) were observed in U. henryi and U. rockii at the species level (Table 2). We detected a high level of total genetic diversity in U. henryi (H T = 0.959 for cpDNA, H T = 0.989 for nrDNA) and in U. rockii (H T = 0.891 for cpDNA, H T = 0.900 for nrDNA), which was higher than other species of Ranunculaceae, such as Aconitum gymnandrum (H T = 0.739 for cpDNA) and Clematis sibirica (H T = 0.496 for nrDNA) [4,73]. Genetic diversity of species is closely related to environmental factors (e.g., climate, topography) and life history characteristics (i.e., life cycle, breeding system). U. henryi and U. rockii possess a relatively long life cycle and sexual reproductive period and high offspring mortality, which may lead to low levels of genetic diversity [74]. It has been found that habitat fragmentation can reduce genetic diversity due to restricted gene flow, genetic erosion and random genetic drift in isolated populations [15,19,75,76]. However, other studies have shown that genetic variation could be maintained or even increased in fragmented populations [23,24]. Individuals of U. henryi and U. rockii are confined to steep and karstic limestones in ravines, which are naturally fragmented. These conditions are found especially in the Yungui Plateau, as high mountains, deep valleys and fast-flowing rivers acting as natural barriers enhancing isolation, drift and mutation [31]. Isolation of U. henryi and U. rockii individuals was also supported by the IBD model and Neutrality test (S3 Fig, S6 and S9 Tables), indicating a significant isolation-by-distance pattern between populations of these two species and populations of U. henryi. Thus, the high genetic diversity observed in these two species may be closely related to their fragmented habitats. Despite a high level of cpDNA genetic diversity detected at the species level, the average genetic diversity within-population was low (H S ). We found a high coefficient of genetic differentiation ( Table 4). These results indicate high genetic differentiation exists between populations. Akin to our study, the combination of low genetic diversity and high genetic differentiation has been reported for other plants in Yungui Plateau and adjacent regions, such as Cercidiphyllum japonicum, Tetrastigma hemsleyanum and Cardiocrinum giganteum [77][78][79].
Limited gene flow may be a crucial factor resulting in low genetic diversity and high genetic differentiation. It is believed that the breeding system of plant (i.e., flowering and seed proliferation) is an important biological characteristic, which strongly influences the spatial-temporal distribution of genetic variation [80] and affects processes that lead to speciation or extinction [81]. Although the reproduction system and pollination mechanism of U. henryi has been little studied, previous research of U. rockii reproduction has indicated the seeds are extremely small and are dispersed into rock gaps due to mechanical strain when the follicle cracks spontaneously after seeds have matured [29]. However, most seeds are washed away by rainwater, and only a few seeds fall to the base of cliff where they may germinate. Petit et al. [82] found that species with seeds dispersed by gravity tended to show higher differentiation between populations than species with wind-dispersed seeds. Interestingly, U. rockii and U. henryi are entomophilous plants and bloom in the winter when lower temperatures reduce the visiting frequency of pollinators, which may result in weak pollen flow. The poor gene flow mediated by seeds and pollen was also supported by our results (N M ) ( Table 5 and S6 Table). Thus, limited gene flow in such fragmented habitats could lead to significant genetic differentiation between populations [1].
Human activities have significantly influenced the distribution of these two species. Over the past few decades, many hydroelectric dams and tourist attractions have been built in areas where the wild populations of U. henryi and U. rockii are located leading to the extensive removal of their habitat. The number of individuals has also sharply declined due to excessive collection for their medicinal value. All of these human activities have reduced the population size, increased fragmentation and isolation, and enhanced population differentiation.

Demographic history of U. henryi
The time of origin obtained from the molecular clock estimation is generally congruent with that from points calibration. The molecular clock estimated that the populations of U. henryi in the Yungui Plateau (Clade I) diverged approximately at 8.86 (6.24-11.0) Mya, which is in accordance with the rapid uplift time of the QTP [83]. It is believed that the QTP reached an elevation similar to the present at about 8.0 Mya, but decreased following extensive faulting. The most recent rapid uplift of the OTP occurred at around 3.6 Mya [84]. Given the geologic close relationship between the Yungui Plateau and the QTP, the orogenic events of the Yungui Plateau were similarly violent during the continuous uplift of the QTP [85]. Therefore, we suggest that rapid QTP uplift and subsequent Yungui Plateau contortions significantly contributed to the differentiation of U. henryi. Besides, global climate also fluctuated dramatically at that time and the prevailed East Asian monsoon brought plenty of rainfall [86].
In the cpDNA tree, two clades of U. henryi were revealed (Clade I and clade II) (Fig 3 and  S4 Fig), which were consistent with geographical distribution of U. henryi (Fig 2). This grouping was reflected by a striking differentiation between populations. U. henryi populations were further divided into four groups when K = 5 in structure analysis, and this roughly corresponded with geographical regions. Mountains such as Wuling, Dalou and Xuefeng. extend from northeast to southwest with an average altitude of 2,000 m [87], and acted as geographical barriers between Yungui Plateau (Clade I) and its adjacent populations (clade II) of U. henryi. Previous phylogeographical studies have identified that Mt. Wuling and Mt. Xuefeng as the major barriers to gene flow [31,88]. In addition, the significant increase in geological and ecological variabilities during rapid uplift of the Yungui Plateau has promoted rapid divergence in other small and isolated populations such as Ligularia-Cremanthodium-P arasenecio complex, Babina pleuraden, Dipentodon and Eurycorymbus cavaleriei [30,31,89,90]. Therefore, the habitat of U. henryi was fragmented due to geographical barriers and fluctuating climate conditions (i.e., unstable rainfall), and populations of U. henryi were separated causing high population differentiation.

The allopatric divergence of U. rockii
To investigate the phylogenetic relationship between U. rockii and U. henryi, we identified distinct morphological traits of petals, leaf epidermis and pollen grains (Fig 1 and S1 Fig). Results of SSRn and nrDNA (Figs 4 and 5, S6 Fig) showed that the U. rockii and U. henryi are distinctly separate species and the genus Urophysa was a monophyly, which were consistent with previous research [70]. The cpDNA phylogeny demonstrated that U. rockii was at the end of the cpDNA phylogenetic tree with high bootstrap support values and exhibited significant differentiation with U. henryi.
Habitat fragmentation is likely to significantly influence the divergence and allopatric speciation of plants [1]. We estimated the origin time of U. rockii was 3.16 Mya. At that time, orogeny of the Yungui Plateau triggered habitat fragmentation and significant geographic isolation between populations. U. rockii grows exclusively on cliffs or fissures of rocks in China and only a few populations have been found. Connectivity with U. henryi populations was hindered by geographic barriers, the Sichuan Basin and Yangzi River, and thus there was no gene flow (Fig 2 and S6 Table). The effect of barriers on gene flow in this region has been documented for Myotis pilosus [91] and Tapiscia sinensis [92]. The significant isolation caused U. rockii to undergo allopatric divergence. Additionally, We found considerable genetic differentiation and limited gene flow between U. rockii and U. henryi (Table 5 and S6 Table). This high genetic differentiation between populations is likely due to isolated populations with restricted gene flow [93]. This is also supported by the Mantel test between the pairwise F ST /(1-F ST ) and geographic distance (S3 Fig). Moreover, local adaptation may have played an important role in driving allopatric speciation of U. rockii. As expected, we identified extensive footprints of local adaptation from its specialized morphologies. These specialized morphologies included unusual floral organs such as petaloid sepals that can display various colors in different flower phases, the petals with a nectar spur that could attract more pollinators, and a mass of small seeds (thousand seed weight is 0.6684 ± 0.0038g) [29]. These specialized morphologies likely contributed enhanced reproductive efficiency.
The divergence of the U. rockii was dated to approximately 1.48 Mya, coinciding with the frequent climatic oscillation during the Pleistocene, which is considered as one of the most important periods for genetic diversification [2]. It is believed that dramatic climate fluctuation during the Pleistocene resulted in massive changes in plant population distributions, with many species shifting to more suitable habitat and other spatial-temporal adjustments [94]. Mismatch analysis indicated that expansion did not occur in populations of U. rockii (S2E and S2H Fig), and that U. rockii populations were presumably impacted by the climatic oscillation of the Pleistocene. However, the mountain system located in north of Sichuan Basin was in favor of preserving plant species [95] The mountains of Qinling, Daba and Micang could also have acted as barriers to reduce the effects of the Pleistocene climate [96]. These mountains have also been regarded as a key glacial refuge for other plants including Pinus massoniana, Liriodendron chinense and Rhinolophus ferrumequinum [79,95,97]. Overall, we believe that long-term geographical isolation, limited gene flow caused by habitat fragmentation and specialized morphologies probably contributed to the allopatric divergence of U. rockii. The climatic oscillation of the Pleistocene further promoted population divergence and resulted in the current distribution of U. rockii.

Implications for conservation
U. rockii is an endangered species of China that has a small geographic range, few distinct populations and extremely low germination rate (less than 2%) [29], and each population has a limited number of individuals (all five populations >2,000) [98]. Genetic drift is likely to have occurred, which could have led to the observed genetic diversity decline and increased population differentiation. The correlation estimation of the genetic differentiation [Fst/(1 -Fst)] and geographical distances (isolation by distance) showed that genetic drift was much more influential than gene flow on the distribution of genetic variability. Consequently, populations of U. rockii are not at equilibrium, which may have resulted from its strict habitat requirement or narrow distribution. A high value of F ST indicates high genetic drift load, while a low value of H S signifies high inbreeding [22,99]. In our study, genetic drift and inbreeding were both high ( Table 3 and Table 4), which was mirrored by the F IS from SSR data (Table 5), indicating a risk of extinction.
Although U. rockii and U. henryi demonstrated high genetic diversity at the species level, high genetic drift and inbreeding can be deleterious to species survival. Similarly, species with small population sizes and fragmented distributions are vulnerable to extinction especially with high genetic drift and inbreeding [100]. Protection of in situ populations and their habitat, including removing threats (i.e., human disturbance) is essential for U. rockii's survival in the wild. Ex situ breeding programs can also be implemented using seeds of local provenance to propagate different genotypes to ensure diversity and provide a source of seedling for planting in the wild. Although U. henryi is not (yet) endangered, frequent human activities have fragmented and degraded its habitat. It is possible that U. henryi will be threatened in the future due to persistent human disturbance, restricted range and isolated populations. It is necessary to adopt clear and practicable measures to restrict anthropogenic disturbances now, as genetic diversity of wild populations can be maintained and conserved, rather than implement measures after future declines.  Table. Information of our haplotype sequences deposited in the GenBank. Numbers in the brackets is the sequence number of each haplotype included. Ã : represent the specific haplotype of each population. (DOC) S2 Table. Microsatellite markers used in this study. For each primer pair, forward (F) and reverse (R) primer sequence, repeat motif (Repeat), size of cloned allele (bp), optimal PCR annealing temperature (Ta). (DOC)