Effects of Mountain Uplift and Climatic Oscillations on Phylogeography and Species Divergence in Four Endangered Notopterygium Herbs

Mountain uplift and climatic fluctuations are important driving forces that have affected the geographic distribution and population dynamics history of organisms. However, it is unclear how geological and climatic events might have affected the phylogeographic history and species divergence in high-alpine herbal plants. In this study, we analyzed the population demographic history and species differentiation of four endangered Notopterygium herbs on the high-altitude Qinghai–Tibetan Plateau (QTP) and adjacent areas. We combined phylogeographic analysis with species distribution modeling to detect the genetic variations in four Notopterygium species (N. incisum, N. franchetii, N. oviforme, and N. forrestii). In total, 559 individuals from 74 populations of the four species were analyzed based on three maternally inherited chloroplast fragments (matK, rbcL, and trnS-trnG) and one nuclear DNA region (internal transcribed spacer, ITS). Fifty-five chloroplast DNA (cpDNA) and 48 ITS haplotypes were identified in the four species. All of the cpDNA and ITS haplotypes were species-specific, except N. franchetii and N. oviforme shared one cpDNA haplotype, H32. Phylogenetic analysis suggested that all four species formed a monophyletic clade with high bootstrap support, where N. franchetii and N. oviforme were sisters. In addition, each Notopterygium species generated an individual clade that corresponded to their respective species in the ITS tree. Population dynamics analyses and species distribution modeling showed that the two widely distributed herbs N. incisum and N. franchetii exhibited obvious demographic expansions during the Pleistocene ice ages. Molecular dating suggested that the divergence of the four Notopterygium species occurred approximately between 3.6 and 1.2 Mya, and it was significantly associated with recent extensive uplifts of the QTP. Our results support the hypothesis that mountain uplift and Quaternary climatic oscillations profoundly shaped the population genetic divergence and demographic dynamics of Notopterygium species. The findings of this and previous studies provide important insights into the effects of QTP uplifts and climatic changes on phylogeography and species differentiation in high altitude mountainous areas. Our results may also facilitate the conservation of endangered herbaceous medicinal plants in the genus Notopterygium.


INTRODUCTION
Geological events and climatic fluctuations are considered to have profoundly shaped the distribution and population dynamics history of species in mountain areas (Hewitt, 2004;Hickerson et al., 2010). Thus, during glacial periods, most species experienced adverse weather conditions in high altitude mountains, where they contracted into refugia in low latitudes and then their ranges expanded again after the ice ages, thereby leading to species divergence or secondary contact evolution (Hewitt, 2004;Nybom, 2004;Ohsawa and Ide, 2008). In addition, the distribution patterns and population genetic structures of some species were reshaped due to geographic barriers and climatic oscillations. Moreover, the population size, mating system, and bio-characteristics of species had important effects on the divergence and evolutionary history of populations of species (Stewart et al., 2010;Zhang et al., 2015). For example, some studies have suggested that the locations of ice age refugia for plants were determined mainly by the adaptability of species to the external environment (Stewart et al., 2010;Liao et al., 2015;Wang et al., 2015).
During the Quaternary ice periods, many species experienced extinction events due to repeated bottlenecks and genetic drift, which led to further divergent evolution and genetic isolation within species (Soltis et al., 2006;Hickerson et al., 2010;Stewart et al., 2010;Keppel et al., 2012). Repeated environmental changes may also have promoted the fragmentation of habitats, as well as causing exotic distributions of different species or intraspecific genetic changes (Hickerson et al., 2010;Jia et al., 2012). Studies of alpine trees have shown that populations were crossed during mountain uplift processes, whereas the exchange of genes among populations was restricted due to climatic and geographic barriers (McLachlan et al., 2005;Soltis et al., 2006;Birks and Willis, 2008;Avise, 2009;Gao et al., 2012). Some alpine species experienced deep lineage divergence due to climatic changes and environmental isolations (Qiu et al., 2011;Liu et al., 2012;Xu et al., 2015).
In the high latitudes of Europe and North America, studies suggested that plant species could have survived in high elevation areas ("invisible refugia") during the ice age periods Provan and Bennett, 2008;Parducci et al., 2012;Stewart and Stringer, 2012;Ortego et al., 2012Ortego et al., , 2015Guichoux et al., 2013;Allen et al., 2015;Cavender-Bares et al., 2015). The presence of Juniperus species in the Qinghai-Tibetan Plateau (QTP) region also supports the existence of invisible refugia (Opgenoorth et al., 2010). The QTP is the largest and highest plateau in the world, with a mean altitude of more than 4000 m. Various endangered species and high levels of global diversity are present on this plateau (Mittermeier et al., 2005). Studies suggest that extensive uplifts of the QTP occurred in the Miocene-Pliocene era between 3.6 and 1.7 Mya (Li and Fang, 1999;Zhang et al., 2000;Zhou et al., 2006). The lifting of mountains triggered species divergence and changed the genetic structure to affect the evolution of high-alpine plants Wen et al., 2014;Ickert-Bond and Renner, 2016). In particular, the geological effects of the QTP on the genetic structure, geographic distribution, and species differentiation of plants have been clearly defined in this area (Wang et al., 2010;Xu et al., 2010;Liu et al., 2013;Sun et al., 2014;Favre et al., 2015;Hughes and Atchison, 2015). However, most of these previous studies focused on the response patterns of tree or shrub species to mountain uplifts and climatic oscillations on the QTP (Liu et al., 2006;Wang et al., 2009;Mao et al., 2010;Xu et al., 2010;Tian et al., 2011;Qiu et al., 2011;Wen et al., 2014;Favre et al., 2015;Hughes and Atchison, 2015), whereas little is known about the effects of mountain uplifts and climate events on cold-tolerant herbal species in the high altitude QTP and adjacent regions.
The genus Notopterygium H. de Boissieu (Apiaceae) comprises perennial and endangered herbaceous medicinal plants, which are mainly distributed in the QTP and its surrounding highaltitude areas. According to records in the Flora of China, this genus comprises six species: N. incisum C. C. Ting ex H. T. Chang, N. oviforme R. H. Shan, N. franchetii H. de Boissieu, N. forrestii H. Wolff, N. tenuifolium M. L. Sheh and F. T. Pu, and N. pinnatiinvolucellum F. T. Pu and Y. P. Wang. N. incisum and N. franchetii have wide distribution ranges at altitudes of 3200-5100 m and 1700-4500 m, respectively. N. oviforme occurs in the eastern part of the QTP at altitudes of 1700-3200 m. The other three species, i.e., N. forrestii (4000-4300 m), N. tenuifolium (4300 m), and N. pinnatiinvolucellum (3400 m), have very limited distributions among the highalpine shrubs and meadows in the west region of China. These herb species provide an excellent model for detecting the effects of the QTP uplifts and Quaternary climatic oscillations on the genetic structure and species divergence of plants. However, in recent years, due to high market demand, the wild resources of these Notopterygium species have decreased rapidly because of human over-exploitation (Zhou et al., 2010). The Notopterygium species are now listed as endangered herb species in the IUCN Red List, and their management and conservation are urgently required (Wu et al., 2005). Information regarding geographic distributions and genetic diversity is vital for formulating effective conservation strategies for wild plant resources. However, most of the previous studies of the Notopterygium species have focused mainly on their phylogenetic evolutionary relationships (Pu et al., 2000;Yang et al., 2017), morphological and physiological characteristics (She and Pu, 1996;Wang et al., 1996;Jiang et al., 2005), and comparative transcriptome analysis (Jia et al., 2017), whereas little is known about their genetic divergence and population demographic history.
In the current study, we sampled four species, i.e., N. incisum, N. oviforme, N. franchetii, and N. forrestii, across their entire geographic distributions in the high-altitude QTP and adjacent areas. We detected the genetic variations in three chloroplast DNA (cpDNA) markers and a nuclear DNA fragment in order to characterize the population histories and species divergence of these endangered herb plants. Our aims were: (1) to determine the genetic structure and population evolutionary history of four Notopterygium species; (2) to identify the phylogenetic relationships among these species and their phylogeographic history; (3) to explore the effects of QTP uplifts and climatic changes in the Quaternary on the divergence and phylogeography of these species; and (4) to propose reasonable conservation FIGURE 1 | Geographic distribution of cpDNA haplotypes for the four Notopterygium species. Each circle represents a population and each color represents each haplotype. The colored outlines of the circles distinguish the four species, where green indicates N. incisum, yellow indicates N. franchetii, blue indicates N. oviforme, and red indicates N. forrestii. and management strategies for the endangered Notopterygium species.

Sample Collection
In this study, in order to obtain information about genetic variation over a wide area, 559 individuals from 74 populations were collected for the four Notopterygium species in Sichuan, Shaanxi, Gansu, Qinghai, and Shanxi provinces in the highaltitude QTP and adjacent areas. These samples covered the complete geographic distribution ranges of the four species in the QTP and surrounding areas. From 2 to 17 individuals were sampled from each population, where all of the samples collected were separated from each other by at least 100 m. Detailed information about the latitude, longitude, and altitude for all of the populations is provided in Figure 1 and Table 1. All of the materials and documents have been deposited in the College of Life Sciences, Northwest University. In addition, two species from the genus Pleurospermum, i.e., P. prattii and P. franchetianum, as well as Heracleum moellendorffii were used as outgroups.

DNA Extraction and Sequencing
Total DNA was extracted using the modified CTAB method (Doyle and Doyle, 1987) or with a plant DNA extraction kit (Tiangen, Beijing, China). We used 1% agarose gels to check the quality of the DNA extracted from the Notopterygium species. To screen for suitable primers, we first randomly selected 50 individuals (one individual from each population) to amplify the universal cpDNA primers and nDNA primers recommended by the Consortium for the Barcode of Life (CBOL) (CBOL Plant Working Group, 2009). Finally, three highly variable cpDNA primers, i.e., trnS-trnG, matK, and rbcL, and one nDNA internal transcribed spacer (ITS) primer were selected to determine the genetic variations in the genus Notopterygium after initial tests with six loci (Supplementary Table S1). Thus, two monomorphic cpDNA loci (trnL-trnF and rpl36-infA; Bai et al., 2010;Soumaya et al., 2014) were excluded from all of the subsequent analyses.
PCR amplification was performed in a volume of 25 µL containing 2 µL DNA template (10-50 ng/µL), 12.5 µL PCR MIX (Xi'an Runde, China), 0.75 µL of each primer (20 ng/µL), and 9 µL double-distilled H 2 O. The PCR reaction conditions were as described in Supplementary Table S2. All of the high quality PCR products were sequenced using the amplified forward and reverse primers with an ABI 3730 XL genetic analyzer (Applied Biosystems, Foster City, CA, United States). All of the

Proofreading and Alignment of DNA, and Data Analysis
BioEdit v 7.0.9.0 (Hall, 1999) software was used for manual proofreading and checking the variable sites. MEGA v 7.0 (Sudhir et al., 2008) was used to remove low quality sequences and only high quality sequences were analyzed. For the ITS sequences, we visualized the possible color spectrum of the overlapping peaks at any one variable site. If a strong signal peak was more than half of a weak signal peak, then we used the strong peak for phrasing. If both the peaks overlapped, we used the following phrases for each variable site instead of both peaks as phrases: R: A+G, Y: C+T, M: A+C, K: G+T, S: G+C, W: A+T. DnaSP v 5.0 software was used for dividing the heterozygous loci into double sequence series (Librado and Rozas, 2009).

Genetic Variation and Genetic Structure Analysis
The genetic diversity of the cpDNA and ITS sequences were analyzed in all four Notopterygium species using PERMUT v 1.0 software, where we calculated the genetic diversity within the population of each species (h S ), total genetic diversity (h T ), and population genetic differentiation coefficients G ST and N ST (Grivet and Petit, 2002). In addition, ARLEQUIN v 3.5 (Excoffier and Lischer, 2010) software was used to perform analysis of molecular variance (AMOVA) for the cpDNA and ITS sequences. AMOVA partitioned the genetic differentiation among the populations F ST , within a population F SC , and among species F CT .

Phylogenetic Analysis
Phylogenetic analyses of the cpDNA and ITS sequences were performed with MEGA v 7.0. JModeltest v 3.06 (Posada and Crandall, 1998) was used to filter the best evolutionary model (GTR+G). One-thousand bootstrap replicates were performed for the maximum likelihood (ML) and maximum parsimony (MP) models to obtain the phylogenetic tree. MrBayes v 3.2.3 was also used to conduct phylogenetic analyses of the cpDNA and ITS sequences based on the Bayesian criterion (Huelsenbeck and Ronquist, 2001). We set the random tree rotation as 10,000,000 generation, where each 1000 generations were kept to construct a phylogenetic tree, with a burn-in of 2500.
NETWORK v 5.0.0 (Polzin and Daneshmand, 2003) was used to construct median-joining networks of the cpDNA and ITS sequences. ArcGIS v 10.2 (Bader, 2005) was employed to draw the haplotype distribution map. BEAST v 1.7.5 (Drummond and Rambaut, 2007) was used to estimate the divergence times of the cpDNA haplotypes where we used the cpDNA evolutionary rates (1.0-3.0 × 10 −9 s/s/y) recorded for other angiosperms to calibrate our datasets due to the lack of fossil evidence for Notopterygium plants (Wolfe et al., 1987). We employed the loose molecular clock method with an uncorrected lognormal distribution for the branch lengths. After a burn-in of 5,000,000 steps, all of the parameters were collected once every 1000 steps up to 50,000,000 Markov chain Monte Carlo (MCMC) algorithm steps. The convergence of the MCMC results was verified by using the Tracer v 1.5 program to check that the chain was balanced, where we then used the Tree Annotator v 1.7.5 program to obtain the best tree merging and Figtree v 1.3.1 (Rambaut, 2009) was employed to view the resulting tree.

Population Dynamics Analysis
DnaSP v 5.0 was used to analyze the genetic diversity parameters, including the haplotype diversity (H d ) (Nei and Tajima, 1981), nucleotide diversity (π) (Nei and Li, 1979), and number of haplotypes (H). We also used DnaSP v 5.0 to detect the mismatched distributions (Schneider and Excoffier, 1999) of cpDNA sequences in the four Notopterygium species. Population demographic expansions were tested using Arlequin v 3.5 (Excoffier and Lischer, 2010) and Tajima's D (Tajima, 1989), Fu's F S (Fu, 1997), and Fu and Li's F * (Fu and Li, 1993) tests. We used the sum of the squared deviations between the observed and expected mismatches as well as Harpending's raggedness index values (Rag) (Harpending, 1994) to determine the validity and significance level of the expansion model. According to the formula: τ =2ut (τ is the mismatch equilibrium expansion   variable) (Rogers and Harpending, 1992), we calculated the expansion time t, where u is the mutation rate per generation calculated using the formula u = 2µkg, where µ is the mutation rate per nucleotide per year, k is the total length of a cpDNA sequence, and g is the generation time. According to our field investigations, the generation time for Notopterygium species was 3 years. In order to further determine the signs of demographic growth in the four Notopterygium species, we used LAMARC v 2.1.8 (Kuhner, 2006) to calculate the population growth parameter g. The MCMC algorithm was run for 100,000 generations and sampled every 200,000 steps, where the first 25% of the sampled trees were discarded as the burn-in.

Species Distribution Modeling
We used MaxEnt v 3.3.3k (Phillips et al., 2006;Phillips and Dudík, 2008) to predict the current, last glacial maximum (LGM), last interglacial (LIG), and future distributions of two widespread Notopterygium species: N. incisum (148 distribution sites) and N. franchetii (80 distribution sites). The distribution sites of Notopterygium species were collected from previous studies as well as websites containing climate data and plant distributions. We also obtained some distribution sites based on field investigations. Bio-climatic environment data were downloaded from the WorldClim website 1 at a resolution of 2.5 arc-minutes. Six bioclimatic environmental variables (Supplementary Table S7) with significant effects on N. incisum and N. franchetii were used to detect changes in the distribution ranges of plants. We set the number of replicates to 10 and the maximum number of iterations to 500 for MaxEnt modeling. The accuracy of the model's performance was assessed based on the area under the receiver operating characteristic curve (AUC) (Fawcett, 2006).

RESULTS cpDNA Variations and Haplotype Distributions
Three chloroplast fragments (matK, rbcL, and trnS-trnG) were used to analyze 559 individuals from 74 populations of the four Notopterygium species. The total length of the fragments was 1605 bp, and the lengths of the matK, rbcL, and trnS-trnG regions were 669, 668, and 268 bp, respectively, which included 21, seven, and eight nucleotide mutation sites (Supplementary Table S3). The cpDNA regions were uniparental inherited markers so we combined the three chloroplast fragments in the subsequent population genetics analysis.

ITS Sequence Variation
The total length of the sequenced ITS region was 593 bp and 48 haplotypes were identified with 66 nucleotide mutation sites (Figure 2 and Supplementary Table S4). All of the ITS haplotypes were species-specific in the four Notopterygium species. The total haplotype diversity (H d ) and π values for N. incisum, N. franchetii, and N. oviforme were 0.71 and 0.0025, 0.55 and 0.00364, and 0.69 and 0.0024, respectively. N. incisum populations from the southeast part of the QTP (E, G, H, J, U, W, X, HB and HC) had the highest haplotype diversity in this species, and haplotypes H1 and H7 had the highest distribution frequencies. N. franchetii populations from KG, KI, KJ, KK, KP, YC, YD, YE, and YF had the highest haplotype diversity in this species, and haplotype H22 had the highest frequency. In addition, N. oviforme and N. forrestii exhibited low haplotype diversity in terms of their ITS sequences ( Table 2 and Supplementary Figure S1).

Genetic Diversity and Structure
The total genetic diversity (h T ) values based on the cpDNA datasets for N. incisum, N. franchetii, N. oviforme, and N. forrestii were 0.939, 0.766, 0.961, and 0.623, respectively, where N. incisum had the highest levels for h S (0.404) and h T , whereas N. forrestii had the lowest level of diversity (h S = 0.167; h T = 0.623) ( Table 3). In addition, we calculated the genetic differentiation coefficients G ST and N ST for the four species. The U statistic (Gaussian test 1000 times) showed that N ST was significantly larger than G ST for N. incisum and N. oviforme (P < 0.05), thereby indicating that these two species exhibited significant phylogeographic structuring (Table 3).
AMOVA analysis of the cpDNA datasets detected genetic variations among the four species (F CT = 0.5804) (Supplementary Table S5 Table S5). In addition, the AMOVA results obtained for the ITS sequences indicated similar genetic differentiation patterns to those based on the cpDNAs, where the differences among species in terms of the variation in the ITS were as high as 92% (F CT = 0.9287) (Supplementary Table S6).

Phylogenetic Relationships
Phylogenetic trees of the cpDNA haplotypes were constructed based on the ML, MP, and Bayesian inference methods, which showed that the topological structures obtained were basically the same using the three methods (Figure 3). The four species of Notopterygium formed a larger monophyletic clade with high bootstrap support, where N. franchetii and N. oviforme were sisters. The median-joining network diagram produced using the cpDNA datasets was consistent with the phylogenetic analysis (Figure 2). The major haplotypes with the highest distribution frequencies (H1, H11, H32, H33, H41, and H53) were located in the central positions of the network. However, the phylogenetic relationships of the ITS sequences differed from those of the cpDNA sequences. No haplotypes were shared among the four Notopterygium species and each species formed its own individual branch (Figure 3) in the ITS tree. In addition, in order to confirm the phylogenetic positions of all four species considered in this study, we analyzed the other two species in the genus Notopterygium, i.e., N. tenuifolium and N. pinnatiinvolucellatum. Phylogenetic analyses based on variations in the chloroplast rbcL sequence showed that the four species considered in this study, i.e., N. incisum, N. franchetii, N. oviforme, and N. forrestii, were more closely related than N. tenuifolium and N. pinnatiinvolucellatum (Supplementary Figure S2).

Population Dynamics History and Divergence Time
Based on the cpDNA sequences, we performed various mathematical analyses to determine the population histories of the four Notopterygium species (Table 4 and Supplementary Figure S3). The mismatch distribution model had a single peak, with negative Tajima's D and Fu's F S values for N. incisum and N. franchetii, which suggested that these two species had experienced rapid range expansions. The larger population growth indexes for N. incisum (g = 809) and N. franchetii (g = 2810.736) were also consistent with rapid population expansions. By contrast, N. oviforme and N. forrestii had bimodal mismatch distributions with positive Tajima's D and Fu's F S values, where these results indicated that they did not experience expansion events. Therefore, we estimated the expansion times for N. incisum and N. franchetii as about 128-43 Kya and 51-17 Kya in the Pleistocene, respectively ( Table 5).

Species Distribution Modeling
In this study, MaxEnt modeling had the highest predictive capacity (AUC > 0.9) for the two widely distributed Notopterygium species (N. franchetii and N. incisum). The distribution ranges predicted for these two species were consistent with the current geographic distributions in the QTP and adjacent areas (Figures 5, 6 and Supplementary  Table S7). Species distribution modeling also showed that the range of N. incisum was limited in the LIG period whereas it expanded very rapidly in the LGM period. However, there were no significant changes in the distribution range from the LGM until the current period for this species. For N. franchetii, MaxEnt modeling suggested that the distribution range of this species increased very rapidly from the LIG until the LGM period. However, it was interesting that the distribution range of N. franchetii did not change greatly from the LGM until the current period.
FIGURE 4 | Chronogram for the four Notopterygium species obtained using BEAST based on the plastid sequences. The turquoise color bar indicates the 95% highest posterior density (HPD) credibility intervals for node ages (million years ago, Mya). Posterior probabilities are labeled above the line, and the mean divergence dates and 95% HPDs are labeled below the line.

DISCUSSION Genetic Diversity and Structure
In the current study, our analysis of the cpDNA sequences showed that N. oviforme had the highest level of genetic diversity (H d = 0.81, π = 0.0013), followed by N. incisum (H d = 0.75, π = 0.00086) and N. forrestii (H d = 0.39, π = 0.0002), whereas N. franchetii had the lowest level of diversity (H d = 0.29, π = 0.00031) ( Table 2). However, the results were different according to the ITS sequence analysis, where N. incisum had the highest diversity (H d = 0.71, π = 0.0025), followed by N. oviforme (H d = 0.69, π = 0.0024) and N. franchetii (H d = 0.55, π = 0.0036), whereas N. forrestii exhibited no variation ( Table 2). In addition, all four Notopterygium species had high levels of genetic differentiation, where the genetic variations in the cpDNA and ITS sequences mainly occurred among the populations within each species (Table 3 and  Supplementary Tables S5, S6). In general, N. franchetii has the most extensive natural geographic distribution range, but we found that its cpDNA and ITS sequences had low diversity. We consider that this low diversity may be due to harvesting and climatic changes, where many of the natural populations of N. franchetii have become extinct because of habitat destruction, thereby causing low diversity and high genetic differentiation (Lowe et al., 2005). N. incisum is another widely distributed species but we found that it had a high level of genetic diversity compared with other three species, which may be explained by the less extensive destruction of the wild populations of this species. According to the field investigations, we found that this species generally occurs in higher altitude areas (≥3000 m) compared with other Notopterygium species, and thus its less frequent harvesting might explain the high genetic variation (Wang G.N. et al., 2014). In addition, the high level of genetic diversity in N. oviforme according to this study might be explained by the lower altitude range of this species (1700-3200 m), which is consistent with a previous report of high species diversity at low altitudes (Lu et al., 2012). The lower genetic diversity of N. forrestii may be due to its narrow geographical distribution, where the smaller localized populations can interbreed and the gene flow is greater, thereby leading to a low level of diversity.
These Notopterygium species may also have been affected by adverse environmental changes in the high altitude QTP and adjacent areas. Thus, repeated climatic oscillations and geological events may have led to genetic drift and the fragmentation of habitats, thereby reducing their diversity (N. oviforme had slightly higher diversity compared with the other three species) and causing a high level of genetic differentiation among the populations of Notopterygium species (Hamrick and Loveless, 1989).

Relationships among Species
Phylogenetic analysis based on the cpDNA and ITS haplotypes showed that all four Notopterygium species formed a monophyletic clade with high bootstrap support (Figure 3 and Supplementary Figure S2). N. franchetii and N. oviforme shared a common branch in the phylogenetic tree based on the cpDNA sequences, and they also shared cpDNA haplotype H32. However, there were no shared ITS haplotypes among the four Notopterygium species, where each species formed an individual clade in the ITS tree. Thus, the ITS marker could identify the species at greater resolution than the cpDNA fragments. In general, cpDNA is a uniparentally inherited region whereas nuclear ITS fragments are biparentally inherited markers in most angiosperms, so ITS markers are superior for discriminating lineages and species than cpDNA fragments (Fan et al., 2011;Wang et al., 2012).
In addition, the shared cpDNA H32 haplotype was found in the parapatric populations (LF and YA) of N. franchetii and N. oviforme. These parapatric geographic distributions may have provided the opportunity for interspecific gene flow and hybridization among the two species. According to the field observations, we found that these two species have overlapping flowering times, which may have facilitated genetic introgression among these species. Previous studies have also suggested the occurrence of hybridization among species distributed in the same geographic regions and subsequent backcrosses with one of the parental species, where these processes resulted in high levels of shared plastid genotypes (Hamzeh et al., 2006;Wang Z. et al., 2014). However, it is also possible that incomplete lineage sorting could have lead to the sharing of cpDNA haplotypes among species. The perennial herb Notopterygium species have large population sizes and long generation times, which may have led to the sharing of ancestral polymorphisms among species.

Species Divergence and Population Dynamics History
Mountain barriers may play a key role in speciation and diversification because their topographic complexity can lead to ecological stratification and environmental heterogeneity (Fjeldså et al., 2012). In the present study, we estimated the divergence time of the four Notopterygium species based on three cpDNA fragments, which showed that their divergence occurred between about 3.6 Mya (95% HPD, 2.1-5.3 Mya) and 1.2 Mya (95% HPD, 0.67-1.8 Mya) in the Pliocene and Pleistocene periods. The divergence of N. forrestii and N. incisum was estimated as occurring between 2.24 Mya (95% HPD, 1.14-3.4 Mya) and 0.75 Mya (95% HPD, 0.4-1.14 Mya) in the early to middle Pleistocene period (Figure 4). We suggest that the divergence of the four Notopterygium species was significantly related to the uplift of the QTP. Previous studies and geological data indicated that the uplift of the QTP started in the Oligocene to Miocene (25-17 Mya), middle of the Miocene (15-13 Mya), late Miocene (8-7 Mya), or in the Pliocene to early Pleistocene period (3.6-1.8 Mya) (Harrison et al., 1992;Coleman and Hodges, 1995;Shi et al., 1998;Spicer et al., 2003). During the uplift of the QTP and adjacent Himalayan mountains, long-term geological events generated great environmental differences, which might have triggered the diversification of species in the genus Notopterygium. Other studies have also shown that the recent extensive uplift of the QTP and adjacent mountains triggered the lineage divergence and evolution of many herb species due to geographical isolation and climatic changes (Lei et al., 2007(Lei et al., , 2015. In addition, dramatic variations in the environment and climate might have affected the genetic structure and geographic distributions of the Notopterygium species. MaxEnt modeling showed that the two cold-tolerant species comprising N. incisum and N. franchetii exhibited significant range expansions from the LIG to the LGM period (Figures 5, 6). The mismatch analysis, neutrality test, and population growth index results also supported similar expansions by these two species. Therefore, we estimated the expansion times for N. incisum and N. franchetii as about 128-43 and 51-17 Kya, respectively, during the late Pleistocene (Table 5). We showed that the geographic ranges of these two species increased significantly during the ice ages in the Pleistocene. Demographic expansions of cold-tolerant tree species during the glacial periods have also been reported in high altitude areas of the QTP (Qiu et al., 2011;Wen et al., 2014). Moreover, repeated founder and bottleneck effects during the expansion processes may explain the low genetic variation in the two species. By contrast, we found that the populations of N. incisum and N. franchetii had high levels of genetic diversity in the southeast part of the QTP. For example, some N. incisum populations in Qinghai (population G), Sichuan (populations J, Q, U, and X), and Gansu (population Z) had high diversity and many more unique haplotypes. The N. franchetii populations in Sichuan (populations KG and YF) and Gansu (populations KO and KZ) also had high genetic diversity and a rich abundance of haplotypes ( Table 2). In addition, the LE population of N. oviforme and the LCA population of N. forrestii had high levels of haplotype diversity. These areas may have provided important glacial refugia for these endemic perennial herb species. Mountain areas at low latitudes can also provide relatively stable environmental conditions according to the "ecological stability hypothesis" (Qu et al., 2014;Lei et al., 2015), which implies that these populations should have high genetic diversity and a rich diversity of haplotypes (Tzedakis et al., 2002;Abbott and Brochmann, 2003;Petit et al., 2003;Schonswetter et al., 2005). Similar results have been obtained for other organisms, such as birds (Lei et al., 2007), mammals (Wu et al., 2005), spiders (Meng et al., 2008), and aphids (Huang et al., 2010).

Conservation Strategies for Endangered Notopterygium Species
The genus Notopterygium comprises unique perennial herbaceous plants with medicinal applications in China (Zhou et al., 2010). These species have high economic value so the market demand is great, especially for N. incisum and N. franchetii. However, in recent years, due to their continuous harvesting, slow growth rate, and low reproductive capacity, the natural populations of these Notopterygium species have been greatly depleted (Zhou et al., 2010). According to field investigations, we found that many of the previously recorded natural populations of species in the genus Notopterygium were extinct, and thus these important resources require urgent conservation and management.
According to the results of our population genetics analysis, we propose that the natural populations of wild Notopterygium species should be protected in situ, especially in the natural refugia areas (i.e., populations J, Q, U, and X of N. incisum; populations KO and KZ of N. franchetii; population LE of N. oviforme; and population LCA of N. forrestii). In addition, it is necessary to control all activities that deplete the sizes of the populations (e.g., illegal harvesting) and genetic fragmentations (e.g., habitat loss) (Hamilton et al., 2012).
In order to conserve the populations of a species, it is necessary to understand the genetic diversity and population structure of the natural populations (Schmitt, 2007). In this study, we found that the genetic variability and haplotype diversity were low for the widely distributed species, thereby indicating that the habitats have been destroyed or fragmented for these species, where interbreeding has occurred with nearby populations of individuals, thereby reducing the haplotype diversity. It is also necessary to protect the populations in different regions in order to increase the genetic links among populations. Finally, the mature seeds from each population should be collected and artificially planted with other populations in order to improve the habitats and to strengthen the gene exchange among populations (Cabrera-Toledo et al., 2012).

AUTHOR CONTRIBUTIONS
Z-HL designed and conceived the study. KS and YJ performed the experiments. Z-HL, KS, YJ, F-LC, and UZ contributed materials/analysis tools. KS and Z-HL wrote the manuscript. KS, YJ, and Z-HL revised the manuscript. All of the authors finally approved the manuscript.

ACKNOWLEDGMENTS
This study was supported by the National Natural Science Foundation of China (31470400) and by the Program for Changjiang Scholars and Innovative Research Team in University (PCSIRT, No. IRT_15R55).