Next Article in Journal
New Insights into the Taxonomy of Malacopsylloidea Superfamily (Siphonaptera) Based on Morphological, Molecular and Phylogenetic Characterization of Phthiropsylla agenoris (Malacopsyllidae) and Polygenis (Polygenis) rimatus (Rhopalopsyllidae)
Next Article in Special Issue
Forest Fragmentation and Developmental Stability of Wood Mice Apodemus sylvaticus: A Food-Mediated Effect?
Previous Article in Journal
Diversity and Endemism of Southern African Gekkonids Linked with the Escarpment Has Implications for Conservation Priorities
Previous Article in Special Issue
Determining Extinction for Small Cryptic Species: The Morro Bay Kangaroo Rat
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Phylogeography of the Plateau Pika (Ochotona curzoniae) in Response to the Uplift of the Qinghai-Tibet Plateau

1
School of Geographic Science, Qinghai Normal University, Xining 810008, China
2
Medical College, Qinghai University, Xining 810008, China
3
School of Life Science, Qinghai Normal University, Xining 810008, China
*
Author to whom correspondence should be addressed.
Diversity 2023, 15(2), 307; https://doi.org/10.3390/d15020307
Submission received: 22 December 2022 / Revised: 13 February 2023 / Accepted: 18 February 2023 / Published: 20 February 2023
(This article belongs to the Special Issue Adaptive Evolution and Management in Small Mammals)

Abstract

:
The evolution and current distribution of species on the Qinghai-Tibet Plateau have been significantly impacted by historical occurrences, including the uplift of the plateau and the Quaternary climate upheaval. As a remnant species, the plateau pika (Ochotona curzoniae) is a great model for researching historical events. In this study, 302 samples from 42 sample sites were utilized to analyze the impact of historical events on the evolution and distribution pattern of plateau pikas. The genetic diversity, patterns of differentiation, and historical dynamics of the plateau pika were investigated using molecular markers that included four mitochondrial genes (COI, D-loop, Cytb, and 12S rRNA) and three nuclear genes (GHR, IRBP, and RAG1). The results showed that: (1) The genetic diversity of the plateau pika was high in the Tibetan Plateau (Hd = 0.9997, π = 0.01205), and the plateau pika evolved into five lineages that occupied different geographical areas, with lineage 1 (Group 1) in the south of the Yarlung Zangbo River, lineage 2 (Group 2) in the hinterland of the plateau, lineage 3 (Group 3) in the northeastern part of the plateau, lineage 4 (Group 4) in the Hengduan Mountains, and lineage 5 (Group 5) in the eastern part of the plateau. (2) The gene flow among the five lineages was low, and the differentiation level was high (Nm < 0.25; Fst > 0.25), indicating that the geographical barriers between the five lineages, such as the Yarlung Zangbo River, the Qaidam-Ghuong-Guide Basin, and the Lancang River, effectively promoted the population differentiation of the plateau pika. (3) The plateau pika first spread from the Hengduan Mountains to the entire Qinghai-Tibet Plateau and then conducted small-scale migration and dispersal in several refuges across the plateau in response to climate changes during the glacial and interglacial periods. (4) Except for Group 1 and Group 4, all the other populations exhibited a rapid expansion between 0.06 and 0.01 Mya, but the expansion was considerably delayed or halted by the effects of climate change during the last glacial maximum (0.02 Mya). Overall, the plateau pika on the Qinghai-Tibet Plateau exhibits high genetic diversity, and topographic obstacles, including mountains, valleys, and basins, created by the uplift of the plateau and climatic changes since the Quaternary period have played an important role in the differentiation and historical dynamics of the plateau pika population.

1. Introduction

The Qinghai-Tibet Plateau, which spans an area of over 2.5 million km2 and has an average elevation of 4000 m, is situated in Southwest China [1]. The Himalaya Mountains, Gangdise Mountains, Nianqing Tanggula Mountains, Tanggula Mountains, and Kunlun Mountains are the principal mountain ranges that run from south to north. The Yellow River, Yarlung Zangbo River, Jinsha River, and Salween River are the principal rivers. Global climate change and Cenozoic geological occurrences resulted in the area’s high elevation and complicated terrain. According to Shi et al. [2], the Qinghai-Tibet Plateau experienced two uplifting and leveling cycles in the Tertiary period. The Gangdisi Mountain uplift resulted from the first uplifting, which took place between 45–38 Mya, and the Himalayan Mountain uplift resulted from the second uplifting, which took place 30–21.8 Mya. Following the uplift, the Tibetan Plateau’s average elevation was about 1500–2000 m during this time, but after a protracted period of denudation and leveling, the elevation of the Qinghai-Tibet plateau was only around 1000 m, although at this time it was still a tropical and subtropical environment [3]. The Qinghai-Tibet Plateau experienced a third dramatic uplift and complex climate changes during the late Mesozoic and Quaternary periods. The changes in this period were the main causes of its current complex topography. The third uplift has been divided into three stages by geologists. The Qinghai-Tibet Plateau rose to a height of 3000–3500 m during the first stage of the Qinghai-Tibet Movement (3.4–1.3 Mya) and the second stage of the Kung-Huang Movement (1.1–0.6 Mya). During the third stage, it rose to a height of 4000 m, which is known as the Republic Movement (0.15 Mya to today). During the dramatic uplift of the plateau in the Quaternary period, the plateau experienced several glaciations, which resulted in the gradual succession of the species-rich tropical forest to sparse grassland or even bare land. In addition, the circulation of the glacial and interglacial periods led to active water erosion in the plateau area during this period, and the flows of the Yarlung Zangbo River, Yangtze River, Yellow River, and other rivers were cut down sharply, forming many gorges [2]. These high mountains and deep gorges formed obvious geographical barriers and significantly affected the inheritance and differentiation of the relict species in the plateau areas [4,5,6]. Yu and Zhang [7] summarized the results of phylogeographical studies conducted on 40 plant species and found that most of the plants that are distributed in the northeastern part of the Qinghai-Tibet Plateau experienced an expansion to the plateau hinterland after the last glacial age. The main clades of most of the species that are distributed in the southeastern part of the Qinghai-Tibet Plateau diverged in the early Quaternary due to the uplift of the plateau, and they survived in several refuges on the edge of the plateau or expanded into small areas in response to climate change during the glacial period. Some plants also showed completely different responses to climate change during the ice age. For example, a few hardy plants, such as Potentilla glabra [8], did not retreat but expanded from the eastern part of the plateau to the hinterland under the low-temperature environment during the ice age. Therefore, different species have different responses to geoclimatic events, according to their genetic, physiological, and life-history characteristics. Amphibians, birds, and rodents are the most studied plateau animals. There are three main geographical distribution patterns in this area. (1) For two lineages, including Pseudopodoces humilis [9,10], or blood pheasant (Ithaginis cruentus) [11], and the alpine pygmy frog (Nanorana parkeri) [12], the differentiation was mainly caused by the high mountain barrier and the 400 mm precipitation boundary [12]. (2) The differentiation of four branches, including Eospalax fontanierii and the plateau pika (Ochotona curzoniae) [13,14], occurred mainly due to geographical barriers, such as mountains and basins. (3) Multiple differentiation clades, including Phrynocephalus vlangalii [15], showed no obvious differentiation patterns in the plateau region. Examples of these include the Blanford snow finch (Pyrgilauda blanfordi), the red-necked snow finch (Pyrgilauda ruficollis), the white-rumped blood finch (Onychostruthus taczanowskii), and the Tibetan antelope (Pantholops hodgsonii) [16,17,18,19].
The plateau pika is a small, non-hibernating mammal, also known as the black-lipped pika, that belongs to the family Ochotonidae [20], which is an important high-altitude model animal that is broadly distributed in the area 3000–6000 m above sea level [21]. Studies have shown that the survival footprint of the pika in China can be traced back to the early Eocene, about 53 Mya ago [22]. It was a relict species in the Quaternary glacial period and is an ideal organism for studying the relationship between historical events and biological evolution [23]. Most research on the lineage differentiation and distribution patterns of the pika has focused on the North American pika (Ochotona princeps) and the spotted necked pika (Ochotona collaris) [24,25,26,27,28].
In recent years, with the further development of pedigree geography, the differentiation of pikas on the Qinghai-Tibet Plateau has also received extensive attention. For example, Yu et al. [29], who constructed a phylogenetic tree using mitochondrial genes (Cytb and ND4), found that pikas could be differentiated into three subgenera: Pika, Ochotona, and Conothoa. In the work of Niu et al. [30], the phylogenetic relationships of 27 species within the genus Ochotona were reconstructed based on the mitochondrial cytochrome b gene, with results showing that pikas are divided into five major species groups: the northern group, the surrounding Qinghai-Tibet Plateau group, the Qinghai-Tibet Plateau group, the Huanghe group, and the Central Asia group. Melo-Ferreira et al. [31] studied 12 nuclear genes from 11 representative species of pikas and concluded that there were four subgenera of pikas. Koju et al. [32] used two mitochondrial (Cytb and COI) and five nuclear gene segments (RAG1, RAG2, TTN, OXAIL, and IL1RAPL1) to infer the phylogeny of 14 taxa, and their results supported the O. syrinx group as a distinct lineage beyond the four recognized subgenera, leading the authors to conclude that there are five subgenera of pikas. In the same year, Liu et al. [33] conducted a comprehensive study of 27 existing species of pikas in China, combining morphological and molecular markers, and their results also supported the claim that the pikas are divided into five subgenera, namely, Alienauroa, Conothoa, Ochotona, Lagotona, and Pika. A study by Lissovsky [34] represents the most systematic molecular phylogenetic study of this topic. He consulted pika specimens from major museums in Europe and the United States, carried out molecular systematic studies, and combined these findings with morphological information, and believed that there were 28 species of pika in the world and performed consistent taxonomic revisions for all widely distributed species. In addition to the above studies on the intra-genus differentiation of pikas, a few scholars have also studied the intraspecific differentiation of Plateau pikas using molecular markers. For instance, Ci et al. [35] divided 32 plateau pika populations into 10 units according to geographical barriers such as mountains, rivers, and basins, and used mitochondrial genes to study the intra-specific differentiation of these 10 units. They found that the genetic diversity of these 10 units was high, and they were divided into six groups under the influence of climate and geographical barriers during the ice age. The central group was located in the hinterland of the plateau, with the five surrounding groups located around it. Then, He et al. [14] constructed a Bayesian phylogenetic tree using mitochondrial genes and divided 37 geographical populations of plateau pikas into four lineages, namely, east, west, north, and south. They found that the east, west, and north lineages had expanded in the past, whereas the southern lineages (located south of the Yarlung Zangbo River) had experienced a sharp decline. They speculated that this may have been related to the climate changes that occurred during the Quaternary glacial and interglacial periods.
Currently, there are many studies available on the differentiation of the pika family but few studies have focused on the evolution within the plateau pika population. Moreover, the existing intra-species research has been carried out using mitochondrial genes as molecular markers because mitochondrial genes are widely used in molecular phylogeography due to their unique maternal genetic characteristics. However, Avise et al. [36] suggested that since nuclear genes offer richer genetic information, they should be gradually adopted, and research that combines the use of multiple nuclear gene loci represents a promising direction for future phylogeographical research. Therefore, some scholars, such as Chen et al. [37], have studied the phylogeography of Sigmella biguttata using mitochondrial (COI, COII, and ND1) and nuclear genes (ITS). In addition, Jeremy et al. [38] studied the origins of family members in New Zealand using mitochondrial and nuclear genes. Moreover, Ding et al. [39] used mitochondrial (COI, COIII, Cytb, and D-loop) and nuclear genes (vWF) to study the intraspecific and intermediate differentiation of Tibetan hamsters. However, no studies hve been conducted on the differentiation of the plateau pika population using both nuclear and mitochondrial genes. To fill this gap, in this study, we used four mitochondria and three nuclear genes to investigate the intraspecies differentiation characteristics of the plateau pika. We aimed to determine the environmental factors that drove the plateau pika’s evolution and diversity by understanding the genetic diversity, differentiation, and gene flow between the different populations; the findings of this study will provide a scientific basis for the genetic management and resource conservation of the plateau pika population.

2. Materials and Methods

2.1. Sample Sources

A total of 302 samples (Table 1) were collected from 42 sample sites based on their geographical distances and geographical barriers (mountains, rivers, and basins) (Figure 1). Each individual was anesthetized on site after measuring its basic biological information, such as its body length and weight, and leg muscle tissue was removed with a sterile scalpel. The sample was placed in an Eppendorf tube and stored in liquid nitrogen.

2.2. Genomic DNA Extraction, Polymerase Chain Reaction Amplification, and Sequencing

The genomic DNA was extracted from the plateau pika samples using a TaKaRa DNA Kit, and its concentration and integrity were detected using 2% agar-agar gel electrophoresis. The samples were stored at −80 °C. The primer design and sources are shown in Table 2. The polymerase chain reaction (PCR) volume was 50 μL, which contained 5.0 μL 10 × Taq buffer (Takara, Beijing China), 2.5 U Taq DNA polymerase, 1.5 mM MgCl2, 0.1 mM deoxynucleotide triphosphates, 0.4 uM of each primer, and 0.5 μL total DNA. The rest of the volume was made up to 50 μL with sterilized ultra-purified water. The PCR conditions were as follows: predenaturation at 95 °C for 8 min; 35 cycles of denaturation at 95 °C for 45 s, annealing at 52 °C–64 °C for 45 s, and extension at 72 °C for 1 min; and a final extension at 72 °C for 7 min. The PCR products were photographed on a 1% agar gel to detect the target fragment, and they were sent to Sangon Biotech (Shanghai) Co., Ltd. for purification and sequencing.

2.3. Population Genetic Structure and Divergence

The nucleic acid sequences that were obtained via sequencing were spliced using DNAMAN 7 [14], and MEGA 7 [40] was used for sequence alignment and correction and to determine the base composition statistics. DnaSP 5.0 [41] was used to calculate the variable sites, parsimony-informative sites, the number of haplotypes (Nh), haplotype diversity (Hd), and nucleotide diversity (π).
SAMOVA 1.0 [42] was used for spatial analysis of the molecular variance, and the K value was set at 2–42. Each run was repeated 100 times. Then, an analysis of molecular variance (AMOVA) was conducted in Arlequin 3.5 [43] to estimate the distribution of the population’s genetic distance and genetic variation and to calculate the population’s genetic differentiation index (Fst). The maximum gene flow (Nm) was calculated using Formula (1) [37,44].
N m = 1 F s t 4 F s t
Govindajaru [45] divided gene flow into three levels according to the Nm value: when Nm > 1, the gene flow among the various populations is considered to be high; when 0.25 < Nm < 0.99, it is at a medium level; and when Nm < 0.25, it is low.

2.4. Phylogenetic Analysis and Haplotype Network

The North American pika (Ochotona princeps) and European rabbit (Oryctolagus cuniculus) were used as the outgroups. The maximum likelihood (ML) and Bayesian inference (BI) methods were used to construct a phylogenetic tree with 302 tandem DNA sequences (12S rRNA, COI, D-loop, Cytb, GHR, IRBP, and RAG1). The maximum likelihood phylogenetic tree was implemented with MEGA 7. The best nucleotide replacement models (12S rRNA, COI, D-loop, Cytb, and RAG1: GTR + I + G; GHR: GTR + G; IRBP: HKY + I) were identified using the corrected Akaike information criterion [46], which was implemented in MrModeltest 3.7 [47]. Bayesian analysis was conducted using MrBayes 3.2.6, which was run with twenty million generations of Markov chain Monte Carlo (MCMC) simulation and sampled every 1000 generations. MCMC was run using the default model parameters, starting from a random tree. The first 25% were discarded as a conservative burn-in, and the remaining samples were used to generate a 50% majority-rule consensus tree. Here, a Bayesian posterior probability equal to or above 0.95 was considered to indicate strong relationships [48]. Subsequently, the tree file was visualized using FigTree 1.4.7 (http://tree.bio.ed.ac.uk/software/figtree/ accessed on 5 October 2022). Finally, the haplotype network diagram was constructed using a median-joining network by means of PopART 1.7 [49].

2.5. Nucleotide Mutation Rate and Estimation of the Divergence Time

To acquire as accurate a divergence time as possible, the number of nucleotide substitutions per site (d) was estimated basd on comparisons between the focal species and outgroup species using the following Formula (2):
d = (tv + tvR)/m
where tv is the number of transversions between Ochotona curzoniae and the outgroup taxa, R is the transition/transversion ratio within the Ochotona curzoniae complex, and m is the sequence length [50,51]. The transition and transversion values were calculated in the program MEGA 7.0 [52]. The rates of nucleotide substitutions per site, per lineage, and per year were calculated using the formula λ = d/2T when an estimate of d was obtained, where T is the divergence time between the ingroup and outgroup species [51]. The mutation rate per nucleotide site and per generation was calculated with the formula μ = λg, where g is the generation time (g = 0.46 years for Ochotona curzoniae) [53,54]. In this study, d was 0.0172 (tv = 17, R = 6.0, m = 6909) for the combined genes (12S rRNA, COI, D-loop, Cytb, GHR, IRBP, and RAG1). The rate of nucleotide substitution per site, per lineage, and per year (λ) was about 0.86 × 10−8 (T = 10.65 Mya, obtained from Rui et al. [55]), and the mutation rate per generation (μ) was about 0.396% Mya−1 for the concatenated DNA sequence of Ochotona curzoniae.
The divergence times among the different populations of Ochotona curzoniae were estimated using BEAST 2.4.7 [56]. There were 6069 bp in the combined genes (12S rRNA, COI, D-loop, Cytb, GHR, IRBP, and RAG1). The best base replacement model was the GTR + I + G model, and the molecular clock model used was the strict clock model. The Markov chain Monte Carlo (MCMC) analysis was run twice, the chain lengths were 50 million, with sampling occurred every 1000 generations, and the first 10% was discarded as burn-in. The program TRACER 1.7 [57] was used to test the validity of the sampling results for MrBayes and BEAST (effective sample size [ESS] > 200). Finally, the tree files were visualized using iTOL (https://itol.embl.de/itol.cgi accessed on 15 October 2022).

2.6. Demographic History

Arlequin 3.5 [43] was used to calculate the mismatch distribution, Tajima’s D, and Fu’s Fs values to check whether the population had experienced historical expansion. The goodness of fit between the nucleotide mismatch distribution and the expected distribution under the population expansion model was detected using the sum of squared deviation (SSD) and Hapending’s raggedness index (Hri). The Bayesian MCMC method in BEAST 2.4.7 [56] was used to speculate regarding the time of divergence for the different pedigrees, and the historical population dynamics were analyzed for the combined plateau pika DNA sequences using a Bayesian skyline plot (BSP). In this study, a strict clock was used in BEAST and the substitution rate was based on the estimated results of the nucleotide mutation rate, as described above. The same nucleotide substitution model was used for the Bayesian phylogenetic analysis. The BSP was run twice for MCMC chain lengths of 100 million generations, it was sampled every 1000 generations, and both the MCMC convergence and the ESS (ESS > 200) were tested using the program TRACER 1.7 [57]. Finally, R 3.6.3 was used to determine the composition.

3. Results

3.1. Analysis of the Gene Sequence Variation

For the seven gene sequences in this study, the variation analysis results are shown in Table 3. Among them, the COI gene exhibited the most total variation, with 150 variable loci, and the IRBP gene exhibited the least total variation, with 20 variable loci. The D-loop gene had the least single mutation sites (2), whereas the GHR gene had the most single mutation sites (26). In addition, the C + G content in the IRBP gene sequence was the highest (61.6%), and in the other genes, it was lower than the A + T content. The average conversion/transposition values for the seven genes were all greater than one. There were also only a few IRBP gene haplotypes, whereas the D-loop gene had the most haplotypes, with totals of 22 and 108, respectively. The concatenated sequence was 6909 bp in length and contained 773 mutation sites, 690 of which were parsimony-informative sites and 83 were single mutation sites. The base composition had a larger percentage of A + T (52.5%) than that of C + G (47.6%). Additionally, there were 290 different haplotypes, and the π and Hd values were 0.9997 and 0.01205, respectively.

3.2. Population Genetic Structure Analysis

For the different geographical populations (Table 1), the nucleotide diversity of eight populations (Guide-2, Heka, Qumalai, Maduo, Anduo, AJikehu, Jiangxhu, and Yushu) was significantly higher than that of the other populations (π > 0.007). The SAMOVA analysis based on the total sequences showed that Fct increased between K = 2 and 5 and decreased between K = 6 and 42. When K = 5, Fct reached its maximum value (Fct = 0.50). Therefore, the best groups for plateau pika population differentiation were determined to be Group 1–Group 5 (Table 4) as the genetic diversity of these five lineages was relatively high (Hd > 0.9; π > 0.004).
The results of the AMOVA (Table 5) showed that the genetic variation of the combined plateau pika sequences without grouping revealed significant differentiation among the populations (61.71% of the variation) when compared with the level within the populations (38.29%). The genetic variation among the groups (49.84%) was greater than that within the populations (16.66%) and among the populations (33.50%) when the haplotypes were divided into five groups.
Two key indices were used to gauge the degree of population differentiation, namely, Fst and Nm. When 0 < Fst < 0.05, the degree of differentiation among the populations was small; when 0.05 < Fst < 0.15, it was moderate; when 0.15 < Fst < 0.25, it was high; and when Fst > 0.25, it was extremely high [58]. As shown in Figure 2a, the pairwise Fst among the lineages ranged from 0.554 to 0.798, which was greater than 0.25, and they all reached a significant level (p < 0.05). Furthermore, the pairwise Fst values among lineage 1 (Group 1) were larger than those of the other groups (all greater than 0.7). When compared with the other groups, the Fst value between Group 4 and Group 5 (0.554) was the smallest. Furthermore, the range of Fst values among the 42 geographical populations was 0.187–0.979 (Figure 2b and Appendix A, Figure A1), and the Fst values among the geographical populations within the other lineages were all small except for those of lineage 1 (Group 1) and lineage 4 (Group 4). The statistical results for the Nm values showed that the Nm values among the five lineages ranged from 0.063 to 0.201 (Figure 2c), and this indicates that the maximum gene flow between the five lineages was low. The maximum gene flow within the geographical populations in lineages 2 (Group 2), 3 (Group 3), and 5 (Group 5) was larger, whereas the maximum gene flow within the geographical populations in lineages 1 (Group 1) and 4 (Group 4) was smaller than that in the other lineages (Figure 2d and Appendix A, Figure A2). Moreover, the pairwise Fst values within and among eight populations (Ajikehu, Anduo, Jiangxigou, Guide-2, Qumalai, Yushu, Maduo, and Heka) and various populations were compared with various groups, and they were determined to be small. In addition, the Nm values were all large (Figure 2d and Appendix A, Figure A2), indicating that there was a large amount of gene flow among these eight populations and the other lineages.

3.3. Phylogenetic Tree and Haplotype Network Diagram Analysis Results

The topological structures of the phylogenetic trees that were obtained by means of the ML and BI methods were the same. The 42 geographical populations of the plateau pika were clustered into five lineages (Figure 3a), which was consistent with the SAMOVA analysis results. Among these, lineage 1 diverged the earliest and was separated from the other lineages, whereas lineage 4 and lineage 5 were the last to diverge and had the smallest genetic distances between them. The phylogenetic tree analysis results were consistent with those of the haplotype network diagram (Figure 3b), where the haplotype in Group 4 was located in the center of the haplotype network diagram. Except for the geographical populations of Zhidao and Qumalai, which had three shared haplotypes, the other haplotypes were unique to each geographical population, indicating a high degree of differentiation among the geographical populations.

3.4. Analysis of the Historical Population Dynamics

3.4.1. Neutrality Test

The neutral test results (Table 6) showed that the Tajiam’s D and Fu’s Fs values of Group 2, Group 3, and Group 5 were all negative. Therefore, this indicates that they have experienced significant historical expansion. This is consistent with the results of the BSP curve (Figure 4b) and the single-peak diagram of the nucleotide mismatch distribution (Figure 4c). The BSP curve showed that the populations of Group 2, Group 3, and Group 5 began to expand about 0.06 Mya ago, and the expansion slowed down about 0.02 Mya. Group 3 expanded dramatically for the second time about 0.01 Mya ago. Then, the Tajima’s D and Fu’s Fs values for Group 1 were both greater than zero, the nucleotide mismatch distribution did not show an obvious single peak, and the BSP curve showed no evidence of expansion. Both the Tajima’s D and Fu’s Fs values for Group 4 were negative, but they were not significant. The distribution of the nucleotide mismatch took the form of a multi-peak graph, and the BSP curve did not indicate an obvious expansion trend. Therefore, this indicates that the plateau pika of Group 1 and Group 4 did not undergo any obvious historical expansion.

3.4.2. Estimation of the Historical Divergence Time

As shown in Figure 3a, the plateau pika population experienced four large differentiation events historically. Firstly, the populations south of the Yarlung Zangbo River (Group 1) underwent differentiation around 0.72 Mya (95% confidence interval HPD: 0.88–0.86 Mya). The second differentiation event occurred about 0.48 Mya (95% confidence interval HPD: 0.39–0.57 Mya), and the populations in the central Tibetan Plateau were differentiated (Group 2). The population in the northeastern Qinghai-Tibet plateau (Group 3) was estimated to have diverged 0.47 Mya (95% confidence interval HPD: 0.38–0.56 Mya). Additionally, the fourth differentiation occurred about 0.23 Mya (95% confidence interval HPD: 0.18–0.28 Mya) on both sides of the Lancang River (Group 4 and Group 5).

4. Discussion

4.1. Effects of the Qinghai-Tibet Plateau Geographic Barrier on the Genetic Structure and Differentiation of the Plateau Pika Population

Genetic diversity is an important indicator of the responses of biological populations to environmental change and human disturbance [59]. Intraspecific genetic diversity determines evolutionary capacity. Haplotype diversity and nucleotide diversity are important parameters that can be used to measure the genetic diversity of a population [60]. The results based on the combined sequences showed that the Hd of the plateau pika was 0.9997, which was greater than 0.5, and π was 0.01205, which was greater than 0.5%, indicating that the genetic diversity of the plateau pika was high on the Tibetan Plateau [61]. This is consistent with the research results of He [14], Ci [35], and others. On the one hand, the complex topography of the Qinghai-Tibet Plateau has fragmented the habitat of the plateau pika, and specific environmental factors in the different habitats could have driven the unique variation in its genetics and its evolution [62]. On the other hand, these findings might have been obtained because the plateau pika has a sizable number of geographical populations, with polygyny being its main form of reproduction [63], and dispersal behavior is common [64]. Dispersal improves the likelihood of mating across various families within a single geographic population, increasing genetic diversity. However, due to the relatively small dispersal range of the plateau pika, which basically remains within the range of 106–178 m2 [65], the haplotype analysis results showed that there were only three shared haplotypes between Guide-2 and Tongren (Figure 2). The other 40 geographical populations did not share haplotypes, and each population had a high proportion of unique haplotypes. Thus, these diverse genetic resources are of positive significance to the future evolution of the plateau pika population.
The SAMOVA analysis and the phylogenetic tree showed that the 42 populations were divided into five lineages, which occupied different geographical areas (Figure 1). Lineage 1 (Group 1) was located south of the Brahmaputra River, Group 2 was in the hinterland of the Qinghai-Tibet Plateau, Group 3 was in the northeastern Tibetan Plateau, Group 4 was in the eastern part of the Tibetan Plateau, and Lineage 5 (Group 5) was in the Hengduan Mountain region. The AMOVA showed that the genetic variation was mainly identified within the lineages, at 49.84%, whereas the genetic variation within the populations was only 16.66%, and the genetic variation among the different geographical populations accounted for 33.5%. The Fst and Nm values among the different lineages were all greater than 0.25, and the Nm values were all less than 0.25. The results showed that there was great genetic differentiation among the five lineages [66]. We hypothesized that the external cause of this variation and differentiation was mainly due to geographical isolation, which hindered gene exchange among the populations and resulted in the directional evolution of each lineage in each region.
As shown in Figure 1, the largest geographical barrier between Lineage 1 and the other lineages was the Yarlung Zangbo River, which was formed as early as the Late Miocene (about 8 Mya) [67] and which has one of the highest elevations for a river in China. The largest geographical barrier between Group 3 and the other lineages was the Qaidam-Gonghe-Guide Basin, which has an average elevation of 2600–9000 m and has very little vegetation coverage due to low precipitation and the soil texture. Therefore, it is not suitable for plateau pika survival, and it is the main geographical barrier that blocks gene exchange between Group 3 and the other lineages. The main geographical barrier between Group 4 and Group 5 was the Lancang River. The Lancang River was formed as early as 17 Mya [68], originating in Zadoi County, Qinghai Province. The valley has an elevation of 1000–4500 m and a width of 150–1200 m. It is a typical mountain-gorge landform [69]. Additionally, both the marshes in the river source area and the valleys in the middle and lower reaches are not suitable for plateau pikas, so they have become the biggest geographical barrier that blocks gene exchange between lineage 4 and lineage 5. The barrier effects of the Yarlung Zangbo River and the Qaidam-Gonghe-Guide Basin are reflected in the studies of He [14] and Ci [35] et al., and these may occur because the pika is sensitive to low altitudes and high temperatures [27]. This is because low-altitude river valleys and basins have relatively high oxygen contents and temperatures, and this is inconsistent with the physiology of the plateau pika, which is adapted to survival at low altitudes and low temperatures and has developed over a long period [70]. Thus, these conditions result in a lack of gene exchange across the valleys and basins. Therefore, the barrier effect on plateau pika dispersal was prominent.
In contrast to the findings of studies by He and Ci et al. [14,35], the population differentiation of the plateau pika between the Hengduan Mountains (Group 4) and the eastern Qinghai-Tibet Plateau (Group 5) was obvious in this study, which may have occurred because those authors used mitochondrial genes to analyze population differentiation. In this study, both mitochondrial and nuclear genes were used as molecular markers to analyze the differentiation of the plateau pika population. Mitochondrial genes are maternal molecular markers, whereas nuclear genes are parental molecular markers, and plateau pika males are the main dispersers [62], so the results obtained when using both mitochondrial and nuclear genes were different from those obtained using mitochondrial genes alone. The Fst and Nm values (Figure 2) showed that there was a large amount of gene flow among the populations of lineages 1 (Group 1), 3 (Group 3), and 5 (Group 5), and the differentiation level was low. It is possible that the plateau pika underwent repeated dispersal, centered in a few refuges within each lineage in the Quaternary glacial and interglacial periods, leading to frequent gene exchange between the diverse populations within the lineage. The findings also indicated that geographical barriers such as the Tanggula, Nianqing Tanggula, and Kunlun Mountains, which are in the hinterland of the plateau, and the Bayan Har Mountains, Animaqing Snow Mountain, and Yellow River, which are in the eastern margin of the plateau, did not have significant effects on the isolation of the plateau pika. The gene flow among the four geographic populations within Group 4 was significantly smaller than that among the other three groups, probably because the populations of Group 4 were located in the Hengduan Mountain region, which has more mountains and valleys and more complex terrain. This could have restricted the dispersal of the plateau pika and resulted in less gene exchange among the geographic populations within Group 4. If this isolation continues to play a role, more differentiation will likely occur. This result is consistent with the fact that the Hengduan Mountain region, because of its favorable geographical and climatic conditions, led to rapid biological evolution and high biodiversity.

4.2. Phylogeny and Divergence Time of the Plateau Pika

In a haplotype network diagram, the original haplotype is usually positioned in the center of the stele structure [71]. In Figure 2b, the oldest lineage in the Qinghai-Tibet Plateau was Group 4, which was located in the Hengduan Mountains, a hotspot for biodiversity research. This area is also the place of origin of many species [72]. The Hengduan Mountains may have been the place of origin of the plateau pika, with the species then spreading over the whole plateau, gradually forming the existing distribution pattern. In the early stage of its formation (about 22 Mya), the Yarlung Zangbo River was just a tributary of the Red River [73]. The riverbed was shallow and did not fully form until about 8 Mya, and it experienced several periods of low flow due to drought. Therefore, the plateau pika may have undergone trans-river dispersal through some of the shallow areas during the low-flow period of the Yarlung Zangbo River [74]. Later, with the rise of the Tibetan Plateau, the Yarlung Zangbo River system continuously expanded to form a large canyon [2], blocking the gene flow of the plateau pika. Consequently, the lineage south of the Yarlung Zangbo River (Group 1) first differentiated about 0.72 Mya, which is a similar situation to that implied by the results of He et al. [14]. After several uplifting movements, the plateau’s hinterland rose to 3000–3500 m about 0.62 Mya [2]. The eastern region of the plateau has a relatively low elevation due to the small impact of the uplifting movement. Then, lineage 2 (Group 2) differentiated around 0.48 Mya. Concurrently, the drought intensified in the northwest region of the plateau, and the Qaidam Basin and Gonghe Basin were successively formed, which blocked the gene exchange of the plateau pikas that were distributed on both sides of the basin, resulting in the differentiation of lineage 3 (Group 3) around 0.47 Mya in the northeast region of the Plateau. In the Quaternary (about 2.48 Mya), the plateau uplifted further and experienced repeated cycles of glaciation and interglaciation [2]. These great changes in geology and climate further expanded the major water systems of the Jinsha River, Nujiang River, and Lancang River in the Hengduan Mountain region [68]. In addition, the further uplift of the major mountains, such as Thanyuntaweng, Mangkang, and Nianqing Tanggula, led to the continuous fragmentation of the plateau pika habitat in the Hengduan Mountains, and, finally, Group 4 and Group 5 were divided into two lineages about 0.23 Mya.

4.3. Historical Dynamics of the Plateau Pika Population

Tajima’s D value, Fu’s Fs value, the nucleotide mismatch distribution, SSD, Hri test results, and BSP diagrams are widely used in the analysis of population history dynamics. In general, when the distribution of the nucleotide mismatch presents a smooth single peak, the SSD and Hri tests do not deviate significantly from the population expansion model and the Tajima’s D and Fu’s Fs values are significantly negative. This is considered evidence of a historical population expansion [75,76,77,78,79]. In this study, the neutral testing and BSP analysis showed that Group 2, Group 3, and Group 5 expanded dramatically about 0.06 Mya, whereas Group 1 and Group 4 did not show similar fluctuations. This period was the interglacial period before the last glacial maximum (about 0.06–0.03 Mya) [2,80], where the temperature in the Tibetan Plateau rose, glaciers melted, the amount of vegetation increased, and the suitable area for plateau organisms increased. Consequently, during this time, the Group 2, Group 3, and Group 5 populations all experienced tremendous growth. In addition, the temperature of the Tibetan Plateau during this period was 1 °C higher than that of the Holocene Megathermal Period according to Yafeng et al. [2]. However, due to the low latitude and the annual average temperature in the southern region of the Qinghai-Tibet Plateau, which can be as high as 4–5 °C [81], the temperature in the southern region of the plateau may have reached 5–6 °C. However, the optimal annual average temperature for the plateau pika is −1–0 °C [82], and the species is sensitive to high temperatures [83,84]. Thus, Group 1 and Group 4, which are located in the southern region of the Qinghai-Tibet plateau, may have had a reduced fitness index due to the high temperatures at this time, so they did not undergo expansion like the other lineages. This result was similar to that of Ci et al. [35].
Yujiao et al. [85] showed that other plateau pika lineages, except for those south of the Yarlung Zangbo River, expanded rapidly between 0.076–0.017 Mya; however, after 0.017 Mya, the plateau pika population rapidly declined. During this period, the Qinghai-Tibet Plateau entered the peak of the ice age [2]. The temperature was 9 ℃ lower than that of modern times, the precipitation was only 30–70% of that in modern times, the climate deteriorated, the glacier extent increased, and the vegetation biomass decreased, restricting the expansion of the plateau pika. The results of this study showed that during this period (about 0.02 Mya), the plateau pika did not exhibit a downward trend but the expansion trend slowed down significantly or even stopped temporarily (Figure 4b), indicating that the climate change during the ice age impacted the survival of the plateau pika. At this time, the numbers of most of the plateau species declined, and the organisms in the hinterland of the plateau retreated to the refuges in the east of the plateau [6] but the cold-tolerant and drought-tolerant organisms, such as the golden dew plum and red sand, began to expand during the ice age [86,87,88].
The plateau pika is also a hardy creature, and it may have used its special physiology [70] to resist the cold, which is why it did not decline significantly during the glacial period. This may also have occurred because the populations of each lineage had their own refuges, such as Amdo and Ajikehu in the region where Group 2 is located. Jinguangou and Gude-2 in the region of Group 3 and Yushu, Heka, Maduo, and Qumalai in the region of Group 5 may also have been refuges for each lineage due to their high haplotypes and higher genetic diversity than those for the other regions for the same lineages [89,90]. Similar findings were also found in the genealogies of other plants and animals. For example, Chen et al. [91] studied the genetic evolution of Rhodiola and found that the northern and southern Tanggula Mountains and Hengduan Mountains may be refuges. The Amdo area in this study is located in the southern region of the Tanggula Mountains. Hippophae tibetana [92] also has several refuges in the western, central, and southeastern parts of the plateau; Potentilla glabra [8] has refuges in the upper reaches of Lancang River and Hengduan Mountains. Moreover, the plateau zokor [13] also has four refuges near the origin of the Lancang River, the northeastern area of the Qinghai-Tibet Plateau, and the eastern part of the Qinghai-Tibet Plateau.
The results of the molecular variance analysis in this study also supported the notion that the existing distribution pattern of the plateau pika was formed due to its dispersal into multiple refuges. If the plateau pika retreated to the eastern part of the Qinghai-Tibet Plateau or the Hengduan Mountains during the ice age and gradually spread from the east to the hinterland of the plateau after the ice age (like other organisms) [93], the populations that are located further away from the eastern plateau and Hengduan Mountains would have decreased genetic diversity due to the long-distance dispersal [94] and the gene flow among the populations, and other geographic populations would be small. However, the molecular variance results in this study showed that the gene flow among the geographical populations within each lineage was large, and the gene flow among the different lineages was small. In particular, lineages 1, 2, and 3 had a higher differentiation level and smaller gene flow than that between lineages 4 and 5, which are located in the eastern part of the plateau and the Hengduan Mountain Range, indicating that there was still isolation among these lineages during the Quaternary glacial and interglacial cycles. Therefore, we speculate that the plateau pika first spread to the entire Qinghai-Tibet Plateau from the Hengduan Mountain Range. Then, small-scale back-migration and dispersal were carried out to multiple refuges in various regions of the plateau in response to climatic cycles during the glacial and interglacial periods.

5. Conclusions

In this study, we found that the genetic diversity of the plateau pika was generally high on the Tibetan Plateau (Hd = 0.9997; π = 0.01205), and the plateau pika evolved into five lineages that occupied different geographical areas, namely, lineage 1 (Group 1) in the south of the Yarlung Zangbo River, lineage 2 (Group 2) in the hinterland of the plateau, lineage 3 (Group 3) in the northeastern part of the plateau, lineage 4 (Group 4) in the Hengduan Mountains, and lineage 5 (Group 5) in the eastern part of the plateau. Our results also proved that the geographical barriers between the five lineages, such as the Yarlung Zangbo River, the Qaidam-Ghuong-Guide Basin, and the Lancang River, effectively promoted the population differentiation of the plateau pika. Other geographical barriers, such as the Jinsha River and Yellow River; Tanggula, Tanggula Nianqing, Kunlun, and Bayan Har Mountains; and Animaqing Snow Mountain, were less effective. We also found that the plateau pika first spread from the Hengduan Mountains to the entire Qinghai-Tibet Plateau and then carried out small-scale migration and dispersal into multiple refuges in various regions of the plateau in response to climate changes during the glacial and interglacial periods.
Except for Group 1 and Group 4, all the other populations underwent a rapid expansion between 0.06 and 0.01 Mya but their expansion was considerably delayed or halted by the effects of climate change during the last glacial maximum (0.02 Mya). In conclusion, the plateau pika in the Qinghai-Tibet Plateau exhibits high genetic diversity, and topographical obstacles including mountains, valleys, and basins created by the uplift of the plateau and climatic changes since the Quaternary have played an important role in the differentiation and historical dynamics of the plateau pika population.

Author Contributions

Y.Q. conducted the research, analyzed the data, and wrote the paper; X.P. made suggestions to this paper; Z.L. and D.S. helped to process the data; Z.C. guided the research and performed extensive updating of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This research study was funded by Qinghai Provincial Key Laboratory of Medicinal Animal and Plant Resources of the Qinghai-Tibetan Plateau.

Institutional Review Board Statement

The animals involved in the study were handled in accordance with the National Regulations on the Administration of Laboratory Animals (GB14923-2010).

Data Availability Statement

All data and materials are available upon request.

Acknowledgments

We are particularly indebted to Su Xu from Qinghai Normal University and Delin Qi from Qinghai University for their constructive suggestions regarding an earlier draft of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Figure A1. The fixation index (FST) results for 42 populations. The diagonal indicates the level of significance, * p < 0.05; ** p < 0.01; *** p < 0.001; -: not significant.
Figure A1. The fixation index (FST) results for 42 populations. The diagonal indicates the level of significance, * p < 0.05; ** p < 0.01; *** p < 0.001; -: not significant.
Diversity 15 00307 g0a1
Figure A2. The maximum gene flow (Nm) results for 42 populations.
Figure A2. The maximum gene flow (Nm) results for 42 populations.
Diversity 15 00307 g0a2

References

  1. Zhang, Y.L.; Li, B.Y.; Zheng, D. A discussion on the boundary and area of the Tibetan Plateau in China. Geogr. Res. 2002, 21, 1–9. (In Chinese) [Google Scholar] [CrossRef]
  2. Shi, Y.F.; Li, J.J.; Li, B.Y. Uplift and environmental changes of Qinghai-Tibetan Plateau in the Late Cenozoic. Guangzhou Guangdong Sci. Technol. Press 1999, 54, 10–20. [Google Scholar]
  3. Li, J.J.; Weng, S.X.; Zhang, Q.S.; Wang, F.B.; Zheng, B.X.; Li, B.Y. Discussion on the age, amplitude and form of the uplift of Qinghai-Tibet Plateau. Sci. China 1979, 22, 1314–1328. [Google Scholar]
  4. Liu, Q.; Chen, P.; He, K.; Kilpatrick, C.W.; Liu, S.-Y.; Yu, F.-H.; Jiang, X.-L. Phylogeographic Study of Apodemus ilex (Rodentia: Muridae) in Southwest China. PLoS ONE 2012, 7, e31453. [Google Scholar] [CrossRef] [Green Version]
  5. Qu, J.; Liu, N.; Bao, X.; Wang, X. Phylogeography of the ring-necked pheasant (Phasianus colchicus) in China. Mol. Phylogenet. Evol. 2009, 52, 125–132. [Google Scholar] [CrossRef]
  6. Qiu, Y.X.; Fu, C.X.; Comes, H.P. Plant molecular phylogeography in china and adjacent regions:tracing the genetic imprints of Quaternary cliamte and environmental change in the world’s most diverse temperate flora. Mol. Phylogenet. Evol. 2011, 59, 225–244. [Google Scholar] [CrossRef]
  7. Yu, H.B.; Zhang, Y.L. Advances in Phylogeography of Alpine Plants in the Tibetan Plateau and Adjacent Regions. Acta Bot. Boreal. Occident. Sin. 2013, 33, 1268–1278. [Google Scholar]
  8. Wang, L.Y.; Ikeda, H.; Liu, T.L.; Wang, Y.J.; Liu, J.Q. Repeated range expansion and glacial endurance of potentilla glabra (rosaceae) in the qinghai-tibetan plateau. J. Integr. Plant Biol. 2009, 51, 698–706. [Google Scholar] [CrossRef]
  9. Schafer, E. Ornithologische Ergebnisse zweier Forschungsreisen nach Tibet. J. Ornithol. 1938, 86, 7–340. [Google Scholar] [CrossRef]
  10. Chen, F.G.; Luo, S.Y. Fauna Sinica. Aves: Passeriformes: Bombycillidae-Prunellidae; Science Press: Beijing, China, 1998; pp. 178–182. [Google Scholar]
  11. Zhan, X.J.; Zheng, Y.F.; Wei, F.W.; Bruford, M.W.; Jia, C.X. Molecular evidence for Pleistocene refugia at the eastern edge of the Tibetan Plateau. Mol. Ecol. 2011, 20, 3014–3026. [Google Scholar] [CrossRef]
  12. Zhou, W.W.; Zhang, B.L.; Chen, H.M.; Jin, J. DNA Barcodes and Species Distribution Models Evaluate Threats of Global Climate Changes to Genetic Diversity: A Case Study from Nanorana parkeri (Anura: Dicroglossidae). PLoS ONE 2014, 9, e103899. [Google Scholar] [CrossRef] [Green Version]
  13. Tang, L.Z.; Wang, L.Y.; Cai, Z.Y.; Zhang, T.Z.; Ci, H.X.; Lin, G.H.; Su, J.P. Allopatric divergence and phylogeographic structure of the plateau zokor (Eospalax baileyi), a fossorial rodent endemic to the Qinghai-Tibetan Plateau. J. Biogeogr. 2010, 37, 657–668. [Google Scholar] [CrossRef]
  14. He, Y.J. Study on Population Genetic Structure and Demographic History of Ochotona curzoniae. Ph.D. Thesis, Northwest Institute of Plateau Biology, UCAS, Xining, China, 2018. [Google Scholar]
  15. Jin, Y.T.; Brown, R.P.; Liu, N.F. Cladogenesis and phylogeography of the lizard phrynocephalus vlangalii (Agamidae) on the Tibetam plateau. Mol. Ecol. 2008, 17, 1971–1982. [Google Scholar] [CrossRef]
  16. Cramp, S.; Perrins, C.M. Handbook of the Birds of Europe, the Middle East and North Africa: The Birds of the Western Palearctic. Vol 8: Crows to Finches; Oxford University Press: New York, NY, USA, 1994; pp. 52–66. [Google Scholar]
  17. Fu, T.S.; Song, Y.J.; Gao, W. Fauna Sinica: Aves; Science Press: Beijing, China, 1998; pp. 165–175. [Google Scholar]
  18. Gebauer, A.; Kaiser, M. Biology and behavior of general Asiatic snow finches Montifringilla and mountain-steppe sparrows Pyrgilauda. J. Ornithol. 1994, 135, 55–57. [Google Scholar] [CrossRef]
  19. Zhang, F.F.; Jiang, Z.G.; Xu, A.C.; Yan, Z. Recent geological events and intrinsic behavior influence the population genetic structure of the chiru and Tibetan gazelle on the Tibetan plateau. PLoS ONE 2013, 8, e60712. [Google Scholar] [CrossRef] [Green Version]
  20. Smith, A.T.; Badingqiuying; Wilson, M.C.; Hogan, B.W. Functional-trait ecology of the plateau pika Ochotona curzoniae in the Qinghai-Tibetan Plateau ecosystem. Integr. Zool. 2019, 14, 87–103. [Google Scholar] [CrossRef] [Green Version]
  21. Huan, L.; Rui, Z.; Jianxiao, Z.; Xiaodan, H.; Jiapeng, Q. Environmental ltering increases with elevation for the assembly of gut microbiota in wild pikas. Microb. Biotechnol. 2019, 12, 976–992. [Google Scholar] [CrossRef] [Green Version]
  22. Rose, K.D.; DeLeon, V.B.; Missiaen, P.; Rana, R.S.; Sahni, A.; Singh, L.; Smith, T. Early Eocene lagomorph (Mammalia) from Western India and the early diversification of Lagomorpha. Proc. R. Soc. Lond. Ser. B Biol. Sci. 2008, 275, 1203–1208. [Google Scholar] [CrossRef] [Green Version]
  23. Angelone, C. Family Ochotonidae (Lagomorpha) and its application in biochronology: Some case studies from the Plio-Quaternary of Eurasia. Quat. Int. 2008, 179, 5–8. [Google Scholar] [CrossRef]
  24. Galbreath, K.E.; Hafner, D.J.; Zamudio, K.R. When Cold Is Better: Climate-Driven Elevation Shifts Yield Complex Patterns of Diversification and Demography in an Alpine Specialist (American Pika, Ochotona princeps). Evolution 2009, 63, 2848–2863. [Google Scholar] [CrossRef]
  25. Galbreath, K.E.; Hoberg, E.P. Return to Beringia: Parasites reveal cryptic biogeographic history of North American pikas. Proc. R. Soc. B Biol. Sci. 2012, 279, 371–378. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Lanier, H.C.; Massatti, R.; He, Q.; Olson, L.E.; Knowles, L.L. Colonization from divergent ancestors: Glaciation signatures on contemporary patterns of genomic variation in collared pikas (ochotona collaris). Mol. Ecol. 2015, 24, 3688–3705. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Castillo, J.A.; Epps, C.W.; Davis, A.R.; Cushman, S.A. Landscape effects on gene flow for a climate-sensitive montane species, the american pika. Mol. Ecol. 2014, 23, 843–856. [Google Scholar] [CrossRef]
  28. Hayley, C.L.; Link, E.O. Deep barriers, shallow divergences: Reduced phylogeographical structure in the collared pika (Mammalia: Lagomorpha: Ochotona collaris). J. Biogeogr. 2013, 40, 466–478. [Google Scholar] [CrossRef]
  29. Yu, N.; Zheng, C.; Zhang, Y.-P.; Li, W.-H. Molecular systematics of pikas (Genus Ochotona) inferred from mitochondrial DNA sequences. Mol. Phylogenet. Evol. 2000, 16, 85–95. [Google Scholar] [CrossRef]
  30. Niu, Y.; Wei, F.; Li, M.; Liu, X.; Feng, Z. Phylogeny of pikas (Lagomorpha, Ochotona) inferred from mitochondrial cytochrome b sequences. Folia Zool.-Praha 2004, 53, 141–156. [Google Scholar]
  31. Melo-Ferreira, J.; Matos, A.L.; Areal, H.; Lissovsky, A.A.; Carneiro, M.; Es-teves, P.J. The phylogeny of pikas (Ochotona) inferred from a multilocus coalescent approach. Mol. Phylogenet. Evol. 2015, 84, 240–244. [Google Scholar] [CrossRef] [Green Version]
  32. Koju, N.P.; He, K.; Chalise, M.K.; Ray, C.; Chen, Z.; Zhang, B.; Wan, T.; Chen, S.; Jiang, X. Multilocus approaches reveal underestimated species diversity and inter-specific gene flow in pikas (Ochotona) from southwestern China. Mol. Phylogenet. Evol. 2017, 107, 239–245. [Google Scholar] [CrossRef]
  33. Liu, S.; Jin, W.; Liao, R.; Sun, Z.; Zeng, T.; Fu, J.; Liu, Y.; Wang, X.; Li, P.; Tang, M.; et al. Phylogenetic study of Ochotona based on mitochondrial Cytb and morphology with a description of one new subgenus and five new species. Acta Theriol. Sin. 2017, 37, 1–43. [Google Scholar] [CrossRef]
  34. Lissovsky, A.A.; Yatsentyuk, S.P.; Obolenskaya, E.V.; Koju, N.P.; Ge, D. Diversification in highlands: Phylogeny and taxonomy of pikas of the subgenus Conothoa (Lagomorpha, Ochotonidae). Zool. Scr. 2022, 51, 267–287. [Google Scholar] [CrossRef]
  35. Ci, H.X.; Lin, G.H.; Cai, Z.Y.; Tang, L.Z.; Su, J.P.; Liu, J.Q. Population history of the plateau pika endemic to the Qinghai-Tibetan Plateau based on mtDNA sequence data. J. Zool. 2009, 279, 396–403. [Google Scholar] [CrossRef]
  36. Avise, J.C. Phylogeography: Retrospect and prospect. J. Biogeogr. 2009, 36, 3–15. [Google Scholar] [CrossRef] [Green Version]
  37. Chen, R.; Niu, L.K.; Duanmu, H.N.; Wang, Z.Q.; Che, Y.L. Natural barriers and climatic oscillation in the Quaternary Pleistocene influence the phylogeography of Sigmella biguttata (Blattodea: Ectobiidae) revealed by mitochondrial and nuclear genes. Acta Entomol. Sin. 2022, 65, 235–245. [Google Scholar] [CrossRef]
  38. Jeremy, B.S.; Paul, M.J.; Islam, G.; Mark, I.S.; Eleanor, P.J.; Chrissen, E.C.G.; Carolyn, M.K. The diverse origins of New Zealand house mice. Mice Proc. R. Soc. B 2009, 276, 209–217. [Google Scholar] [CrossRef] [Green Version]
  39. Ding, L.; Liao, J. Phylogeography of the Tibetan hamster Cricetulus kamensis in response to uplift and environmental change in the Qinghai-Tibet Plateau. Ecol. Evol. 2019, 9, 7291–7306. [Google Scholar] [CrossRef] [Green Version]
  40. Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 2016, 33, 1870–1874. [Google Scholar] [CrossRef] [Green Version]
  41. Librado, P.; Rozas, J. DnaSPv5: A software for comprehensive analysisi of DNA polymorphism data. Bioinformatics 2009, 25, 1451–1452. [Google Scholar] [CrossRef] [Green Version]
  42. Dupanloup, I.; Schneider, S.; Excoffier, L. A simulated annealing approach to define the genetic structure of populations. Mol. Ecol. 2002, 11, 2571–2581. [Google Scholar] [CrossRef]
  43. Excoffier, L.; Lischer, H.E.L. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Eol. Ecol. Resour. 2010, 10, 564–567. [Google Scholar] [CrossRef]
  44. Wright, S. The genetical structure of populations. Ann. Eugen 1949, 15, 323–354. [Google Scholar] [CrossRef]
  45. Govindara, J.; Diddahally, R. Estimates of gene flow in forest trees. Biol. J. Linn. Soc. 1989, 37, 345–357. [Google Scholar] [CrossRef]
  46. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control 1974, 19, 716–723. [Google Scholar] [CrossRef]
  47. Posada, D.; Crandall, K.A. MODELTEST: Testing the model of DNA substitution. Bioinformatics 1998, 14, 817–818. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Leaché, A.D.; Reeder, T.W. Molecular systematics of the Eastern Fence lizard (Sceloporus undulatus): A comparison of parsimony, likelihood, and Bayesian approaches. Syst. Biol. 2002, 51, 44–68. [Google Scholar] [CrossRef]
  49. Leigh, J.W.; Bryant, D. POPART: Full-feature software for haplotype network construction. Methods Ecol. Evol. 2015, 6, 1110–1116. [Google Scholar] [CrossRef]
  50. Nei, M. Age of the common ancestor of human mitochondrial DNA. Mol. Biol. Evol. 1992, 9, 1176–1178. [Google Scholar] [CrossRef] [Green Version]
  51. Rooney, A.P.; Honeycutt, R.L.; Derr, J.N. Historical population size change of bowhead whales inferred from DNA sequence polymorphism data. Evolution 2001, 55, 1678–1685. [Google Scholar] [CrossRef]
  52. Tamura, K.; Stecher, G.; Peterson, D.; Filipski, A.; Kumar, S. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 2013, 30, 2725–2729. [Google Scholar] [CrossRef] [Green Version]
  53. Zhang, C.D. The Reproductive Characteristics of Plateau Pika. J. Northeast For. Univ. 2001, 29, 90–92. [Google Scholar] [CrossRef]
  54. Jiapeng, Q.; Ming, L.; Min, Y.; Yanming, Z.; Weihong, J. Reproduction of plateau pika (Ochotona curzoniae) in Guoluo, Qinghai-Tibetan Plateau. Eur. J. Wildl. Res. 2012, 58, 269–277. [Google Scholar] [CrossRef]
  55. Rui, X.T.; Jiao, W.; Yi, F.L.; Cheng, R.Z.; Guan, L.M.; Feng, J.L.; Yue, L.; Megan, P.; Podsiadlowski, L.; Yan, Y.; et al. Genomics and morphometrics reveal the adaptive evolution of pikas. Zool. Res. 2022, 43, 813–826. [Google Scholar] [CrossRef]
  56. Drummond, A.J.; Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 2007, 7, 214. [Google Scholar] [CrossRef] [Green Version]
  57. Rambaut, A.; Drummond, A.J.; Xie, D.; Baele, G.; Suchard, M.A. posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 2018, 67, 901–904. [Google Scholar] [CrossRef] [Green Version]
  58. Hudson, R.R.; Slatkin, M.; Maddison, W.P. Estimation of levels of gene flow from DNA Sequence data. Geneties 1992, 132, 583–589. [Google Scholar] [CrossRef]
  59. Reed, D.; Frankham, R. The correlation between population fitness and genetic diversity. Conserv. Biol. 2003, 17, 230–237. [Google Scholar] [CrossRef]
  60. Vrijenhoek, R.C. Genetic Diversity and Fitness in Small Populations. Conserv. Genet. 1994, 68, 37–53. [Google Scholar] [CrossRef]
  61. Durand, J.D.; Tsigenopoulos, C.S.; Unlu, E.; Berrebi, P. Phylogeny and biogeography of the family cyprinidae in the middle east inferred from cytochrome b DNA—Evolutionary significance of this region. Mol. Phylogenet. Evol. 2002, 22, 91–100. [Google Scholar] [CrossRef]
  62. Pearson, R.G.; Dawson, T.P. Predicting the impacts of climate change on the distribution of species: Are bioclimate envelope models useful? Glob. Ecol. Biogeogr. 2003, 12, 361–371. [Google Scholar] [CrossRef] [Green Version]
  63. Wei, W.H.; Fan, N.C.; Zhou, W.Y.; Yang, S.M.; Cao, Y.F. Aggressive behavior of Plateau pika during breeding period. Curr. Zool. 2000, 46, 278–286. [Google Scholar]
  64. Qu, J.P.; Chen, Q.Q.; Zhang, Y.M. Behaviour and reproductive fitness of postdispersal in plateau pikas (Ochotona curzoniae) on the Tibetan Plateau. Mammal Res. 2018, 63, 151–159. [Google Scholar] [CrossRef]
  65. Qu, J.P.; Li, K.; Yang, M.; Li, W.; Zhang, Y.M.; Smith, A.T. Seasonal dynamic pattern of spacial territory in social groups of plateau pikas (Ochotona curzoniae). Acta Theriol. Sin. 2007, 27, 215–220. [Google Scholar]
  66. Balloux, F.; Lugon, M.N. The estimation of population differentiation with microsatellite markers. Mol. Ecol. 2002, 11, 155–165. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  67. Song, Z.; Wan, S.; Colin, C.; Yu, Z.; Li, A. Paleoenvironmental evolution of south asia and its link to himalayan uplift and climatic change since the late eocene. Glob. Planet. Chang. 2021, 200, 103459. [Google Scholar] [CrossRef]
  68. Nie, J.S.; Ruetenik, G.; Gallagher, K.; Hoke, G.; Garzione, C.N.; Wang, W.T.; Stockli, D.; Hu, X.F.; Wang, Z.; Wang, Y.; et al. Rapid incision of the Mekong River in the middle Miocene linked to monsoonal precipitation. Nat. Geosci. 2018, 11, 944–948. [Google Scholar] [CrossRef]
  69. Li, H.G.; Zhang, P.Q.; Guan, Z. Analysis of runoff evolution characteristics in the upper watershed of Lancang River in recent 30 years. J. Yangtze River Sci. Res. Inst. 2022, 12, 1–7. [Google Scholar] [CrossRef]
  70. Bai, Z.Z.; Wuren, T.; Liu, S.; Han, S.R.; Chen, L.; Mc, C.D.; Ge, R.L. Intermittent cold exposure results in visceral adipose tissue “browning” in the plateau pika (Ochotona curzoniae). Comp. Biochem. Physiol. Part A Mol. Integr. Physiol. 2015, 184, 171–178. [Google Scholar] [CrossRef]
  71. Bao, W.Y.; Zhang, Y.; Lin, P.C.; Nan, P.; Huang, Y.Y.; Jin, H.F.; Zhong, Y. Phylogeography of Gymnadenia conopsea from the Qinghai-Tibet Plateau. Biotechnol. Bull. 2016, 32, 96–102. [Google Scholar] [CrossRef]
  72. Zhang, R.Z.; Zheng, C.L. The geographical distribution of mammals and the evolution of mammals fauna in Qinghai-Xizang Plateau. Acta Geogr. Sin. 1985, 40, 225–231. [Google Scholar]
  73. Clark, M.K.; Schoenbohm, L.M.; Royden, L.H.; Whipple, K.X.; Burchfiel, B.C.; Zhang, X.; Tang, W.; Wang, E.; Chen, L. Surface uplift, tectonics, and erosion of eastern Tibet from large-scale drainage patterns. Tectonics 2004, 23, 6.1–6.20. [Google Scholar] [CrossRef] [Green Version]
  74. Wu, Y.; Long, D.; Lall, U.; Bridget, R.S.; Fu, Q.T.; Xu, D.F.; Jian, S.Z.; Jian, Y.Z.; Hao, W.; Chun, H.H. Reconstructed eight-century streamflow in the Tibetan Plateau reveals contrasting regional variability and strong nonstationarity. Nat. Commun. 2022, 13, 6416. [Google Scholar] [CrossRef]
  75. Tajima, F. The effect of change in population size on DNA polymorphism. Genetics 1989, 123, 597–601. [Google Scholar] [CrossRef]
  76. Tajima, F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 1989, 123, 585–595. [Google Scholar] [CrossRef]
  77. Fu, Y.X.; Li, W. Statistical tests of neutrality of mutations. Genetics 1993, 133, 693–709. [Google Scholar] [CrossRef]
  78. Fu, Y.X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 1997, 147, 915–925. [Google Scholar] [CrossRef]
  79. Rogers, A.R.; Harpending, H. Population growth makes waves in the distribution of pairwise genetic differences. Mol. Biol. Evol. 1992, 9, 552–569. [Google Scholar] [CrossRef]
  80. Schafer, J.M.; Tschudi, S.; Zhao, Z.; Wu, X.; Ivy, O.S.; Wieler, R.; Baur, H.; Kubik, P.W.; Schluchter, C. The limited influence of glaciations in Tibet on global climate over the past 170,000 yr. Earth Planet Sci. Lett. 2002, 194, 287–297. [Google Scholar] [CrossRef]
  81. Shi, Y.F.; Kong, Z.C.; Wang, S.M.; Tang, L.Y.; Wang, R.B.; Yao, T.D.; Zhao, X.T.; Zhang, P.Y.; Shi, S.H. Climate and environment at the peak of the Holocene Megathermic Period in China. Sci. Sin. Chim. 1993, 8, 865–873. (In Chinese) [Google Scholar]
  82. Qi, Y.; Pu, X.; Li, Y.; Li, D.; Huang, M.; Zheng, X.; Guo, J. Prediction of Suitable Distribution Area of Plateau pika (Ochotona curzoniae) in the Qinghai-Tibet Plateau under Shared Socioeconomic Pathways (SSPs). Sustainability 2022, 14, 12114. [Google Scholar] [CrossRef]
  83. Smith, A.T. Distribution and dispersal of pikas: Influences of behavior and climate. Ecology 1974, 55, 1368–1376. [Google Scholar] [CrossRef]
  84. Smith, A.T.; Weston, M.L. Ochotona princeps. Mamm. Species 1990, 352, 1–8. [Google Scholar] [CrossRef] [Green Version]
  85. Yujiao, H.; Gonghua, L.; Haixin, C.; Cuixia, L.; Jianping, S. The past population dynamics of ochotona curzoniae and the response to the climate change. North-West. J. Zool. 2018, 14, 220–225. Available online: http://biozoojournals.ro/nwjz/content/v14n2/nwjz_e171704_Yujiao.pdf (accessed on 22 October 2022).
  86. Sun, Y.; Ikeda, H.; Wang, Y.; Liu, J. Phylogeography of potentilla fruticosa (rosaceae) in the qinghai-tibetan plateau revisited: A reappraisal and new insights. Trans. Bot. Soc. Edinb. 2010, 3, 249–257. [Google Scholar] [CrossRef]
  87. Li, Z.H.; Chen, J.; Zhao, G.F.; Guo, Y.P.; Kou, Y.X.; Zheng, M.Y. Response of a desert shrub to past geological and climatic change: A phylogeographic study of reaumuria soongarica (tamaricaceae) in western china. J. Syst. Evol. 2012, 50, 11. [Google Scholar] [CrossRef]
  88. Kropf, M.; Kadereit, J.W.; Comes, H.P. Differential cycles of range contraction and expansion in Europe an high mountain plants during the Late Quaternary: In sights from Pritzelago alpina (L.) O. Kuntze (Brassicaceae). Mol. Ecol. 2003, 12, 931–949. [Google Scholar] [CrossRef] [PubMed]
  89. Comes, H.P.; Kadereit, J.W. The effect of quaternary climatic changes on plant distribution and evolution. Trends Plant Sci. 1998, 3, 432–438. [Google Scholar] [CrossRef]
  90. Hewitt, G. The genetic legacy of the Quaternary ice ages. Nature 2000, 405, 907–913. [Google Scholar] [CrossRef]
  91. Chen, S.L.; Qing, B.G.; De, J.Z.; Sheng, Y.C.; Yi, Z.D.; Fa, Q.Z.; Ying, H.L. Chloroplast DNA Phylogeography of Rhodiola alsia (Crassulaceae) in the Qinghai—Tibet Plateau. Botany 2009, 87, B09–B059. [Google Scholar]
  92. Hao, W.; Laqiong, K.S.; Fan, L.; Yuguo, W.; Zhiping, S.; Qianhong, W.; Jiakuan, C.; Wenju, Z. Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan Plateau. Mol. Ecol. 2010, 19, 2964–2979. [Google Scholar] [CrossRef]
  93. Li, Z.H.; Zhang, Q.; Liu, J.Q.; Tomas, K.; Martin, L. The pleistocene demography of an alpine juniper of the Qinghai-Tibetan Plateau:tabula rasa, cryptic refugia or something else? J. Biogeogr. 2011, 38, 31–43. [Google Scholar] [CrossRef]
  94. Qing, B.G.; De, J.Z.; Yi, Z.D.; Fa, Q.Z.; Yi, H.L.; Peng, C.F.; Shi, L.C. Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Bot. J. Linn. Soc. 2012, 168, 204–215. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Map of sample locations for plateau pikas.
Figure 1. Map of sample locations for plateau pikas.
Diversity 15 00307 g001
Figure 2. Fst and Nm values among populations and groups of plateau pika based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1). (a) Fst value of the group; the diagonal indicates the level of significance, ** p < 0.01; *** p < 0.001; (b) Fst values of 42 populations, G1–G5: Group 1–Group 5; (c) Nm value of group; (d) Nm value of 42 populations, G1–G5: Group 1–Group 5.
Figure 2. Fst and Nm values among populations and groups of plateau pika based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1). (a) Fst value of the group; the diagonal indicates the level of significance, ** p < 0.01; *** p < 0.001; (b) Fst values of 42 populations, G1–G5: Group 1–Group 5; (c) Nm value of group; (d) Nm value of 42 populations, G1–G5: Group 1–Group 5.
Diversity 15 00307 g002
Figure 3. Maximum likelihood (ML) phylogenetic tree (a) and median-joining network (b) based on the combined dataset (12S rRNA + COI + D-loop + Cytb + GHR + IRBP + RAG1) of plateau pikas. Values near the node indicate the posterior probability values and bootstrap values of the phylogenetic tree. Short bars in the haplotype network represent mutation steps and the black dots represent missing haplotypes. The circle name and color correspond to different regions, and the haplotype circle size denotes the number of sampled individuals.
Figure 3. Maximum likelihood (ML) phylogenetic tree (a) and median-joining network (b) based on the combined dataset (12S rRNA + COI + D-loop + Cytb + GHR + IRBP + RAG1) of plateau pikas. Values near the node indicate the posterior probability values and bootstrap values of the phylogenetic tree. Short bars in the haplotype network represent mutation steps and the black dots represent missing haplotypes. The circle name and color correspond to different regions, and the haplotype circle size denotes the number of sampled individuals.
Diversity 15 00307 g003
Figure 4. Historical dynamic analysis of plateau pika population. (a) Phylogenetic time-tree of plateau pika populations. (b) Bayesian skyline plots (BSPs) for population dynamics with the time of plateau pikas in the middle and lower reaches; dotted lines indicate median values and solid lines indicate 95% upper and lower confidence intervals. (c) Mismatch distribution analysis diagram of each group of plateau pika populations.
Figure 4. Historical dynamic analysis of plateau pika population. (a) Phylogenetic time-tree of plateau pika populations. (b) Bayesian skyline plots (BSPs) for population dynamics with the time of plateau pikas in the middle and lower reaches; dotted lines indicate median values and solid lines indicate 95% upper and lower confidence intervals. (c) Mismatch distribution analysis diagram of each group of plateau pika populations.
Diversity 15 00307 g004
Table 1. Sampling information and genetic parameters based on the combined dataset (Cytb + COI + 12S rRNA + D-loop + RAG1 + IRBP + GHR) of plateau pika geographic populations.
Table 1. Sampling information and genetic parameters based on the combined dataset (Cytb + COI + 12S rRNA + D-loop + RAG1 + IRBP + GHR) of plateau pika geographic populations.
GroupsSample Point
Code
PopulationLongitudeLatitudeAltitude (m)Sample
Size
Number of Haplotypes
Nh
Haplotype
Diversity
Hd
Nucleotide
Diversity
π
Group 1LKZLangkazi90.41729.1094425331.0000.00518
JZJiangzi90.10128.9014660750.9050.00077
Group 2BDQBudongquan39.89735.5224610661.0000.00491
TTHTuotuohe92.4434.2164540771.0000.00407
NQUNaqu91.79731.284600551.0000.0028
ADAnduo91.71832.1574799771.0000.00876
MZGKMozhugongka92.29629.693443412121.0000.00181
BLHBeiluhe92.94234.8624590551.0000.00515
T1Tibet-187.21829.2374481551.0000.00167
T3Tibet-382.56330.5784944320.6670.00048
T4Tibet-485.08929.4934607540.9000.00055
NMNimo90.2729.5023908661.0000.00205
GLDDGeladandong91.65233.5894873441.0000.0023
KEQKKaerqiuka90.75537.0434184551.0000.00477
AJKHAjikehu88.6137.0034250551.0000.01159
TZHTuzihu87.30836.84750551.0000.00466
Group 3MYMenyuan101.27537.69326015151.0000.0062
MLMole100.29937.963379015140.9900.00235
TJTianjun99.10637.245337011111.0000.00496
QLEBQilianebao100.93437.9683429771.0000.00223
QLARQilianarou100.52538.0483031221.0000.00129
RSReshui100.43437.5483520551.0000.00276
GCGangcha100.13437.3253370540.9000.00216
NDNiaodao99.75837.1713158331.0000.00556
JXGJiangxigou100.21136.6213157551.0000.02028
GD1Guide-1102.06737.23725221.0000.00575
GD2Guide-2101.20536.2543650661.0000.00854
Group 4NQNangqian96.50832.19362012121.0000.00452
BSBasu97.20630.6744490551.0000.00109
YLSYelashan97.29530.1874338551.0000.00104
BDBangda97.12930.5294348441.0000.00065
Group 5TRTongren101.71635.5863813551.0000.00247
ZKZeku101.4735.056369012121.0000.0028
KAHeka99.90835.8213890551.0000.00818
MQMaqin100.21234.505372012121.0000.00446
GanDGande100.21834.203421011111.0000.00295
ABAba101.58133.009344013120.9870.00113
SQShiqu98.04732.984440010101.0000.00421
YSYushu96.88633.057384015151.0000.01348
ZDZhiduo95.69633.9394170991.0000.00305
QMLQumalai95.87734.1394390760.9520.00784
MDMaduo98.13334.796425011111.0000.00753
Table 2. Primers used in this study.
Table 2. Primers used in this study.
GenesPrimer Sequences (5′-3′)Amplified Size (bp)SourceAnnealing Temperature (°C)
COIH:ACTACTGGCTTCAATCTACTTCTC
L:AAGACATAGAGGTTATGGAGTTGG
1705Ding, 2020 [34]64
12S rRNAF:AAAGCAAAACACTGAAAATG
R:TTTCATCTTTTCCTTGCGGTA
1149He, 2019 [15]64
CytbH:CGGAATTCCATTTTTGGTTTACAAGAC
L:CGAAGCTTGATATGAAAAACCATCGTTG
1170Kocher, 1989
[36]
64
D-loopF:ATGTTCCGCCCAATCAGCCAAT
R:GTTGCTGGTTTCACGGAGGATGG
655Ci, 2009 [31]59
GHRF:TCAGCCACAGAGGTTAGAAGG
R: CACATAGCCACACGATGAGAG
665This study52
RAG1F:CCGACACCACCAACATTCAA
R:CCTTCACATCTCCACCTTCTTC
1069This study52
IRBPH:GTCCTCTTGGATAACTACTGCTT
L:CTCCACTGCCCTCCCATGTCT
744Ding, 2020 [34]61
Note: 12S rRNA, 12S ribosomal RNA gene; COI, cytochrome oxidase subunit 1 gene; Cytb, cytochrome b gene; D-loop, control region; IRBP, interphotoreceptor retinoid-binding protein gene; RAG1, recombination activating gene 1 gene; GHR, growth hormone receptor gene.
Table 3. Analysis of different gene sequences and total sequence variation.
Table 3. Analysis of different gene sequences and total sequence variation.
Gene12S rRNACytbCOID-loopGHRIRBPRAG1Combined
Sequence
Sequence length/bp1093115817056565737449806909
Variable sites114150277143402029773
Singleton variable sites26882265883
Parsimony informative sites88142269141141521690
Average conversion/transpose value2.610.75.56.32.69.45.46
Haplotype738092108282272290
Haplotype/Sample number/%24.226.530.535.89.37.323.896
A/%35.426.526.73423.219.828.728.1
T/%2326.228.52620.118.722.924.4
C/%24.93427.829.932.430.824.228.8
G/%16.613.31710.124.230.824.118.8
Nuclear genes: GHR, IRBP, RAG1. Mitochondrial genes: 12S rRNA, Cytb, COI, D-loop.
Table 4. Genetic diversity of each group of plateau pika populations based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1).
Table 4. Genetic diversity of each group of plateau pika populations based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1).
GroupPopulations
Included
Number of
Haplotypes
Nh
Haplotype
Diversity
Hd
Nucleotide
Diversity
π
Group 1LKZ JZ80.9560.01027
Group 2NM TZH AJKH KEQK
T1 T3 T4 BLH MZGK
AD NQU TTH BDQ GLDD
780.9990.00625
Group 3MY ML GD1 TJ ND GC RS QLAR QLEB JXG GD2680.9990.00708
Group 4BS BD YLS NQ261.0000.00432
Group 5MQ MD GanD TR ZK KA AB SQ YS ZD QML1141.0000.00699
For population information, see Table 1. The same below.
Table 5. AMOVA analysis results of the populations of plateau pika based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1).
Table 5. AMOVA analysis results of the populations of plateau pika based on the combined dataset (12S rRNA + COI + D-loop +Cytb + GHR + IRBP + RAG1).
Source of VariationD.fSum of SquaresVariance ComponentsVariation
Percentage
Fixation
Indices
One groupAmong populations428578.15826.7191561.71Fst = 0.61714
Within populations2614326.29616.5758538.29
Total30312,904.45443.29499
Five groups based on
lineages
Among groups45752.32524.74701va49.84Fsc = 0.33220
Among populations within groups372774.0288.27352vb16.66Fst = 0.66503
Within populations2604324.29616.63191vc33.50Fct = 0.49840
Table 6. Neutral test and mismatch analysis of different lineages of plateau pika in the Qinghai-Tibet Plateau based on the combined dataset (12S rRNA + COI + D-loop + Cytb + GHR + IRBP + RAG1).
Table 6. Neutral test and mismatch analysis of different lineages of plateau pika in the Qinghai-Tibet Plateau based on the combined dataset (12S rRNA + COI + D-loop + Cytb + GHR + IRBP + RAG1).
GroupTajima’s DFu’s FsSSDHarpending’s Raggedness Index rT (Ma)
Group 10.98023 ns4.00851 ns0.89400 ns0.083333 ns-
Group 2−1.14529 *−24.19805 *0.84891 ns0.00086444 ns0.06
Group 3−1.32911 *−24.41101 *0.84603 ns0.00164939 ns0.06
Group 4−0.83933 ns−8.9991275 ns0.85944 ns0.01231735 ns-
Group 5−1.65000 *−24.16632 *0.88954 *0.00142879 ns0.06
SSD: Sum of squared deviations; ns: not significant; * p < 0.05; -: no increase.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Qi, Y.; Pu, X.; Li, Z.; Song, D.; Chen, Z. Phylogeography of the Plateau Pika (Ochotona curzoniae) in Response to the Uplift of the Qinghai-Tibet Plateau. Diversity 2023, 15, 307. https://doi.org/10.3390/d15020307

AMA Style

Qi Y, Pu X, Li Z, Song D, Chen Z. Phylogeography of the Plateau Pika (Ochotona curzoniae) in Response to the Uplift of the Qinghai-Tibet Plateau. Diversity. 2023; 15(2):307. https://doi.org/10.3390/d15020307

Chicago/Turabian Style

Qi, Yinglian, Xiaoyan Pu, Zhilian Li, Daoguang Song, and Zhi Chen. 2023. "Phylogeography of the Plateau Pika (Ochotona curzoniae) in Response to the Uplift of the Qinghai-Tibet Plateau" Diversity 15, no. 2: 307. https://doi.org/10.3390/d15020307

APA Style

Qi, Y., Pu, X., Li, Z., Song, D., & Chen, Z. (2023). Phylogeography of the Plateau Pika (Ochotona curzoniae) in Response to the Uplift of the Qinghai-Tibet Plateau. Diversity, 15(2), 307. https://doi.org/10.3390/d15020307

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop