Population Structure of Plateau Pika based on Mitochondrial DNA and Microsatellite Analysis

He Yujiao1,2, Xie Huichun1, Lin Gonghua2,3, Zhang Tongzuo2,3, Su Jianping2,3 and Du Yurong1*

1Key Laboratory of Medicinal Plant and Animal Resources in Qinghai-Tibetan Plateau, College of Life Science, Qinghai Normal University, Xining, Qinghai, China

2Key Laboratory of Adaptation and Evolution of Plateau Biota, Northwest Institute of Plateau Biology, Chinese Academy of Sciences, Xining, Qinghai, China

3Qinghai Provincial Key Laboratory of Animal Ecological Genomics, Northwest Institute of Plateau Biology, Chinese Academy of Sciences, Xining, Qinghai, China

ABSTRACT

The plateau pika (Ochotona curzoniae) is widely distributed across the Qinghai-Tibetan Plateau (QTP) and has experienced the quaternary climate and geological events of the plateau. In this study, we assessed the population patterns and the differentiation time in the plateau pika using mitochondrial DNA (mtDNA). Analysis suggests that the differentiation within this species could be related to quaternary glaciation. We also compared the clustering analysis results generated using mtDNA versus microsatellite markers. The most obvious difference between the two markers was reflected in the fact that the populations Yelashan (YLS) and Bangda (BD) represented a separate cluster based on the microsatellite data. Additionally, our studies involving the comparison of the two markers also reveal that a few of the differentiated populations from the middle region have a higher admixture. Due to a faster mutation rate, microsatellite markers offer better insight into newer evolutionary events than mtDNA, likely resulting in the discrepancies we observe between the two markers. We suggest that combining the two different markers may be beneficial in obtaining a more complete assessment of the population genetic structure.

Article Information

Received 25 November 2019

Revised 11 March 2020

Accepted 20 May 2020

Available online 09 February 2021

(early access)

Published 28 December 2021

Authors’ Contribution

HY, SJ and DY involved in planning the study design, laboratory analyses, analyzing the data, and writing the manuscript. XH contributed to the laboratory work. LG contributed to the sampling and analyzing the data. ZT contributed to the analysis of data.

Key words

Ochotona curzoniae, mtDNA, Microsatellite, Population structure, Glaciation

DOI: https://dx.doi.org/10.17582/journal.pjz/20191125131109

* Corresponding author: duyurong123@outlook.com

0030-9923/2022/0002-0503 $ 9.00/0

Copyright 2022 Zoological Society of Pakistan

INTRODUCTION

The Qinghai-Tibet Plateau (QTP) is the largest and highest plateau in the world with an average elevation exceeding 3000m and covering most of the Tibet and Qinghai provinces (Zhang et al., 2002). Since the Pliocene, the plateau has experienced severe uplift events. The altitude of the plateau exceeded 3000m after the Kunlun-Yellow River Movement, causing the plateau to enter the cryosphere (Li and Fang, 1998). Combined with global climate change, the climate of the plateau entered a glacial fluctuation in the Quaternary, which led to the expansion and constriction of fauna living on the plateau (Liu, 1999; Jin and Liu, 2010; Hofmann, 2012). The uplift of the plateau also resulted in complex geographic and topographic features, such as a series of large mountains stretching from the southern to the northern plateau, valleys, basins, rivers and lakes (Pan et al., 2004). These directly affected the distribution and dispersal of the species living on the plateau; as a result, such factors have had a major impact on the genetics of these populations and the resultant phylogeographical patterns that we observe today (Tang and Li, 2001; Chakraborty et al., 2015).

The plateau pika (Ochotona curzoniae) is a diurnal mammal that is widely distributed across the QTP, and is viewed as a keystone species of the plateau (Lai and Smith, 2003). The plateau pika lives in the alpine meadow and alpine steppe at altitudes of 3000–5100m (Smith and Forrin, 1999). About 1.2 million years ago (Mya), the plateau pika separated from its closest living relative, Ochotona nubrica (Yu et al., 2000). We hypothesize that changes in climate and geography of the QTP during the Pleistocene influenced the subsequent migration and dispersal of this species. Liu et al. (2013) had previously analyzed the population differentiation of plateau pika using 12s rRNA and Cyt b. The results revealed that the populations located in the south of Brahmaputra River were clustered together and the populations in the north of Qaidam Basin and Gonghe Basin were clustered together. The authors suggested that the geographical barriers played the most important role in mediating the population differentiation of this species (Liu et al., 2013). Similar results were also reported by Chang (2016) using Cyt b of Ochotona curzoniae, however, they didn’t interpret the observations. We have compiled data on the 12s rRNA, Cyt b, tRNA-Thr/Pro and partial D-loop, and reconstructed the population history dynamics. We found that the effective population size of plateau pika is affected by the Quaternary climatic and geological changes (He et al., 2018). And we expect to find evidence that historical and current ecological conditions have shaped the population genetic structure of the plateau pika in this region.

Interestingly, different genetic markers have varying capabilities of detecting genetic variation and polymorphisms. Phenomena such as independent evolution of mtDNA and nuclear genes could complicate population genetic analyses. In this context, investigations into the population genetic structure of a species often require the use of multiple markers to minimize the effects of the above-mentioned marker-specific complications (Flanders et al., 2009).

Mitochondrial DNA is commonly used in the study of population structure because of its high rates of substitution and maternal inheritance. However, mtDNA only represents a single locus and is an extra-nuclear genetic marker with evolutionary dynamics that can be divergent from nuclear DNA (Flanders et al., 2009). Furthermore, due to the fact that nuclear and mtDNA evolve at different rates, patterns of variation among mtDNA often do not mirror those among nuclear genes (Naidoo et al., 2016). To date, there is no report that has studied the population genetic structure of plateau pika on the basis of microsatellite data. In this study, we have sequenced microsatellite markers of plateau pika from the same locations in the QTP spanning Qinghai, Tibet, Xinjiang, and Sichuan where mtDNA data has already been published and assessed the population genetic structure using both mtDNA and microsatellite markers. We also investigated the evolutionary determinants of genetic structure in the plateau pika, and why we observed discordant patterns from mitochondrial and nuclear genetic data.

MATERIALS AND METHODS

Sample collection and DNA extraction

A total of 226 O. curzoniae samples were collected from 34 locations throughout the species distribution area described in Figure 1 and Table I. The same samples as our previous research (He et al., 2018) were used in this study. The muscle samples obtained from pika were preserved in 95% ethanol. Voucher specimens were kept in the Northwest Institute of Plateau Biology, Chinese Academy of Sciences. Total genomic DNA was extracted from ~15mg of muscle tissue per animal using the TIANamp Genomic DNA Kit (Tiangen Inc., Beijing, China).

 

Table I. Details of the sampling sites.

Sampling site code	Sampling site	Latitude (°N)	Longitude (°E)	Altitude (m)	Sample size
YSL	Yelashan	30.187	97.295	4338	14
BD	Bangda	30.529	97.128	4348	10
JZ	Jiangzi	28.901	90.101	4660	10
JG	Jungong	34.647	100.592	3467	6
HB	Haibei	37.685	101.274	3238	6
NQ	Naqu	31.441	92.277	4420	6
QL	Qilian	37.959	100.255	3645	6
GC	Gangcha	37.503	100.482	3164	6
EL	Eling lake	35.074	97.716	4234	8
T1	Tibet1	29.237	87.218	4481	5
T3	Tibet3	30.578	82.536	4944	8
T4	Tibet4	29.493	85.089	4607	4
TJ	Tianjun	37.174	99.283	3275	6
CD	Chengduo	33.356	97.236	4348	8
ZD	Zhiduo	33.54	96.064	4308	7
KLS	Kunlunshan	35.734	94.310	4110	7
SQ	Shiqu	33.033	98.015	4254	8
TT	Tuotuo river	34.330	92.591	4533	5
NM	Nimu	29.502	90.270	3908	10
AT	Tuzi lake	36.800	87.308	4750	8
AQ	Aqike lake	37.003	88.610	4250	6
AA	Kaerqiuer	37.043	90.755	4184	6
AD	Anduodong	32.309	91.732	4675	5
AX	Anduoxi	32.178	91.593	4586	3
NMC	NamCo	30.721	91.035	4807	6
XDT	Xidatan	35.712	94.058	4532	8
XX	Xingxing sea	34.833	98.135	4194	8
HSX	Huashixia	35.082	98.855	4237	7
BL	Beilu river	34.862	92.942	4590	5
CM	Chumaer river	35.356	93.386	4520	2
LKZ	Langkazi	29.109	90.417	4425	8
GL	Geladandong	33.589	91.652	4873	4
QML	Qumalai	34.142	95.884	4389	7
QHL	Qinghai lake	36.630	100.107	3175	5

 

MtDNA sequencing and microsatellite genotyping

The following mtDNA sequences were analyzed: Cytb (1140bp), tRNA-Thr/Pro (135bp), D-loop (342bp), and 12SrRNA (946bp). All the mtDNA sequences used in this study were derived from previous studies conducted in our laboratory and the GenBank accession numbers of the haplotypes of these sequences are FJ227330–FJ227482, KF038167–KF038213, and NM225702–NM225732 (Ci et al., 2009; Liu et al., 2013; He et al., 2018). Each individual’s sequences were concatenated for further analysis.

The total genomic DNA was used to amplify the nuclear genes. Ten published microsatellite loci (P7, P47, P63, P120, P124, P149, P172, P175, OCP7, OCP8) were analyzed in all samples. Primers for OCP7 and OCP8 were developed from Ochotona princeps sequences (Peacock et al., 2002); and other primers (P7, P47, P63, P120, P124, P149, P172, and P175) were developed from O. curzoniae (Li et al., 2009). Polymerase chain reaction (PCR) amplification was conducted in a total reaction volume of 50 µL containing 5.0 µL 10X Taq buffer (TaKaRa, Dalian, China), 2.5U Taq DNA polymerase (TaKaRa), 1.5mM MgCl2, 0.1mM dNTPs, 0.4µM each primer, 40–60 ng genomic DNA, and sufficient ddH2O to achieve the final volume. PCR was performed as follows: initial denaturation at 95°C for 8 min; followed by 35 cycles of 95°C for 45 s, 52–60°C for 45 s, and 72°C for 80 s; and a final extension at 72°C for 7 min. Fragment length was analyzed using an ABI 3700 sequencer (Applied Biosystems, Carlsbad, CA, USA). Allele calling and sizing was conducted using Genescan v3.1.2.

Genetic analysis

For mtDNA data, the number of haplotypes (M), haplotype diversity (h), and nucleotide diversity (π) were estimated using DnaSPv5 (Rozas et al., 2003). Molecular diversity indices of the microsatellite data, including the number of alleles (Na), effective number of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (FIS), were calculated using the POPGENE software package (Yeh et al., 1999). Populations experiencing heterozygous deficiency were detected using GENEPOP v4.6 through the detection of significant departures from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (Rousset, 2008). Parameters of the Markov Chain employed in estimations of HWE and linkage disequilibrium were 10,000 dememorizations, 1000 batches, and 10,000 iterations per batch.

Estimates of divergence times

To infer the phylogenetic relationships and coalescence date of the mtDNA haplotypes, we used a Bayesian Markov chain Monte Carlo (MCMC) in BEAST v1.8.0 (Drummond and Rambaut, 2007). Because there is no suitable fossil record available for calibrating the mutation rate in the plateau pika, we adopted a substitution rate of 4.8% per million years based on the average mutation rate of the Cytb/D-loop of the American pika (Ochotona princeps) (Drummond and Rambaut, 2007). The best-fitting model of sequence evolution was selected using jModelTest 2.1.4 (Galbreath et al., 2009). We used the HKY+G+I model of nucleotide substitution and a strict clock model with the default parameters and default operators set in BEAUti to estimate the divergence times. Each MCMC sample was based on a run of 100,000,000 generations, and sampled every 1,000 generations. All samples were then examined in Tracer v1.6 (Rambaut and Drummond, 2013) with the first 10% discarded as burn-in to produce the final population trajectory. All runs were combined in LogCombiner 1.4.8 (Drummond et al., 2012), and the resulting tree was visualized using Figtree 1.2 (Rambaut, 2008).

Population structure

Population cluster analysis was conducted with mtDNA data using the spatial Bayesian clustering method in BAPS v6.0 (Cordander et al., 2004). The mixture analysis was conducted using the “spatial clustering of groups” method. The simulation was run from K=2 to K=34. The optimum number of genetic clusters was determined based on the partition with the maximum likelihood and highest posterior probability. The index of differentiation (FST) was calculated using mtDNA data in Arlequin v3.5, and this estimated the degree of divergence between groups based on the result of population cluster analysis (Excoffier and Lischer, 2010).

We further analyzed the genetic structure of the plateau pika based on the microsatellite dataset using a Bayesian model-based clustering method in STRUCTURE v2.3 (Pritchard et al., 2000). The most likely number of genetic clusters (K) was estimated independently in ten independent runs of K= 2-34, based on the admixture model with correlated allele frequencies. The length of the MCMCs was 1,000,000 steps after a burn-in period of 100,000 steps. The most likely value of K was estimated with Structure Harvester (Earl, 2012) using the statistic ∆K (Evanno et al., 2005). Individual assignment probability plots were generated using CLUMPP v1.1.2 (Jakobsson and Rosenberg, 2007) and DISTRUCT v1.1 (Rosengerg, 2004).

Analysis of molecular variance (AMOVA) in Arlequin v3.5 under alternative hierarchical arrangements were conducted for both mtDNA and microsatellite data. The significance of the between-group variance was assessed with 10,000 random permutations.

RESULTS

Genetic diversity at mtDNA and microsatellite loci

The mtDNA sequences for Cytb, tRNA, 12SrRNA, and part of the D-loop were combined for each individual to analyze genetic diversity. Among the 226 samples, we identified a total of 139 haplotypes containing 356 polymorphic sites. Most haplotypes appeared only once, resulting in relatively high haplotype diversity (0.991). Nucleotide diversity for all samples was 0.0184.

A total of 135 alleles were detected among the ten microsatellites that were analyzed across the 34 populations (Table II). The number of alleles ranged from seven at locus P120 to 29 at locus P7. There was no linkage disequilibrium among microsatellite loci after standard Bonferroni correction for multiple comparisons (Rice, 1989). Ho ranged from 0.532 to 0.827 and He ranged from 0.750 to 0.940, and FIS was negative at P7, P47, P63, P120, P124, P149, P175, and OCP7, indicating an abundance of heterozygotes at these loci.

Divergence times

Coalescence analysis of the mtDNA sequences revealed four well-supported lineages (Fig. 2). Molecular estimates of divergence times suggested that the southern clade of O.curzoniae diverged approximately 0.76 (95% highest posterior density (HPD): 0.62–0.91) Mya, followed by the divergence of the eastern clade approximately 0.33 (95% HPD: 0.27–0.40) Mya. Finally, the divergence of the northern and western clades occurred about 0.29 (95% HPD: 0.23–0.36) Mya.

Population structure

Next, we analyzed the population structure of the plateau pika at the QTP. Clustering analysis using the mtDNA data separated the 34 localities into four distinct groups (Fig. 3), which were in agreement with the phylogenetic analysis. The southern clade consisted of the JZ and LKZ populations, while the northern clade consisted of the QHL, TJ, GC, HB, and QL populations. The remaining populations were separated into eastern and western clades following their geographic locations. The pairwise FST estimate for mtDNA indicated that the four clades were highly differentiated, and the lowest FST value was between the eastern and western clades (Table III).

 

Table II. Genetic variation among microsatellite loci.

Marker	Na	Ne	Ho	He	FIS
OCP8	27	9.848	0.540	0.900	0.197
P172	12	4.016	0.532	0.752	0.034
P175	9	4.313	0.628	0.768	-0.075
P120	7	3.047	0.621	0.669	-0.231
P47	16	4.000	0.649	0.750	-0.066
OCP7	22	9.134	0.750	0.893	-0.118
P124	16	9.240	0.782	0.894	-0.079
P149	20	7.202	0.783	0.863	-0.155
P63	15	9.333	0.786	0.895	-0.122
P7	29	16.149	0.827	0.940	-0.058
Mean	17.3	7.628	0.690	0.832	-0.067

 

Na, Number of alleles; Ne, Effective number of alleles; Ho, Observed heterozygosity; He, Expected heterozygosity; FIS, Inbreeding coefficient.

 

Table III. FST between groups using mtDNA.

	I	II	III	IV
I	-	-	-	-
II	0.81305	-	-	-
III	0.73467	0.54529	-	-
IV	0.76557	0.59331	0.46787	-

 

I, southern clade; II, northern clade; III, eastern clade; IV, western clade.

In contrast, Bayesian analysis based on the microsatellite data supposed a model with the highest ΔK value when K=3 (Fig. 4). Populations formed into three clusters, with YLS and BD consisting of one cluster, and the other populations separated into two clusters (K=3, Fig. 5). However, increasing the K-value provided additional information. When K=4, populations HB, QL, GC, TJ, and QHL grouped together, where QHL had higher admixture. When K=5, the JZ and LKZ populations separated into a clusters, and LKZ had higher admixture (Figs. 5 and 6).

 

We conducted an AMOVA analysis to test the variance of the mtDNA data, and found that 58.18% of the variance resulted from differences among groups, and 41.82% were within groups. In contrast, AMOVA analysis of the microsatellite data showed 13.52% and 86.48% of the variance among plateau pika groups and within populations, respectively. However, we have detected significant genetic variance at all three hierarchical levels tested (i.e. among regions, among populations, and within populations; P< 0.001) for both types of marker (Table IV).

 

Table IV. Hierarchical analysis of molecular variance (AMOVA).

Source of variation	Variation (%)	Fixation indices	P-value
Mitochondrial DNA
Among groups	58.18	FCT=0.58184	<0.001
Among populations	12.4	FSC=0.29651	<0.001
Within populations	29.42	FST=0.70583	<0.001
Microsatellites
Among populations	13.52	FSC=0.08782	<0.001
Within populations	86.48	FST=0.10321	<0.001

 

 

 

 

DISCUSSION

In this study, we used two different markers to investigate the population genetic structure of plateau pika in the QTP. We found thirty-four populations of plateau pika were separated into four highly differentiated clades based on the mtDNA data, consistent with previous studies (Liu et al., 2013). In this research, we estimated the divergence times of the four clades, and found that the separation of the southern clade likely occurred about 0.76 Mya, which corresponds to the Kunlun glaciation (Shi, 2002). Geological evidence suggests that this was the largest glaciation event in the QTP, and the ice sheet covered an area 18 times larger than that of the current glacier (Li and Pan, 2002; Cui et al., 2011). The glaciation may have constrained the pika populations to the south of Brahmaputra River and caused differentiation. Our analysis estimated that the separation of the eastern clade as well as the differentiation of the western and northern clades occurred in the middle of the Pleistocene, coinciding with the Guxiang glaciation of the QTP (Cui et al., 2011). The central populations would have retreated to marginal shelters during that period, leading to differentiation among clades.

The most obvious difference between the results of population structure analysis using different markers was that the YLS and BD populations represented a separate cluster based on the microsatellite data. On the contrary, the other populations analyzed on the basis of the microsatellite data, were separated into two clusters, when K=3 in Bayesian clustering analysis. With the increasing K value, the separation trend of the populations excluding YLS and BD, was similar between the microsatellite and mtDNA. However, there is higher admixture in a few populations which are located at the transition zone of different clades when K=5 (Fig. 5).

The sampling locations of populations YLS and BD are between the Lantsang River and Tenasserim chain, and this region is characterized by the longitudinal range-gorge formed in a parallel arrangement between the mountains and rivers after the Gonghe Movement (Ming, 2006). This complex geography may underlie the differentiation of the YLS-BD cluster. We did not identify this cluster using mtDNA, and this may be due to the higher mutation rates in the microsatellites compared to mtDNA (Kato et al., 2015), making microsatellites more informative for recent evolutionary events. The formation of the longitudinal range-gorge likely occurred too recently for the accumulation of genetic differences in mtDNA of the pika populations.

In the mtDNA analysis, populations JZ and LKZ formed a clade; interestingly, these populations are separated from other populations by the Brahmaputra River. However, it is generally recognized that the Brahmaputra River formed before the speciation of plateau pika (Zhu, 2012). Liu et al. (2013) suggested the existence of several beaded river valleys along the Brahmaputra River. These rivers were prone to clogging, thereby, leading to the formation of ancient dammed lakes. These natural dams formed due to temporary land formations, may have provided a channel for the dispersal of the pikas during this period. The microsatellite analysis suggested that the LKZ population was undergoing admixture with the western group, which further supports the hypothesis that there were channels between the two sides of the Brahmaputra River. In contrast, the JZ population samples were collected from a location surrounded by mountains, likely resulting in limited gene flow with other populations.

When assessing population structure using the microsatellite data, populations HB, GC, TJ, QL, and QHL were differentiated as a cluster when K was set to 4 in the STRUCTURE analysis. Although the temperature and precipitation were higher during the interglacial period, the climates of the Qiadam Basin and Gonghe Basin have been dry and cold since the last glacial period, gradually evolving into the current arid habitat (Shi et al., 2004; Sun et al., 2010). The Qiadam and Gonghe Basin does not provide suitable habitat for the plateau pika, and these basins separate pika populations from each other. The QHL population showed higher admixture, which could be due to the recent aridification of Gonghe Basin that resulted in relatively recent separation of this population (Sun et al., 2007). The QHL population was sampled from the south side of Qinghai Lake, while the other four populations in this area were sampled from the north side of the lake. Qinghai Lake thus represents a geographic barrier that enhances the isolation of pika populations in the QTP.

We also found separation between the eastern and western clusters using microsatellite analysis. This continued separation could be caused by the development of vast pan-lakes in the plateau during the interglacial period following the Guxiang glaciation. The Tuotuohe pan-Lake, a northern Tibet pan-lake in the central plateau, likely reduced migration between the eastern and western populations (Wu et al., 2009). Although the vast pan-lakes shrank 45,000 years ago, plenty of rivers and lakes in this area are retained (Yu, 2008), further restricting gene flow. The populations XX, ZD, XDT, KLS, and NQ (K=3) had higher admixture in the microsatellite analysis. These populations were from the middle region of the differentiated populations, and this is likely caused by the population expansion during the interglacial periods and after the retreat of vast pan-lakes.

In our study, we have suggested that the analysis based on the mtDNA provides information about the differentiation of plateau pika in the more distant past, while microsatellite data could reflect more recent evolutionary events. Combining these two different markers, each having different evolutionary patterns, could better expose the population genetic structure of plateau pika and reveal underlying reasons for the evolution of the population genetics. We also found the differentiation times of plateau pika populations coincide with the ancient climate and geology events of Qinghai-Tibet Plateau. Hence, this species not only plays important roles in the plateau ecology, but also could be a good model to study the response to climate and geological changes of the plateau biota.

CONCLUSION

Our data indicates discrepancies between the population genetic structures of plateau pika as inferred from their mtDNA versus their microsatellite markers. Further, exploring the mtDNA revealed that the population divergence time of plateau pika is related to the Quaternary Glaciation. On the other hand, results from the analysis of microsatellites unravel the principal cause of persistent differentiation as well as new differentiation events of different populations of this species. By combining the results obtained when analyzing two different markers, we put forth an extensive history of the differentiation of plateau pika population in the QTP. In conclusion, we propose that it is necessary to use multiple different markers to carry out research on population genetic structure. In the future, we would like to verify our results using more nuclear genes.

ACKNOWLEDGEMENT

This work was supported by the Foundation of Qinghai Science & Technology Department (grant number 2019-ZJ-986Q).

Statement of conflicts of interest

The authors have declared no conflict of interest.

REFERENCES

Chakraborty, D., Ramakrishnan, U. and Sinha, A., 2015. Quaternary climate change and social behavior shaped the genetic differentiation of an endangered montane primate from the southern edge of the Tibetan Plateau. Am. J. Primatol. 77: 271-284. https://doi.org/10.1002/ajp.22343

Chang, Y.B., 2016. Comparative phylogeography of Ochotona curzoniae and O. cansus. Master’s thesis. Hebei University, Shijiazhuang, China.

Ci, H.X., Lin, G.H., Cai, Z.Y., Tang, L.Z., Su, J.P. and Liu, J.Q., 2009. Population history of the plateau pika endemic to the Qinghai-Tibetan Plateau based on mtDNA sequence data. J. Zool., 279: 396-403. https://doi.org/10.1111/j.1469-7998.2009.00635.x

Corander, J., Waldmann, P., Marttinen, P. and Sillanpää, M.J., 2004. BAPS 2: Enhanced possibilities for the analysis of genetic population structure. Bioinformatics, 20: 2363-2369. https://doi.org/10.1093/bioinformatics/bth250

Cui, Z.J., Chen, Y.X., Zhang, W., Zhou, S.Z., Zhou, L.P., Zhang, M. and Li, C.C., 2011. Research history, glacial chronology and origins of Quaternary glaciations in China. Quat. Sci., 31: 749-764.

Darriba, D., Taboada, G.L., Doallo, R. and Posada, D., 2012. jModelTest 2: More models, new heuristics and parallel computing. Nat. Methods, 9: 772. https://doi.org/10.1038/nmeth.2109

Drummond, A.J. and Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol., 7: 214. https://doi.org/10.1186/1471-2148-7-214

Drummond, A.J., Suchard, M.A., Xie, D. and Rambaut, A., 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol., 29: 1969–1973. https://doi.org/10.1093/molbev/mss075

Earl, D.A., 2012. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour., 4: 359. https://doi.org/10.1007/s12686-011-9548-7

Evanno, G., Regnaut, S. and Goudet, J., 2005. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 14: 2611. https://doi.org/10.1111/j.1365-294X.2005.02553.x

Excoffier, L. and Lischer, H., 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour., 10: 564–567. https://doi.org/10.1111/j.1755-0998.2010.02847.x

Flanders, J., Jones, G., Benda, P., Dietz, C. and Zhang, S., 2009. Phylogeography of the greater horseshoe bat, Rhinolophus ferrumequinum: Contrasting results from mitochondrial and microsatellite data. Mol. Ecol., 18: 306-318. https://doi.org/10.1111/j.1365-294X.2008.04021.x

Galbreath, K.E., Hafner, D.J. and Zamudio, K.R., 2009. When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotona princeps). Evolution, 63: 2848-2863. https://doi.org/10.1111/j.1558-5646.2009.00803.x

He, Y.J., Lin, G.H., Ci, H.X., Liu, C.X., Zhang, T.Z. and Su, J.P., 2018. The past population dynamics of Ochotona curzoniae and the response to the climate change. North-West J. Zool., E171704.

Hofmann, S., 2012. Population genetic structure and geographic differentiation in the hot spring snake Thermophis baileyi (Serpentes, Colubridae): Indications for glacial refuges in southern-central Tibet. Mol. Phylogenet. Evol., 63: 396-406. https://doi.org/10.1016/j.ympev.2012.01.014

Jakobsson, M. and Rosenberg, N.A., 2007. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23: 1801. https://doi.org/10.1093/bioinformatics/btm233

Jin, Y.T. and Liu, N.F., 2010. Phylogeography of Phrynocephalus erythrurus from the Qiangtang Plateau of the Tibetan Plateau. Mol. Phylogenet. Evol., 54: 933–940. https://doi.org/10.1016/j.ympev.2009.11.003

Kato, A.B., Hyseni, C., Okedi, L.M., Ouma, J.O. and Aksoy, S., 2015. Mitochondrial DNA sequence divergence and diversity of Glossina fuscipes fuscipes in the Lake Victoria basin of Uganda: implications for control. Paras. Vector, 8: 1-13. https://doi.org/10.1186/s13071-015-0984-1

Lai, C.H. and Smith, A.T., 2003. Keystone status of plateau pikas (Ochotona curzoniae): Effect of control on biodiversity of native birds. Biodivers. Conserv. 12: 1901–1912.

Li, B.Y. and Pan, B.T., 2002. Progress in paleogeographic study of the Tibetan plateau. Geogr. Res., 21: 61-70.

Li, J.J. and Fang, X.M., 1998. The research of uplift and climate change of Qinghai-Tibetan Plateau, Chinese. Sci. Bull., 43: 1569-1574.

Li, K.X., Geng, J.N. and Yang, J., 2009. Isolation and characterization of 13 microsatellite loci in the plateau pika (Ochotona curzoniae). Conserv. Genet., 10: 785-787. https://doi.org/10.1007/s10592-008-9662-6

Liu, C.X., Su, J.P., Zhang, T.Z. and Lin, G.H., 2013. The effect of Qingzang-Tibet Plateau geographical barrier on plateau pika population differentiation. Sichuan J. Zool., 32: 651-657.

Liu, X.D., 1999. Influences of Qinghai-Xizang (Tibet) Plateau uplift on the atmospheric circulation, global climate and enviroment changes. Plt. Meteorol., 18: 321-332.

Ming, Q.Z., 2006. The landform development and environment effects of three parallel rivers in the north longitudinal range-gorge region (LRGR). PhD thesis, Lanzhou University, Gansu, China.

Naidoo, T., Schoeman, M.C., Goodman, S.M. and Taylor, P.J., 2016. Discordance between mitochondrial and nuclear genetic structure in the bat Chaerephon pumilus (Chiroptera: Molossidae) from southern Africa. Mammal. Biol., 81: 115-122. https://doi.org/10.1016/j.mambio.2015.11.002

Pan, B.T., Gao, H.S., Li, B.Y. and Li, J.J., 2004. Step-like landforms and uplift of the Qinghai–Xizang Plateau. Quat. Sci., 24: 50-56.

Peacock, M.M., Kirchoff, V.S. and Merideth, S.J., 2002. Identification and characterization of nine polymorphic microsatellite loci in the North American pika, Ochotona princeps. Mol. Ecol. Notes, 2: 360-362. https://doi.org/10.1046/j.1471-8286.2002.00249.x

Pritchard, J.K., Stephens, M. and Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics, 155: 945.

Rambaut, A., 2008. FigTree version 1.2. Computer program available from website. http://tree.bio.ed.ac.uk/software/figtree.

Rambaut, A. and Drummond, A.J., 2013. Tracer v1.6. http://tree.bio.ed.ac.uk/software/ tracer/

Rice, W.R., 1989. Analyzing tables of statistical tests. Evolution, 43: 223-225. https://doi.org/10.1111/j.1558-5646.1989.tb04220.x

Rosenberg, N.A., 2004. DISTRUCT: A program for the graphical display of population structure. Mol. Ecol. Notes, 4: 137. https://doi.org/10.1046/j.1471-8286.2003.00566.x

Rousset, F., 2008. GENEPOP’ 007: A complete re-implementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Res. 8: 103–106. https://doi.org/10.1111/j.1471-8286.2007.01931.x

Rozas, J., Sanchez-DelBarrio, J.C., Messeguer, X. and Rozas, R., 2003. DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics, 19: 2496-2497. https://doi.org/10.1093/bioinformatics/btg359

Shi, W., Ma, Y.S., Wu, M.L., Du, J.J. and Zhang, X.J., 2004. Quaternary sporopollen assemblages and environmental evolution of the Gonghe Basin on the northeastern margin of the Qinghai-Tibet plateau. J. Geomech., 10: 310-318.

Shi, Y.F., 2002. A suggestion to improve the chronology of Quaternary glaciations in China. J. Glaciol. Geocryol., 24: 687-692.

Smith, A.T. and Foggin, J.M., 1999. The plateau pika (Ochotona curzoniae) is a keystone species forbiodiversity on the Tibetan plateau. Anim. Conserv., 2: 235-240. https://doi.org/10.1111/j.1469-1795.1999.tb00069.x

Sun, F.F., Zhang, W.Y., Gong, J.C. and Zhang, C.J., 2010. The palaeoenvironmental reconstruction on pollen proxy in the Qaidam Basin since late Pliocene. Geol. Rev., 56: 621-628.

Sun, Y.G., Fang H.B., Zhang K., Zhao, F.Y. and Liu S.Y., 2007. Step-like landform system of the Gonghe basin and the uplift of Qinghai-Tibet plateau and development of the Yellow River. Geol. China, 34: 1141-1147.

Tang, L.Y. and Li, C.H., 2001. Temporal-spatial distribution of the Holocene vegetation in the Tibetan Plateau. J. Glaciol. Geocryol., 23: 367-374.

Wu, Z.H., Wu, Z.H., Hu, D.G., Zhou, C.J., Ye, P.S. and Zhang, Y.L., 2009. Vast paleo-lakes, planation surface and topographic evolution of the Tibetan plateau. Geoscience, 23: 993-1002.

Yeh, F.C., Yang, R.C., Boyle, T., Ye, Z.H. and Mao, J., 1999. POPGENE version 1.32, the user-friendly shareware for population genetic analysis. Molecular biology and biotechnology centre, Canada, University of Alberta, http://www. ualberta.ca/fyeh/.

Yu, J., 2008. The Quaternary environmental evolution of typical lakes in central Tibetan Plateau. Master’s thesis. Chinese Academy of Geological Sciences, Beijing, China.

Yu, N., Zheng, C., Zhang, Y.P. and Li W.H., 2000. Molecular Systematic of pikas inferred from mitochondrial DNA sequences. Mol. Phylogenet. Evol., 16: 85-95. https://doi.org/10.1006/mpev.2000.0776

Zhang, Y.L., Li, B.Y. and Zheng, D., 2002. A discussion on the boundary and area of the Tibetan Plateau in China. Geogr. Res., 21: 1-8.

Zhu, S., 2012. River landform and geology environment evolution in the Yarlung Zangbo River vally. Master’s thesis. Academy of Geological Science, Beijing, China.