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 high-altitude 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.

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₂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 sequences were deposited in GenBank under accession numbers MF787385–MF787579.

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).

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 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 log-normal 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.

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.

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¹ 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).

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.

In total, 55 cpDNA haplotypes were detected in the four Notopterygium species (Figures 1, 2). Most of the haplotypes were species-specific, except the H32 haplotype was shared by N. franchetii and N. oviforme. N. incisum contained 31 haplotypes, where haplotypes H1–H7, H12, H18, and H26 were shared among populations, and the remainder were unique to each population. Populations from the southeast part of the QTP (G, J, Q, U, Z, HB, and HC; see Table 1 for the site codes) had the highest haplotype diversity. N. franchetii had nine haplotypes with nine mutation sites. Populations from the west part of China (KG, KO, KZ, and YF) also had the highest haplotype diversity. N. oviforme had 12 haplotypes with 16 mutation sites, where the LE population had the highest haplotype diversity for this species. N. forrestii had three haplotypes with two mutation sites, where H55 was unique to the LCB population, and haplotypes H53 and H54 were shared by the other populations.

Notopterygium oviforme had the highest levels of genetic diversity and nucleotide 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).

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).

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). In the four individual species, the genetic variations among populations (N. incisum, F_ST = 0.8196; N. franchetii, F_ST = 0.8391; N. oviforme, F_ST = 0.8474; and N. forrestii, F_ST = 0.7585) were significantly higher than those within populations (Supplementary 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 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).

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).

We estimated the divergence times between the four species of Notopterygium based on a range of mutation rates (1.0–3.0 × 10^-9 s/s/y). The first divergence among the four species occurred between approximately 3.6 Mya (95% highest posterior density (HPD), 2.1–5.3 Mya) and 1.2 Mya (95% HPD, 0.67–1.8 Mya), whereas the estimated divergence between N. forrestii and N. incisum occurred between 2.24 Mya (95% HPD, 1.14–3.4 Mya) and 0.75 Mya (95% HPD, 0.4–1.14 Mya). In addition, the divergence between the major lineages of N. franchetii and N. oviforme occurred between 1.3 Mya (95% HPD, 0.54–2.2 Mya) and 0.42 Mya (95% HPD, 0.18–0.73 Mya) (Figure 4).

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.

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).

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; Li et al., 2013; 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.

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, 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).

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).