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 Article Information Received 25 November 2019 Revised 11 March 2020 Accepted 20 May 2020 Available online 09 February 2021


INTRODUCTION
T he 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., O n l i n e F i r s t A r t i c l e 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.

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

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.

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 (F IS ), 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 (F ST ) 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 Population Structure of Plateau Pika O n l i n e

F i r s t A r t i c l e
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.

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 F IS was negative at P7, P47, P63, P120, P124, P149, P175, and OCP7, indicating an abundance of heterozygotes at these loci.

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 F ST estimate for mtDNA indicated that the four clades were highly differentiated, and the lowest F ST value was between the eastern and western clades (Table III).  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).

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 O n l i n e

F i r s t A r t i c l e
Population Structure of Plateau Pika 7 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 panlakes 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.

O n l i n e F i r s t A r t i c l e O n l i n e F i r s t A r t i c l e
Population Structure of Plateau Pika