Parnassius (Lepidoptera: Papilionidae) is a genus of attractive butterflies mainly distributed in the mountainous areas of Central Asia, the Himalayas, and western China. In this study, we used the internal transcribed spacer (ITS1 and ITS2) sequence data as DNA barcodes to characterize the genetic differentiation and conduct the phylogenetic analysis and divergence time estimation of the 17 Parnassius species collected in China. Species identification and genetic differentiation analysis suggest that the ITS barcode is an effective marker for Parnassius species identification; additionally, a relatively high level of genetic diversity and low level of gene flow were detected in the five Parnassius species with diverse geographic populations. Phylogenetic analysis indicates that the 17 species studied were clustered in six clades (subgenera), with subgenus Parnassius at the basal position in the phylogenetic trees. Bayesian divergence time estimation shows that the genus originated about 18 million years ago during the early Miocene, correlated with orogenic events in the distribution region, probably southwestern China about 20–10 million years ago. Our estimated phylochronology also suggests that the Parnassius interspecific and intraspecific divergences were probably related with the rapid rising of the Qinghai-Tibet Plateau, the Tibet Movement, the Kunlun-Yellow River Tectonic Movement, and global cooling associated with intensified glaciation in the region during the Quaternary Period.

The Parnassius are one of the most charming butterfly groups of subfamily Parnassiinae within family Papilionidae, often called “Apollos.” Species in the genus, totaling about 60 worldwide, are distributed mostly in high-altitude mountainous areas of the Himalayas, Central Asia, western China, and other parts of northern Eurasia (Condamine 2018, Katoh et al. 2005, Omoto et al. 2009), with 33 species reported in China (Chou 1999). Several recent attempts to resolve the phylogenetic relationships within Parnassius based on molecular and morphological data sets have provided considerable resolutions, but some intersubgeneric relationships have not been well supported while others were mutually contradicted among studies (Condamine et al. 2013, 2018b; Katoh et al. 2005; Michel et al. 2008; Omoto et al. 2004; Zheng et al. 2018). For example, Omoto et al. (2004) used the analysis of mitochondrial ND5 gene to propose the division of eight subgenera. Michel et al. (2008) obtained similar results based on four mitochondrial DNA segments (COI, ND1, ND5, 16S rRNA), and named these eight subgenera (Parnassius, Kailasius, Koramius, Kreizbergia, Lingamius, Driopa, Sachaia, and Tadumia); however, Katoh et al. (2005) and Zheng et al. (2018) acquired different results about the subgeneric relations, based on mitochondrial genes. Additionally, contradictory results were reported by Condamine et al. (2013) and Condamine et al. (2018b) using similar criteria.

The elevation of the Qinghai-Tibet Plateau experienced dramatic uplift during the frequent orogeny of the Cenozoic epoch, and its southern part likely reached an elevation comparable to present-day elevation during the Miocene, transforming the local environments, as well as causing the aridification of Central Asia (Du et al. 2019, Favre et al. 2015, Najman et al. 2010), and its high-elevation environments also experienced dramatic physical changes during the cyclical expansion and retreat of glacial ice sheets in the Quaternary (Shi 2002, Xu et al. 2010, Zheng et al. 2002). These consequences of such extraordinary geological and environmental changes for animal and plant life in alpine ecosystems included large geographic range shifts, long periods of isolation, and, in some cases, recent, rapid diversifications (Schoville and Roderick 2009). Because the species abundance of Parnassius is much higher than those of other genera in the subfamily Parnassiinae, these “Apollo” butterflies may have experienced a unique evolutionary process in recent geological history (Omoto et al. 2004, Rebourg et al. 2006). Previous studies attempted to associate the Parnassius diversification with geological events (such as the rise of the Qinghai-Tibet Plateau and Himalayas, and related orogenies during Cenozoic) and climate changes (such as the Quaternary glaciations) (Condamine et al. 2018b, Favre et al. 2015, Lei et al. 2014, McLean et al. 2018, Omoto et al. 2009). However, it was suggested that current biodiversity has been caused by the joint effect of multiple factors (biogeography, species traits, environmental drivers, and species extinction) rather than a single factor, and that addressing the origin and evolution of the Parnassius requires a reliable and accurate phylogenetic framework for pinpointing significant phylogenetic events (Condamine et al. 2018a, 2018b). Although a lot of studies have investigated their phylogeny, the phylogenetic backbone of Parnassius has not yet been resolved, and many evolutionary hypotheses were only founded on poorly supported phylogenetic reconstructions (Condamine et al. 2013, 2018b; Nazari et al. 2007; Omoto et al. 2009).

For species identification, the recently developed DNA barcoding technique has been proven useful to corroborate the traditional morphological approaches (Hebert et al. 2003). The method assumes that the genetic variation between two species exceeds that within the species for selected DNA segment (Badotti et al. 2018). COI gene is commonly used as the standard barcode for animals and rbcL+matK for plants; the nuclear internal transcribed spacer (ITS) sequences is also recommended as a candidate marker for plant, animal, and fungus species identifications (Dentinger et al. 2011, Xu 2016, Yao et al. 2010, Zhu et al. 2017). Although the ITS in eukaryotes contains two separate regions (ITS1 and ITS2) (Hillis and Dixon 1991), recent studies suggest concerted evolution among them and the ITS can be treated as a single gene (Rampersad 2014). Furthermore, it is shown that in both yeasts and fruit flies, ITS is essential for the formation of ribosomal subunits and the evolutionary rate is relatively fast (Morgan and Blair 1998, Rampersad 2014). Therefore, ITS has been frequently utilized as a marker for phylogenetic analyses at the generic and specific levels (Brown et al. 2000, Coleman 2003, Poczai and Hyvönen 2010).

In the present study, we sequenced ITS1 and ITS2 from 267 individuals of 17 Parnassius species collected from Qinghai-Tibet Plateau and its neighboring areas. We tested the applicability of ITS in species identification of Parnassius by analyzing their inter- and intraspecific genetic variation. We used ITS1 and ITS2, for the first time, to reconstruct their phylogeny with multiple methods, and to estimate their divergence times using multiple calibrations with relaxed molecular clock methods, in order to establish a reference framework to evaluate subgenus-level relationships within Parnassius and assess the possible geological and palaeoenvironmental events that are likely to be the driving forces for the divergences of Parnassius in the region.

Sampling and DNA sequencing. We collected 267 adults of 17 Parnassius species from various locations in China (Fig. 1; Table 1). Species identification was initially based on morphological traits. Samples used in this study were preserved in the Laboratory of Molecular Evolution and Biodiversity, College of Life Sciences, Anhui Normal University, Wuhu, and Molecular Paleobiology Lab at Nanjing Institute of Geology and Palaeontology, Chinese Academy of Sciences, Nanjing. Genomic DNA was extracted from the leg or thorax tissues using Rapid Animal Genomic DNA Isolation Kit (Sangon Biotech. Co. Ltd, Shanghai, China). The DNA was stored in Tris-EDTA buffer at –20°C.

Fig. 1

Sampling localities of the 17 Parnassius species used in this study; the code information in the figure refers to Table 1.

Fig. 1

Sampling localities of the 17 Parnassius species used in this study; the code information in the figure refers to Table 1.

Close modal
Table 1

List of sample localities of the Parnassius species and two other Parnassiinae species used in this study.

List of sample localities of the Parnassius species and two other Parnassiinae species used in this study.
List of sample localities of the Parnassius species and two other Parnassiinae species used in this study.
Table 1

Continued.

Continued.
Continued.
Table 1

Continued.

Continued.
Continued.
Table 1

Continued.

Continued.
Continued.
Table 1

Continued.

Continued.
Continued.
Table 1

Continued.

Continued.
Continued.
Table 1

Continued.

Continued.
Continued.

The amplification of the ITS1 regions was conducted using polymerase chain reaction (PCR) with the forward primer 18sF1 (5′-TACACACCGCCCGTCGCTAC TA-3′) and reverse primer 5.8sB1d (5′-ATGTGCGTTCRAAATGTCGATGTTCA-3′). ITS2 regions were amplified using the forward primer 5.8sFc (5′-TGAACATCGA CATTTYGAACGCACAT-3′) and reverse primer 28sB1d (5′-TTCTTTTCCTCC SCTTAYTRATATGCTTAA-3′) (Ji et al. 2003), in 50 µl reagents containing 6.0 µl 10× PCR buffer, 8.0 µl MgCl2, 1.5 µl dNTPs, 2.0 µl each primer, 1.5 µl DNA template, 1.0 µl Taq DNA polymerase (1.0 U), and 28 µl ddH2O. The thermal cycle parameters were: an initial denaturation at 95°C for 5 min; followed by 35 cycles: denaturation at 95°C for 50 s, annealing at 59.5°C (ITS1) and 54°C (ITS2) for 1 min, and extension at 72°C for 2 min; and a final extension at 72°C for 10 min. PCR products were purified using a DNA purification kit (Sangon Biotech. Co. Ltd, Shanghai, China) and sequenced from both directions by General Biosystems Co. Ltd, Anhui, China.

Data analysis. For all ITS1 and ITS2 sequences obtained in this study, multiple sequence alignment was conducted by software MAFFT Version 7.1 using the FFT-NS-i algorithm (Katoh and Standley 2013). The nucleotide composition of the sequences was calculated using MEGA version 7.0 (Kumar et al. 2016). The Kimura 2-parameter (K2P) inter- and intraspecific genetic distances were estimated using MEGA version 7.0. The haplotype diversity (Hd), nucleotide diversity (Pi), genetic differentiation index (Fst), and gene flow (Nm) (Grant and Bowen 1998, Slatkin and Maddison 1989, Wright 1965) were analyzed using DnaSP version 6.0 (Rozas et al. 2017). Analysis of molecular variance (AMOVA) was performed by using Arlequin version 3.1 (Excoffier et al. 2005) with 1,000 permutations to compare levels of genetic diversity within and among populations.

Species identification analysis. Six distance parameters were calculated using MEGA version 7.0 for inter- and intraspecific variation, including three parameters for interspecific divergences: (1) average interspecific distance between all species; (2) average theta prime (the mean pairwise distance between all samples); (3) minimum interspecific distance; and three parameters for intraspecific variation: (1) average intraspecific distance among all samples within each species; (2) theta (the mean pairwise distance within each species with at least two representatives); and (3) average coalescent depth (the maximum intraspecific distance within each species with at least two individuals) (Lahaye et al. 2008, Meier et al. 2008, Meyer and Paulay 2005). Distribution of the pairwise inter- and intraspecific distances for ITS1 and ITS2 was calculated with the K2P method using the software TaxonDNA version 1.8 (Meier et al. 2006). The barcoding gaps were graphed by the distribution of pairwise inter- and intraspecific distances of the two DNA markers. Species identification ability was evaluated using the “Best match,” “Best close match,” and “All species barcodes” strategies in TaxonDNA, based on the K2P method (Meier et al. 2006). For the “Best match,” a query is assigned the species name of its best-matching barcode sequences, regardless of how similar the query and barcode sequences are; while for the “Best close match,” a threshold similarity value is required to define how similar a barcode match needs to be before it can be identified. The “All species barcodes” is the most rigorous application for identifying queries that are assigned a species name only if the query is followed by all known barcodes for a particular species and only if there are at least two conspecific matches (Meier et al. 2006).

Phylogenetic analysis.Parnassius phylogeny was reconstructed with the Bayesian inference (BI) and maximum likelihood (ML) methods based on the concatenated ITS1 and ITS2 sequence data using Sericinus montelus Gray (MN129452 and MN129720) and Luehdorfia chinensis Leech (AB071924.1) as the outgroups. Bayesian analysis was performed in MrBayes version 3.2 (Ronquist et al. 2012) under the general time reversible GTR+G nucleotide substitution models determined by PartitionFinder version 1.1.1 (Lanfear et al. 2012). Two independent Markov chain Monte Carlo (MCMC) runs were allowed to go for 4 million generations with sampling each 200th generation. Each run had four chains, one cold and three heated. Convergence of the Bayesian runs was ensured by checking the average standard deviation of split frequencies (StdDev), and the potential scale reduction factor (PSRF) values, and by examining the effective sample size (ESS) of all parameters in Tracer version 1.7.1 (Rambaut et al. 2018). To reach a good convergence, the standard deviation value should be below 0.01, the PSRF value close to 1.00, and the ESS value larger than 200. The first 25% of the sampled generations were discarded as burn-in samples. The resultant posterior probability (BPP) was obtained as the supporting values of each tree node. ML analysis was conducted with software IQ-TREE version 1.6.8 under the GTR+F+R3 models ascertained by the ModelFinder (Kalyaanamoorthy et al. 2017, Nguyen et al. 2014). In the ML analysis, 5,000 ultrafast bootstraps were performed to obtain the ultrafast bootstrap support values (BS) of each node.

Divergence time estimation. The calibrations for the molecular dating were based on two butterfly fossils: Praepapilio colorado Durden (Papilionidae) from the Green River Shale of Colorado (USA) of mid-Eocene early Lutetian age (41.2–47.8 Ma) and Thaites ruminiana Scudder (subfamily Parnassiinae) from Aix-en-Provence (southern France) of late Oligocene Chattian age (23.03–28.1 Ma) (Condamine et al. 2018b, Durden and Rose 1978, Jong 2017); thus, the crown group divergence of the Parnassiinae was constrained between 23.03 Ma and 47.8 Ma. The earliest divergence time of five Parnassius subgenera exclusive of subgenus Parnassius was set to be 16–37 Ma according to the divergence time estimates of their host plant Corydalis (subfamily Fumarioideae) and the insect–host plant coevolutionary scenario (Pérez-Gutiérrez et al. 2015, Su et al. 2017).

Our phylochronological analysis here was conducted using BEAST version 1.8.3 (Drummond et al. 2012) under a relaxed clock model with an uncorrelated lognormal distribution (Drummond et al. 2006). Lognormal priors with soft bounds and 95% confidence intervals (CI) were used for each fossil constraint. The tree prior was set to the birth–death process with incomplete sampling (Stadler 2009), and the nucleotide substitution model was set to the GTR+G model; the MCMC chains were run for 60 million generations with sampling every 1,000 generations; convergence was assessed by Tracer version 1.7.1 (Rambaut et al. 2018), determining whether the ESS of all parameters was larger than 200 as recommended; the nodal heights and maximum credibility tree were generated with TreeAnnotator version 1.8.3 (Rambaut and Drummond 2016), with the first 6,000 trees being discarded as burn-ins. Finally, the maximum credibility tree was obtained by Figtree version 1.4.3 (Rambaut 2009).

DNA sequence data analysis. The ITS1 and ITS2 regions of 267 individuals from 17 Parnassius species were successfully sequenced, and annotated and defined by the NCBI database, submitted to and deposited into GenBank (accession No. MN129185–MN129451 and MN129453–MN129719). Our analyses show substantial interspecific variations with the longest (Parnassius acdestis Grum-Grshimailo) and shortest (P. hide Koiwaya) of ITS1 sequences being 680 bp and 557 bp, and the longest (P. acdestis) and shortest (P. glacialis Butler) of ITS2 sequences being 657 bp and 485 bp in length, respectively (Table 2). On the other hand, intraspecific variations are relatively low, but in a few cases, are significant; for example, the size variation of P. acdestis ITS1 reached up to 63 bp. Such variations may be attributed to replication slippage, unequal crossing over, and biased gene conversion, as well as geographical isolations and other factors (Hlinka et al. 2002, Platas et al. 2004).

Table 2

Characteristics of the ITS1 and ITS2 sequences of 17 Parnassius species determined in this study.

Characteristics of the ITS1 and ITS2 sequences of 17 Parnassius species determined in this study.
Characteristics of the ITS1 and ITS2 sequences of 17 Parnassius species determined in this study.

Average AT content of the ITS2 regions is lower than the ITS1 regions (52.8% versus 57.7%). For both ITS1 and ITS2 sequences, mean AT contents vary remarkably: the highest mean AT content of ITS1 being 61.5% (P. glacialis) and the lowest 54.9% (P. nomion Fischer von Waldheim), while the highest and the lowest of ITS2 are 55.5% (P. hide) and 50.2% (P. imperator Oberthür), respectively. The total length of the aligned ITS1 sequences are 792 bp long, with 344 variable sites and 284 parsimony informative sites; while the aligned ITS2 sequences are 836 bp long, with 315 variable sites and 262 parsimony informative sites.

Species identification. The resulting six distance parameters are shown in Table 3. The results revealed that both ITS1 and ITS2 exhibited a relatively higher interspecific and lower intraspecific divergences. For example, the average interspecific and intraspecific distances were 0.2163 and 0.0025 (ITS1), and 0.1850 and 0.0024 (ITS2), respectively. Although the minimum interspecific distance (0.0023) of the ITS2 sequences was somewhat less than the maximum intraspecific distance (0.0055), the minimum interspecific distance (0.0062) of the ITS1 sequences was greater than the maximum intraspecific distance (0.0049). The barcoding gaps between inter- and intraspecific distances are shown in Fig. 2. There was overlap between inter- and intraspecific distances in both ITS1 and ITS2 regions. The species identification efficiency of the ITS1 and ITS2 regions were tested, as shown in Table 4. The results indicated that under the “Best match” and “Best close match” criteria, the two regions had the same species identification rate of 98.12%; however, both of the region's species identification rates were 97.37% and 83.89%, respectively, using the “All species barcodes” standard.

Table 3

Inter- and intraspecific variations of ITS1 and ITS2 sequences of all individuals of 17 Parnassius species.

Inter- and intraspecific variations of ITS1 and ITS2 sequences of all individuals of 17 Parnassius species.
Inter- and intraspecific variations of ITS1 and ITS2 sequences of all individuals of 17 Parnassius species.
Fig. 2

Relative distribution of the inter- and intraspecific distances for ITS1 (a) and ITS2 (b) regions.

Fig. 2

Relative distribution of the inter- and intraspecific distances for ITS1 (a) and ITS2 (b) regions.

Close modal
Table 4

Species identification efficiency of ITS1 and ITS2 regions in Parnassius.

Species identification efficiency of ITS1 and ITS2 regions in Parnassius.
Species identification efficiency of ITS1 and ITS2 regions in Parnassius.

Genetic differentiation and phylogenetic analysis. The BI and ML trees (Fig. 3) from the ITS1+ ITS2 data sets have the same topology, with slight differences of node supporting values, and all the node supporting values were relatively strong (BS > 86%, BPP > 0.85; except node A: BS = 68%, BPP = 0.57). It is shown that the 17 Parnassius species are divided into six clades (subgenera) with their relationship as follows: (([Driopa + Kreizbergia] + [Kailasius + Tadumia] + Koramius)) + Parnassius). The five Parnassius species (P. stubbendorfii Ménétriés, P. nomion, P. simo Gray, P. orleans Oberthür, P. acdestis) with different populations all contain two major geographic branches (Fig. 3, red branches).

Fig. 3

Bayesian inference (BI) and maximum likelihood (ML) trees inferred from 17 Parnassius species based on ITS1 and ITS2. Samples of different colors represent six different subgenera; red and black bars represent different populations within species. Asterisk (*): ML ultrafast bootstrap support values and BI posterior probability values of 100 and 1, respectively.

Fig. 3

Bayesian inference (BI) and maximum likelihood (ML) trees inferred from 17 Parnassius species based on ITS1 and ITS2. Samples of different colors represent six different subgenera; red and black bars represent different populations within species. Asterisk (*): ML ultrafast bootstrap support values and BI posterior probability values of 100 and 1, respectively.

Close modal

The interspecific genetic distances among 17 Parnassius species and the genetic distance within and among populations of the five Parnassius species are shown in Tables 5 and 6, respectively. The results revealed that the interspecific genetic distances ranged from 0.006 (between P. cephalus Grum-Grshimailo and P. choui Huang and Shi) to 0.310 (between P. actius Eversmann and P. acco Gray), and the genetic distance within populations of the five Parnassius species were significantly smaller than those among populations. The haplotype diversity (Hd) of the total population of the five Parnassius species ranged from 0.442 to 0.909, and the nucleotide diversity (Pi) ranged from 0.00058 to 0.00494 (Table 6). AMOVA showed that the genetic variation among populations of the five Parnassius species was significantly greater than those within population (P ≤ 0.003) (Table 7). The total population genetic differentiation index (Fst) and gene flow (Nm) of the five Parnassius species ranged from 0.6667 to 1.0000 and 0 to 0.13 (Table 6), respectively.

Table 5

Kimura 2-Parameter genetic distance among 17 Parnassius species in this study.

Kimura 2-Parameter genetic distance among 17 Parnassius species in this study.
Kimura 2-Parameter genetic distance among 17 Parnassius species in this study.
Table 6

Analysis of intraspecific genetic differentiation* of the five Parnassius species.

Analysis of intraspecific genetic differentiation* of the five Parnassius species.
Analysis of intraspecific genetic differentiation* of the five Parnassius species.
Table 7

Analysis of molecular variance of genetic variance among different populations of the five Parnassius species.

Analysis of molecular variance of genetic variance among different populations of the five Parnassius species.
Analysis of molecular variance of genetic variance among different populations of the five Parnassius species.

Divergence time analysis. As shown in our relaxed molecular clock analysis using C1 (divergence of the subfamily Parnassiinae) and C2 (divergence of the genus Parnassius exclusive of subgenus Parnassius) as the calibration points (Fig. 4), the diversification time of genus Parnassius (i.e., the splitting of subgenus Parnassius and other subgenera) started about 18.41 Ma (95% CI: 23.83–14.13 Ma) in early Miocene. The divergence of two major clades (subgenera Driopa + Kreizbergia versus Tadumia + Kailasius + Koramius) within genus Parnassius occurred about 17.39 Ma (95% CI: 22.11–13.49 Ma) during early Miocene to mid-Miocene, followed by the splitting of subgenus Koramius from the clade of Kailasius + Tadumia about 14.41 Ma (95% CI: 19.02–10.45 Ma), and Kailasius from Tadumia about 12.15 Ma (95% CI: 16.34–8.48 Ma). Subgenus Kreizbergia diverged from Driopa about 12.07 Ma (95% CI: 16.78–8.07 Ma). According to our divergence data analysis, the interspecific divergences of the 17 Parnassius species and intraspecific differentiations of the five Parnassius species mostly concentrated from about 3.5 Ma and earlier during late Pliocene and Quaternary time, when the global temperature continually declined (Fig. 4, lower graph) associated with widespread glaciation.

Fig. 4

Phylochronology of Parnassius species based on ITS1 and ITS2 data under a Bayesian relaxed clock analysis (BEAST, version 1.8.3). C1, C2: calibration points (see text section: Divergence time estimation). Numbers and numbers in parentheses at nodes: median time estimates and 95% confidence interval (CI) in million years ago (Ma), respectively. The red and black bars: different populations within species. Lower graph: The main global temperature curve (Zachos et al. 2001) from the late Oligocene to the Pleistocene with palaeoclimatological information (monsoon and related Central Asian aridification) (from Favre et al. 2015).

Fig. 4

Phylochronology of Parnassius species based on ITS1 and ITS2 data under a Bayesian relaxed clock analysis (BEAST, version 1.8.3). C1, C2: calibration points (see text section: Divergence time estimation). Numbers and numbers in parentheses at nodes: median time estimates and 95% confidence interval (CI) in million years ago (Ma), respectively. The red and black bars: different populations within species. Lower graph: The main global temperature curve (Zachos et al. 2001) from the late Oligocene to the Pleistocene with palaeoclimatological information (monsoon and related Central Asian aridification) (from Favre et al. 2015).

Close modal

Species identification. Previous studies have shown that ideal candidate DNA barcodes should, first, have sufficient variations to discriminate among species and also need to be sufficiently conserved so that there is less variability within species than between species. Second, ideal DNA barcodes should harbor a “barcode gap,” where the distribution of intraspecific variation and interspecific divergence have discrete distributions and no overlap. Third, they should have a high success rate of species identification (Hartvig et al. 2015, Hebert et al. 2003, Meier et al. 2008, Meyer and Paulay 2005, Yao et al. 2010). In the present study, our ITS1 and ITS2 regions presented the most promising universal DNA barcodes for authenticating Parnassius species as assessed by several criteria. First, determination of genetic divergence using six distance parameters confirmed that the ITS1 and ITS2 regions possessed high interspecific divergences and low intraspecific variations (Table 3). Second, DNA barcoding gaps analyses indicated that there existed slight overlaps between inter- and intraspecific distances for the two evaluated regions (Fig. 2). Third, species identification efficiency via three common criteria (“Best match,” “Best close match,” and “All species barcodes”) as suggested by Meier et al. (2006) indicated that both ITS1 and ITS2 harbored significantly high species identification rates (98.12%) and were able to correctly identify 261 out of 267 individuals using the “Best match” and “Best close match” criteria; however, under the “All species barcodes” standard, the correct identification rate of ITS2 region (83.89%) was slightly low and ITS1 higher (97.37%) (Table 4). Overall, our newly determined ITS1 and ITS2 markers are useful and reasonably efficient to Parnassius species identification.

Genetic differentiation and phylogenetic analysis. In this study, the basal position of subgenus Parnassius is congruent with results obtained by Omoto et al. (2004), Rebourg et al. (2006), Michel et al. (2008), and Condamine et al. (2018b). However, the relationships among other subgenera resolved in this study differ from previous studies, prompting further studies with more comprehensive data sets to resolve.

The species relationships within each subgenus obtained in this study based on ITS data essentially conform to the traditional morphological classifications, as also shown in Omoto et al. (2004), Katoh et al. (2005), Michel et al. (2008), and Zheng et al. (2018), based on mitochondrial or mitochondrial plus nuclear gene sequence data. Our analyses also suggest that P. acdestis belongs to subgenus Kailasius instead of subgenus Koramius, despite the similar morphology in wing venation and sphragis, verifying the analysis by Omoto et al. (2009) based on mitochondrial data.

The haplotype diversity (Hd) and the nucleotide diversity (Pi) analysis of the total population of the five Parnassius species showed that they all harbored a relatively low level of nucleotide diversity (Pi < 0.005) and a high level of haplotype diversity (Hd > 0.5) except for P. nomion (Table 6). AMOVA indicated that their genetic variations came mainly from different geographical populations (Table 7), and the results correspond to the those of the genetic distance. The total population genetic differentiation index (Fst) and gene flow (Nm) analysis suggested that their genetic differentiation levels were significantly high (Fst > 0.25) with relatively limited gene flows (Nm < 1) (Table 6).

Divergence time analysis. Our phylochronological results (Fig. 4) generally agree with those of Condamine et al. (2013), who determined the earliest Parnassius divergence at about 17 Ma (95% CI: 22–13 Ma), but significantly differed from the results obtained by Nazari et al. (2007), Omoto et al. (2009), and Condamine et al. (2018b), who provided dates of the same event at about 39–34 Ma (late Eocene), 24.3 Ma (late Oligocene), and 13.4 Ma (95% CI: 16.6–10.5 Ma) (middle Miocene), respectively. We found that the major source of differences in dating the tree came from calibration points and their distribution in the tree. In this analysis, we adopted, in our calibration system (Fig. 4), high-quality fossil dates of Papilionidae and Parnassiinae and relevant dates from the host plants (Corydalis) based on the coevolution perspective, which are likely better proxies approaching the true reference time frame.

The uplift of the Qinghai-Tibetan Plateau and Himalayas caused a dramatic climatic and ecological shift. The forests were replaced by grasslands, and the climate gradually became drier, colder, and windier; glaciers started to develop, deserts formed (Wu et al. 2001); and accordingly, an area of worldwide importance for biodiversity in the Qinghai-Tibetan Plateau gradually developed, due to the unique geomorphologic configuration, complex land and climate conditions, as well as the distinct geological history that gave rise to the endemic, specialized montane animal and plant species (Lei et al. 2014).

In this study, the estimated divergences of Parnassius coincide and are likely related with geological events in the distribution regions of the butterflies, especially the progressive uplifting of the Himalayas and the Qinghai-Tibet Plateau during the early to middle Miocene (20–10 Ma) (Favre et al. 2015, Lu and Guo 2014, Molnar and Stock 2009). Studies show that these geological events remarkably changed the atmospheric circulations that gave rise to the intensified Asian monsoon and Central Asian aridification (Guo et al. 2008, Miao et al. 2012, Sun and Wang 2005, Wan et al. 2007). Afterwards, the two profound strengthening of these two events were accompanied by large-scale cooling of the Asian continent, which occurred at about 8 Ma and 3 Ma in the late Miocene and late Pliocene, respectively; and the latter event (3 Ma) may have been related to the Parnassius interspecific divergences for the extant species in our samples (Guo et al. 2011, Herbert et al. 2016, Wan et al. 2007). Meanwhile, the Tibet Movement between 3.6 and 1.7 Ma may also be responsible for interspecific genetic differentiation and speciation in the genus (Li et al. 1996).

The climatic oscillations in the Quaternary had profound effects on the organisms now inhabiting alpine ecosystems (Hewitt 2000). Numerous studies have documented effects of climatic cycles on population genetic diversifications, as well as both recent and ancient speciation events in alpine organisms (Chiocchio et al. 2017, DeChaine and Martin 2006, Sandel et al. 2011, Schoville and Roderick 2009). For example, the systematics and biogeography studies of Parnassius phoebus complexes in North America showed that P. smintheus Doubleday and P. behrii Edwards differentiated in the Pleistocene due to alpine glaciers (Schoville and Roderick 2009). The Qinghai-Tibetan Plateau region was more strongly affected by the widespread Quaternary glaciations than other regions of the world (Lei et al. 2014, Owen and Dortch 2014, Yang et al. 2008, Zhou et al. 2006), and the related geological events such as the Kunlun-Yellow River Tectonic Movement (1.1–0.6 Ma) (Wu et al. 2001, Zhou et al. 2006) are likely responsible for the habitat fragmentation and isolation of the Parnassius species and intraspecific differentiations of this study, judging from their space and time agreement between the butterfly phylogenesis and the paleoenvironmental events in the region.

This work was supported by the National Natural Science Foundation of China (41972029, 41472028), and the funds from the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB26000000).

Badotti,
F.,
P.L.C.
Fonseca
,
L.M.R.
Tomé
,
D.T.
Nunes
and
A.
Góes-Neto
.
2018
.
ITS and secondary biomarkers in fungi: Review on the evolution of their use based on scientific publications.
Braz. J. Bot.
41
:
471
479
.
Brown,
B.,
R.
Emberson
and
A.
Paterson
.
2000
.
Phylogenetic relationships within the genus Wiseana (Lepidoptera: Hepialidae).
N. Z. J. Zool.
27
:
1
14
.
Chiocchio,
A.,
R.
Bisconti
,
M.
Zampiglia
,
G.
Nascetti
and
D.
Canestrelli
.
2017
.
Quaternary history, population genetic structure and diversity of the cold-adapted Alpine newt Ichthyosaura alpestris in peninsular Italy.
Sci. Rep.
7
:
2955
.
Chou,
I.
1999
.
Monographia Rhopalocerorum Sinensium.
Henan Scientific and Technological Publishing House
,
Zhenzhou, China
.
191
pp.
Coleman,
A.W.
2003
.
ITS2 is a double-edged tool for eukaryote evolutionary comparisons.
Trends Genet.
19
:
370
375
.
Condamine,
F.L.
2018
.
Limited by the roof of the world: Mountain radiations of Apollo swallowtails controlled by diversity-dependence processes.
Biol. Lett.
14
:
20170622
.
Condamine,
F.L.,
B.
Nabholz
,
A.L.
Clamens
,
J.R.
Dupuis
and
F.A.H.
Sperling
.
2018a
.
Mitochondrial phylogenomics, the origin of swallowtail butterflies, and the impact of the number of clocks in Bayesian molecular dating.
Syst. Entomol.
43
:
460
480
.
Condamine,
F.L.,
J.
Rolland
,
S.
Höhna
,
F.A.H.
Sperling
and
I.
Sanmartín
.
2018b
.
Testing the role of the Red Queen and Court Jester as drivers of the macroevolution of Apollo butterflies.
Syst. Biol.
67
:
940
964
.
Condamine,
F.L.,
F.A.H.
Sperling
and
G.J.
Kergoat
.
2013
.
Global biogeographical pattern of swallowtail diversification demonstrates alternative colonization routes in the Northern and Southern hemispheres.
J. Biogeogr.
40
:
9
23
.
DeChaine,
E.G.
and
A.P.
Martin
.
2006
.
Using coalescent simulations to test the impact of Quaternary climate cycles on divergence in an alpine plant–insect association.
Evolution
60
:
1004
1013
.
Dentinger,
B.T.M.,
M.Y.
Didukh
and
J.M.
Moncalvo
.
2011
.
Comparing COI and ITS as DNA barcode markers for mushrooms and allies (Agaricomycotina).
PLoS ONE
6
:
e25081
.
Drummond,
A.J.,
S.Y.W.
Ho
,
M.J.
Phillips
and
A.
Rambaut
.
2006
.
Relaxed phylogenetics and dating with confidence.
PLoS Biol.
4
:
e88
.
Drummond,
A.J.,
M.A.
Suchard
,
D.
Xie
and
A.
Rambaut
.
2012
.
Bayesian phylogenetics with BEAUti and the BEAST 1.7.
Mol. Biol. Evol.
29
:
1969
1973
.
Du,
K.,
L.R.
Yang
,
R.
Zhang
,
L.P.
Yue
,
H.J.
Guo
,
Y.X.
Zhang
,
J.X.
Li
and
H.J.
Gong
.
2019
.
Cenozoic tectonics and landform evolution in the Qilian Mountains and adjacent areas.
Int. Geol. Rev.
doi: 10.1080/00206814.2019.1627588.
Durden,
C.J.
and
H.
Rose
.
1978
.
Butterflies from the middle Eocene: The earliest occurrence of fossil Papilionoidea (Lepidoptera).
Pearce-Sellards Ser. Tex. Mem. Mus.
29
:
1
25
.
Excoffier,
L.,
G.
Laval
and
S.
Schneider
.
2005
.
Arlequin (version 3.0): An integrated software package for population genetics data analysis.
Evol. Bioinform. Online
1
:
47
50
.
Favre,
A.,
M.
Päckert
,
S.U.
Pauls
,
S.C.
Jähnig
,
D.
Uhl
,
I.
Michalak
and
A.N.
Muellner-Riehl
.
2015
.
The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas.
Biol. Rev.
90
:
236
253
.
Grant,
W.S.
and
B.W.
Bowen
.
1998
.
Shallow population histories in deep evolutionary lineages of marine fishes: Insights from sardines and anchovies and lessons for conservation.
J. Hered.
89
:
415
426
.
Guo,
X.G.,
X.
Dai
,
D.
Chen
,
T.J.
Papenfuss
,
N.B.
Ananjeva
,
D.A.
Melnikov
and
Y.Z.
Wang
.
2011
.
Phylogeny and divergence times of some racerunner lizards (Lacertidae: Eremias) inferred from mitochondrial 16S rRNA gene segments.
Mol. Phylogenet. Evol.
61
:
400
412
.
Guo,
Z.T.,
B.
Sun
,
Z.S.
Zhang
,
S.Z.
Peng
,
G.Q.
Xiao
,
J.Y.
Ge
,
Q.Z.
Hao
,
Y.S.
Qiao
,
M.Y.
Liang
and
J.F.
Liu
.
2008
.
A major reorganization of Asian climate by the early Miocene.
Clim. Past
4
:
153
174
.
Hartvig,
I.,
M.
Czako
,
E.D.
Kjær
,
L.R.
Nielsen
and
I.
Theilade
.
2015
.
The use of DNA barcoding in identification and conservation of rosewood (Dalbergia spp. ).
PLoS ONE
10
:
e0138231
.
Hebert,
P.D.N.,
A.
Cywinska
,
S.L.
Ball
and
J.R.
Dewaard
.
2003
.
Biological identifications through DNA barcodes.
Proc. R. Soc. Lond. B Biol. Sci.
270
:
313
321
.
Herbert,
T.D.,
K.T.
Lawrence
,
A.
Tzanova
,
L.C.
Peterson
,
R.
Caballero-Gill
and
C.S.
Kelly
.
2016
.
Late Miocene global cooling and the rise of modern ecosystems.
Nat. Geosci.
9
:
843
847
.
Hewitt,
G.
2000
.
The genetic legacy of the Quaternary ice ages.
Nature
405
:
907
913
.
Hillis,
D.M.
and
M.T.
Dixon
.
1991
.
Ribosomal DNA: Molecular evolution and phylogenetic inference.
Q. Rev. Biol.
66
:
411
453
.
Hlinka,
O.,
A.
Murrell
and
S.C.
Barker
.
2002
.
Evolution of the secondary structure of the rRNA internal transcribed spacer 2 (ITS2) in hard ticks (Ixodidae, Arthropoda).
Heredity
88
:
275
279
.
Ji,
Y.J.,
D.X.
Zhang
and
L.J.
He
.
2003
.
Evolutionary conservation and versatility of a new set of primers for amplifying the ribosomal internal transcribed spacer regions in insects and other invertebrates.
Mol. Ecol. Notes
3
:
581
585
.
Jong,
R.D.
2017
.
Fossil butterflies, calibration points and the molecular clock (Lepidoptera: Papilionoidea).
Zootaxa
4270
:
1
63
.
Kalyaanamoorthy,
S.,
B.Q.
Minh
,
T.K.F.
Wong
,
A.V.
Haeseler
and
L.S.
Jermiin
.
2017
.
ModelFinder: Fast model selection for accurate phylogenetic estimates.
Nat. Methods
14
:
587
589
.
Katoh,
K.
and
D.M.
Standley
.
2013
.
MAFFT multiple sequence alignment software version 7: Improvements in performance and usability.
Mol. Biol. Evol.
30
:
772
780
.
Katoh,
T.,
A.
Chichvarkhin
,
T.
Yagi
and
K.
Omoto
.
2005
.
Phylogeny and evolution of butterflies of the genus Parnassius inferences from mitochondrial 16S and ND1 sequences.
Zool. Sci.
22
:
343
351
.
Kumar,
S.,
G.
Stecher
and
K.
Tamura
.
2016
.
MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets.
Mol. Biol. Evol.
33
:
1870
1874
.
Lahaye,
R.,
M.V.D.
Bank
,
D.
Bogarin
,
J.
Warner
,
F.
Pupulin
,
G.
Gigot
,
O.
Maurin
,
S.
Duthoit
,
T.G.
Barraclough
and
V.
Savolainen
.
2008
.
DNA barcoding the floras of biodiversity hotspots.
Proc. Natl. Acad. Sci. USA
105
:
2923
2928
.
Lanfear,
R.,
B.
Calcott
,
S.Y.
Ho
and
S.
Guindon
.
2012
.
PartitionFinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses.
Mol. Biol. Evol.
29
:
1695
1701
.
Lei,
F.M.,
Y.H.
Qu
and
G.
Song
.
2014
.
Species diversification and phylogeographical patterns of birds in response to the uplift of the Qinghai-Tibet Plateau and Quaternary glaciations.
Curr. Zool.
60
:
149
161
.
Li,
J.J.,
X.M.
Fang
,
H.Z.
Ma
,
J.J.
Zhu
,
B.T.
Pan
and
H.L.
Chen
.
1996
.
Geomorphological and environmental evolution in the upper reaches of the Yellow River during the late Cenozoic.
Sci. China Ser. D
39
:
380
390
.
Lu,
H.Y.
and
Z.T.
Guo
.
2014
.
Evolution of the monsoon and dry climate in East Asia during late Cenozoic: A review.
Sci. China Earth Sci.
57
:
70
79
.
McLean,
B.S.,
B.
Nyamsuren
,
A.
Tchabovsky
and
J.A.
Cook
.
2018
.
Impacts of late Quaternary environmental change on the long-tailed ground squirrel (Urocitellus undulatus) in Mongolia.
Zool. Res.
39
:
364
372
.
Meier,
R.,
K.
Shiyang
,
G.
Vaidya
and
P.K.L.
Ng
.
2006
.
DNA barcoding and taxonomy in Diptera: A tale of high intraspecific variability and low identification success.
Syst. Biol.
55
:
715
728
.
Meier,
R.,
G.Y.
Zhang
and
F.
Ali
.
2008
.
The use of mean instead of smallest interspecific distances exaggerates the size of the “barcoding gap” and leads to misidentification.
Syst. Biol.
57
:
809
813
.
Meyer,
C.P.
and
G.
Paulay
.
2005
.
DNA barcoding: Error rates based on comprehensive sampling.
PLoS Biol.
3
:
e422
.
Miao,
Y.F.,
M.
Herrmann
,
F.L.
Wu
,
X.L.
Yan
and
S.L.
Yang
.
2012
.
What controlled Mid-Late Miocene long-term aridification in Central Asia?—Global cooling or Tibetan Plateau uplift: A review.
Earth-Sci. Rev.
112
:
155
172
.
Michel,
F.,
C.
Rebourg
,
E.
Cosson
and
H.
Descimon
.
2008
.
Molecular phylogeny of Parnassiinae butterflies (Lepidoptera: Papilionidae) based on the sequences of four mitochondrial DNA segments.
Ann. Soc. Entomol. Fr.
44
:
1
36
.
Molnar,
P.
and
J.M.
Stock
.
2009
.
Slowing of India's convergence with Eurasia since 20 Ma and its implications for Tibetan mantle dynamics.
Tectonics
28
:
TC3001
.
Morgan,
J.A.T.
and
D.
Blair
.
1998
.
Trematode and monogenean rRNA ITS2 secondary structures support a four-domain model.
J. Mol. Evol.
47
:
406
419
.
Najman,
Y.,
E.
Appel
,
M.
Boudagher-Fadel
,
P.
Bown
,
A.
Carter
,
E.
Garzanti
,
L.
Godin
,
J.T.
Han
,
U.
Liebke
and
G.
Oliver
.
2010
.
Timing of India-Asia collision: Geological, biostratigraphic, and palaeomagnetic constraints.
J. Geophys. Res.
115
:
B12416
.
Nazari,
V.,
E.V.
Zakharov
and
F.A.H.
Sperling
.
2007
.
Phylogeny, historical biogeography, and taxonomic ranking of Parnassiinae (Lepidoptera, Papilionidae) based on morphology and seven genes.
Mol. Phylogenet. Evol.
42
:
131
156
.
Nguyen,
L.T.,
H.A.
Schmidt
,
A.V.
Haeseler
and
B.Q.
Minh
.
2014
.
IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies.
Mol. Biol. Evol.
32
:
268
274
.
Omoto,
K.,
T.
Katoh
,
A.
Chichvarkhin
and
T.
Yagi
.
2004
.
Molecular systematics and evolution of the “Apollo” butterflies of the genus Parnassius (Lepidoptera: Papilionidae) based on mitochondrial DNA sequence data.
Gene
326
:
141
147
.
Omoto,
K.,
T.
Yonezawa
and
T.
Shinkawa
.
2009
.
Molecular systematics and evolution of the recently discovered “Parnassian” butterfly (Parnassius davydovi Churkin, 2006) and its allied species (Lepidoptera, Papilionidae).
Gene
441
:
80
88
.
Owen,
L.A.
and
J.M.
Dortch
.
2014
.
Nature and timing of Quaternary glaciation in the Himalayan-Tibetan orogen.
Quat. Sci. Rev.
88
:
14
54
.
Pérez-Gutiérrez,
M.A.,
A.T.
Romero-García
,
M.C.
Fernández
,
G.
Blanca
,
M.J.
Salinas-Bonillo
and
V.N.
Suárez-Santiago
.
2015
.
Evolutionary history of fumitories (subfamily Fumarioideae, Papaveraceae): An old story shaped by the main geological and climatic events in the Northern Hemisphere.
Mol. Phylogenet. Evol.
88
:
75
92
.
Platas,
G.,
C.
Ruibal
and
J.
Collado
.
2004
.
Size and sequence heterogeneity in the ITS1 of Xylaria hypoxylon isolates.
Mycol. Res.
108
:
71
75
.
Poczai,
P.
and
J.
Hyvönen
.
2010
.
Nuclear ribosomal spacer regions in plant phylogenetics: Problems and prospects.
Mol. Biol. Rep.
37
:
1897
1912
.
Rambaut,
A.
2009
.
FigTree version 1.4: Tree figure drawing tool.
Institute Evol. Biol.
,
Univ. Edinburgh
,
Edinburgh, Scotland
.
http://tree.bio.ed.ac.uk/ (accessed 01 November 2019).
Rambaut,
A.
and
A.
Drummond
.
2016
.
TreeAnnotator version 1.8.3.
http://beast.bio.ed.ac.uk/TreeAnnotator (accessed 01 November 2019).
Rambaut,
A.,
A.J.
Drummond
,
D.
Xie
,
G.
Baele
and
M.A.
Suchard
.
2018
.
Posterior summarization in Bayesian phylogenetics using Tracer 1.7.
Syst. Biol.
67
:
901
904
.
Rampersad,
S.N.
2014
.
ITS1, 5.8 S and ITS2 secondary structure modelling for intra-specific differentiation among species of the Colletotrichum gloeosporioides sensu lato species complex.
SpringerPlus
3
:
684
.
Rebourg,
C.,
F.
Péténian
,
E.
Cosson
and
E.
Faure
.
2006
.
Patterns of speciation and adaptive radiation in Parnassius butterflies.
J. Entomol.
3
:
204
215
.
Ronquist,
F.,
M.
Teslenko
,
P.V.D.
Mark
,
D.L.
Ayres
,
A.
Darling
,
S.
Höhna
,
B.
Larget
,
L.
Liu
,
M.A.
Suchard
and
J.P.
Huelsenbeck
.
2012
.
MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space.
Syst. Biol.
61
:
539
542
.
Rozas,
J.,
A.
Ferrer-Mata
,
J.C.
Sánchez-DelBarrio
,
S.
Guirao-Rico
,
P.
Librado
,
S.E.
Ramos-Onsins
and
A.
Sánchez-Gracia
.
2017
.
DnaSP 6: DNA sequence polymorphism analysis of large data sets.
Mol. Biol. Evol.
34
:
3299
3302
.
Sandel,
B.,
L.
Arge
,
B.
Dalsgaard
,
R.G.
Davies
,
K.J.
Gaston
,
W.J.
Sutherland
and
J.C.
Svenning
.
2011
.
The influence of Late Quaternary climate-change velocity on species endemism.
Science
334
:
660
664
.
Schoville,
S.D.
and
G.K.
Roderick
.
2009
.
Alpine biogeography of Parnassian butterflies during Quaternary climate cycles in North America.
Mol. Ecol.
18
:
3471
3485
.
Shi,
Y.F.
2002
.
Characteristics of late Quaternary monsoonal glaciation on the Tibetan Plateau and in East Asia.
Quat. Int.
97–98
:
79
91
.
Slatkin,
M.
and
W.P.
Maddison
.
1989
.
A cladistic measure of gene flow inferred from the phylogenies of alleles.
Genetics
123
:
603
613
.
Stadler,
T.
2009
.
On incomplete sampling under birth-death models and connections to the sampling-based coalescent.
J. Theor. Biol.
261
:
58
66
.
Su,
C.Y.,
Q.H.
Shi
,
X.Y.
Sun
,
J.Y.
Ma
,
C.X.
Li
,
J.S.
Hao
and
Q.
Yang
.
2017
.
Dated phylogeny and dispersal history of the butterfly subfamily Nymphalinae (Lepidoptera: Nymphalidae).
Sci. Rep.
7
:
8799
.
Sun,
X.J.
and
P.X.
Wang
.
2005
.
How old is the Asian monsoon system?—Palaeobotanical records from China.
Palaeogeogr., Palaeoclimatol., Palaeoecol.
222
:
181
222
.
Wan,
S.M.,
A.C.
Li
,
P.D.
Clift
and
J.B.W.
Stuut
.
2007
.
Development of the East Asian monsoon: Mineralogical and sedimentologic records in the northern South China Sea since 20 Ma.
Palaeogeogr., Palaeoclimatol., Palaeoecol.
254
:
561
582
.
Wright,
S.
1965
.
The interpretation of population structure by F-statistics with special regard to systems of mating.
Evolution
19
:
395
420
.
Wu,
Y.Q.,
Z.J.
Cui
,
G.N.
Liu
,
D.K.
Ge
,
J.R.
Yin
,
Q.H.
Xu
and
Q.Q.
Pang
.
2001
.
Quaternary geomorphological evolution of the Kunlun Pass area and uplift of the Qinghai-Xizang (Tibet) Plateau.
Geomorphology
36
:
203
216
.
Xu,
J.P.
2016
.
Fungal DNA barcoding.
Genome
59
:
913
932
.
Xu,
L.B.,
X.J.
Ou
,
Z.P.
Lai
,
S.Z.
Zhou
,
J.
Wang
and
Y.
Fu
.
2010
.
Timing and style of Late Pleistocene glaciation in the Queer Shan, northern Hengduan Mountains in the eastern Tibetan Plateau.
J. Quat. Sci.
25
:
957
966
.
Yang,
F.S.,
Y.F.
Li
,
X.
Ding
and
X.Q.
Wang
.
2008
.
Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan Plateau and its correlation with the Quaternary climate change.
Mol. Ecol.
17
:
5135
5145
.
Yao,
H.,
J.Y.
Song
,
C.
Liu
,
K.
Luo
,
J.P.
Han
,
Y.
Li
,
X.H.
Pang
,
H.X.
Xu
,
Y.J.
Zhu
and
P.G.
Xiao
.
2010
.
Use of ITS2 region as the universal DNA barcode for plants and animals.
PLoS ONE
5
:
e13102
.
Zachos,
J.,
M.
Pagani
,
L.
Sloan
,
E.
Thomas
and
K.
Billups
.
2001
.
Trends, rhythms, and aberrations in global climate 65 Ma to present.
Science
292
:
686
693
.
Zheng,
B.,
Y.L.
Wang
,
C.C.
Xia
,
D.Y.
Huang
,
Y.
Cao
,
J.S.
Hao
and
C.D.
Zhu
.
2018
.
The complete mitochondrial genome of Parnassius actius (Lepidoptera: Papilionidae: Parnassinae) with the related phylogenetic analysis.
Zool. Syst.
43
:
1
17
.
Zheng,
B.X.,
Q.Q.
Xu
and
Y.P.
Shen
.
2002
.
The relationship between climate change and Quaternary glacial cycles on the Qinghai-Tibetan Plateau: Review and speculation.
Quat. Int.
97–98
:
93
101
.
Zhou,
S.Z.,
X.L.
Wang
,
J.
Wang
and
L.B.
Xu
.
2006
.
A preliminary study on timing of the oldest Pleistocene glaciation in Qinghai-Tibetan Plateau.
Quat. Int.
154–155
:
44
51
.
Zhu,
R.W.,
Y.C.
Li
,
D.L.
Zhong
and
J.Q.
Zhang
.
2017
.
Establishment of the most comprehensive ITS2 barcode database to date of the traditional medicinal plant Rhodiola (Crassulaceae).
Sci. Rep.
7
:
10051
.