Phylogeography of the Tibetan hamster Cricetulus kamensis in response to uplift and environmental change in the Qinghai–Tibet Plateau

Abstract Aim The evolutionary process of an organism provides valuable data toward an understanding of the Earth evolution history. To investigate the relationship between the uplift of the Qinghai–Tibet Plateau (QTP) and mammalian evolution since the late Cenozoic, the geographic distribution of genetic variations in the Tibetan hamster Cricetulus kamensis was investigated using phylogeographical methods. In particular, population divergence, demographic history, genetic variation, and the prediction of species distribution area were investigated. Location The Qinghai–Tibet Plateau. Methods A total of 53 specimens, representing 13 geographic populations, were collected from the QTP. The phylogeographical pattern and demographic history of C. kamensis were analyzed, and the probable factors in the QTP uplift and the Quaternary glacial periods were inferred from one nuclear and four mitochondrial genes. Furthermore, the species distribution model (SDM) was used to predict changes in potentially suitable habitats since the last Interglacial. Results Phylogenetic analysis demonstrated that two major genetic differentiations of the C. kamensis population occurred during the Early Pleistocene that were influenced by the Qing‐Zang tectonic movement from the Middle Pliocene to the Early Pleistocene. Genetic distance between two major clades indicated low genetic divergence. Demographic history analysis showed that the C. kamensis population was affected by the Quaternary glacial period. SDM analysis indicated that C. kamensis was endemic to the QTP and the suitable habitat was affected by climate change, especially during the Last Glacial Maximum (LGM). Main conclusion Our results indicated that the QTP uplift led to the population divergence of C. kamensis, and vicariance well accounted for the geographic distribution of genetic variation in C. kamensis as a result of genetic divergence and lack of gene flow. The genetic distance shows that C. alticola may be a subspecies of C. kamensis. Demographic history analysis suggests that the QTP was affected by the last glacial period. SDM analysis supports that almost the entire QTP is covered by a huge ice sheet during the LGM.


| INTRODUC TI ON
Phylogeography is an integrative discipline that focuses on the principles and processes underlying the geographic distributions of genealogical lineages, namely those that clarify geological history and phylogenetic components of the spatial distributions of gene lineages (Avise, 2000). Molecular phylogenetic reconstructions of biotas can help to test biogeographic hypotheses, illustrate evolutionary relationships between spatial and temporal population structures of species, and support speculations of speciation mechanisms using genetic data across the distribution range of species under the phylogeographical approach. The Qinghai-Tibet Plateau (QTP) is a relatively young geological structure with a complicated evolutionary history, complex landform configuration, unique climatic environment, and diversified habitats (Favre et al., 2015;Lei, Qu, & Song, 2014). During recent years and by using a phylogeographical approach, a large number of studies indicated that the distribution patterns and population genetic structures of plateau-dwelling organisms are related to the tectonic configuration and climate change in the QTP (Deng, Wang, Wang, Li, & Hou 2015;Favre et al., 2015;Lei et al., 2014;Yang, Dong, & Lei, 2009).
Historically, the uplift of the QTP happened in several steps and followed an asynchronous uplift with several partitioned blocks, showing a progressive uplift from south to north (Mulch & Chamberlain, 2006;Wang et al., 2008a). However, the rapid uplift during the late Cenozoic plays a major role in shaping the landform and climate of the QTP, as well as exerting a remarkable effect on the biota Shi, Li, & Li., 1998). Its process of geographic formation and the evolution of biotas remains insufficiently understood despite its outstanding geographic extent so far (Favre et al., 2015). Several recent studies have suggested that the uplifting of the QTP plays a major role in shaping plateau biodiversity due to strong ecological gradients in temperature, precipitation, and topography which contribute to the generation of a variety of terrestrial habitats Favre et al., 2015;Lei et al., 2014). A large number of studies demonstrated that the phylogeny, speciation, and demography of organisms are associated with the uplift of the QTP and simultaneous changes of climate, for example in frogs, lizards, birds, and mammals (James et al., 2003;Jin, Brown, & Liu, 2008;Liu et al., 2015;Tseng et al., 2014). With respect to mammals, Yu et al. (1992) proposed that the diversification of pika correlates with the uplift of the QTP in China; furthermore, geographic isolation created numerous diversified microhabitats for pika survival. Molecular systematics of the genus Ochotona also indicate that the speciation of pikas is consistent with the historical episodes of geologic and climatic changes of the QTP (Niu, Wei, Ming, Liu, & Feng, 2004;Yu, Zheng, Zhang & Li, 2000); however, studies on the intraspecific divergence of pika on the QTP have not been reported. The population divergence (between 2.4 and 4.4 Mya) of the Tibetan gazelle (Procapra picticaudata) was influenced by the uplift of the QTP (Zhang & Jiang, 2006), whereas its dispersal capacity far exceeds that of muroid rodents. In addition, a similar result was reported in a study of the plateau zokor (Eospalax baileyi) (Tang et al., 2010). Its life history is completely different from the extant mice. The hamster is an ancient animal in comparison to the aforementioned mammals, and their fossil data can be traced back to the middle Eocene (approximately 40 Mya). The hamster has been considered to originate from Asia (Huang, 2004;Storer, 1988;Tong, 1992). Herein, the phylogeographical study of the hamster living on the QTP is more likely to reveal the geographic and climatic evolution of the QTP. Small mammals, such as rodents, are generally good models to test the relationships between geology and climatic events as well as for the study of the evolution and adaption of animals, such as their localized habitat requirements and low dispersal abilities. These animals are generally limited to persist on a particular site even under extreme environmental conditions (Nava-García, Guerrero-Enríquez, & Arellano, 2016;Wang, Zhao, Fang, Liao, & Liu, 2014). Given their tendency toward a short generation time and rapid mitochondrial DNA substitution, speciation seems to be their main strategy to adapt to the original environment despite habitat constraints (Ben Faleh et al., 2013).
The Tibetan hamster Cricetulus kamensis (Satunin, 1902) is endemic to the QTP and is a small to medium sized hamster, similar to C. longicaudatus with restricted distribution in the QTP (Feng, Cai, & Zheng, 1986;Luo et al., 2000). To date, there have been relatively few publications about C. kamensis, which might be the reason why this species is poorly known outside of China. The classification of C. kamensis has long been controversial. Historically, Tibetan hamster has been classified into four species (U. kamensis, C. kozlovi, C. kama, and C. alticola) from different biogeographic populations based on their respective morphological characteristics (Bonhote, 1905;Satunin, 1902;Thomas, 1917). As such, Smith and Xie (2013) persist in this view on Tibetan hamster classification. Argyropulo (1933) sorted the Tibetan hamster specimens and divided them into three species (C. kamensis, C. kozlovi, and C. lama), while Lebedev and Potapova (2008) pointed out that C. kozlovi is a subspecies of C. migratorius based on the pseudosciuromorphous structure of the zygomatic plate. Corbet (1978) admitted that both C. kamensis and C. alticola are valid species and other two species (C. kozlovi and C. lama) are synonymous of C. kamensis. Subsequent taxonomic publications on Tibetan hamster adopt the classification conclusion of two separate species (Don & DeeAnn, 2005;Lebedev et al., 2018). However, several studies on Tibetan hamster classification suggested three species (C. kozlovi, C. lama and C. alticola) K E Y W O R D S climatic change, Cricetulus kamensis, demographic history, phylogeography, Quaternary glaciations, species distribution model are subspecies of C. kamensis rather than that of a separate species (Feng et al., 1986;Luo et al., 2000;Wang & Cheng, 1973;Zheng, 1979). In summary, the classification of Tibetan hamster has not been determined.
Prior phylogenetic analysis suggests that C. kamensis is the first species to differentiate among Cricetulus species with a long evolutionary history (Ding, Li, & Liao, 2016). This plateau-dwelling rodent species inhabited high-altitude environments for a long time, which is related to the geological history and climatic change of the QTP. Hence, this species can provide an ideal model organism for the study of the historical biogeography of this region. Herein, we hypothesized that the current distribution of C. kamensis was influenced by physical events (such as the tectonic configuration) and climatic change (such as cold, dry, and glacial periods) of the QTP. In the present study, one nuclear and four mitochondrial genes were used to (a) detect the relationship between the phylogeographical pattern of C. kamensis and the uplift of the QTP, (b) examine both the population genetic structure and demographic history of C. kamensis, (c) assess the likely cause for the decrease of habitats suitable for the survival of C. kamensis during the Late Pleistocene glacial period, and (d) discuss the taxonomical issue of C. kamensis by using phylogenetic analysis and genetic distance.

| Sample collection and ethics statement
A total of 53 C. kamensis were trapped in the QTP area from 13 geographic locations (Table 1; Figure 1). The specimens were identified in the field based on external morphological characteristics (Feng et al., 1986;Luo et al., 2000). Muscle tissue of C. kamensis was collected from these individuals; then, both tissue and voucher specimens were preserved in absolute ethanol and formalin solutions, respectively. Neither sample collections nor the study conflicted with the ethical guidelines, religious beliefs, or legal requirements of China.
All animal experiments conformed to the guidelines for the care and use of laboratory animals and were approved by the Committee of Laboratory Animal Experimentation at Lanzhou University (Lanzhou, China).

| Laboratory protocols
The total genomic DNA was extracted from the muscle samples of C. kamensis using the TIANamp Genomic DNA Kit (Tiangen).
Polymerase chain reaction (PCR) was used to amplify both mitochondrial DNA (mtDNA) and nuclear DNA (nuDNA) fragments of C. kamensis, in which four primer pairs were used to amplify four mtDNA markers (cytochrome c oxidase subunit 1 (Cox1), cytochrome c oxidase subunit 3 (Cox3), cytochrome b (CYTB), and control region (D-loop), Table 2). Two primer pairs were used to amplify partial sequences of von Willebrand Factor exon 28 (vWF) using both standard PCR and nested PCR (Table 2) (Porter, Goodman, & Stanhope, 1996). Nested PCR was used to amplify vWF fragments that required two pairs of primers (outside and inside primer pairs), the outside primer of which was used to enrich the target region from complex genome and the inside primer was used to amplify the target region from the first round PCR products (Shen, Liang, Feng, Chen, & Zhang, 2013). The inside primer pair of vWF was specifically designed for C. kamensis based on other conservative hamster homologous sequences of vWF (Table 2). PCR was employed to amplify DNA fragments using a 25 μl reaction system that contained 2.5 μl 10 × PCR buffer, 1.5 μl 25 mM/L MgCl 2 , 0.5 μl 2.5 mM/L dNTP, 1 μl of each primer pairs (10 μM each), 0.5 μl Taq

| Phylogenetic analysis
A total of 38 haplotypes (53 individuals) were used to build the phylogenetic relationships between 13 geographic populations.
Each concatenated sequence had a length of 5,078 bp, including outgroup taxa. A concatenated sequence consisted of four mtDNA genes (Cox1 + Cox3 + CYTB + D-loop) and one nuclear gene (vWF) fragment, and each marker was separately aligned and these sequences were then manually concatenated using the program Sequence Matrix v1.7.8 (Vaidya, Lohman, & Meier, 2011 (Porter et al., 1996) Out-vWF-L TCGGGGGAGCGTCTCAAAGTCCTGGATGA In this study In-vWF-L GTGACCATGTAGACCAGATTAGG and P. sungorus (Cricetinae, Phodopus) were used as outgroups (Lebedev et al., 2018). The model with the best fit of nucleotide substitution and the best partition schemes for the DNA dataset was identified under the Corrected Akaike Information Criterion (AICc) (Akaike, 1974), which was implemented in the MrModeltest 3.7 program (Posada & Crandall, 1998). Model test result showed that the best-fitted substitution model for the concatenated data was GTR + I + G. 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 1,000 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 or above 0.95 was considered to indicate strong relationships (Leaché & Reeder, 2002). In addition, to represent relationships between haplotypes, the program Network 4.6.0.0 (Bandelt, Forster, & Rohl, 1999) was run to construct clustered relationships for the concatenated mtDNA haplotypes of C. kamensis with a median-joining network approach.

| Nucleotide mutation rate and estimation of divergence time
To Divergence times among the different populations of C. kamensis were estimated using BEAST version 1.7.4 (Drummond & Rambaut, 2007). There were 5,078 bp in the concatenated genes (mtDNA + nuDNA) within the C. kamensis complex. Phylogenetic relationships were reconstructed with the Yule speciation process (Steel & McKenzie, 2001), with enforced monophyly of the ingroup.
MCMC analysis was run twice, chain lengths were 50 million, sampling every 1,000 generations, and the first 10% were discarded as burn-in. One calibration point was used as prior for the time for most recent common ancestor of P. roborovskii and P. sungorus. The age of the P. roborovskii/P. sungorus divergence (5.69 Mya) was used as calibration interval (Lebedev et al., 2018

| Population genetic structure and divergence
Population genetic parameters were calculated for combined mtDNA genes. Haplotypes were assigned to groups based on the phylogenetic relationships between geographic populations. The general parameters of diversity in the combined dataset were calculated for each clade using the program DnaSP v5, including numbers of haplotypes (H) and unique haplotypes (UH), haplotype diversity (Hd), nucleotide diversity (π), and the average number of pairwise nucleotide differences (K). The percentages of unique haplotypes per population were calculated via unique haplotypes divided by the number of haplotypes in a given population. Here, three sample sites were removed because only one sample was available or the nucleotide diversity was zero ( Table 1). Analysis of molecular variance (AMOVA) was conducted using the program Arlequin version 3.5 (Excoffier, Laval, & Schneider, 2007) to identify variation among individuals, populations, and clades, and the significance was tested by 10,000 permutations. The derived populations were grouped based on the phylogeographic structure obtained in the phylogenetic analyses. The population differentiation indexes (F ST ) (Hudson, Boos, & Kaplan, 1992b) between populations were calculated in Arlequin and gene flow (Nm) between populations was calculated using DNASP (Hudson, Slatkin, & Maddison, 1992a;Librado & Rozas, 2009). The genetic distance between observed clades was calculated using MEGA 6 (Tamura et al., 2013) with the Kimura-2-parameter (K2P) model of nucleotide substitution (Kimura, 1980).

| Demographic history
Three methods were used in demographic analyses of C. kamensis.
Potential signals of population growth were detected using Fu's FS (Fu, 1997), Tajima's D (Tajima, 1989), and R2 (RamosOnsins & Rozas, 2002) under a model of sudden expansion, and these parameters were calculated in the program Arlequin and DnaSP. Furthermore, the pairwise mismatch distribution of the combined sequences (mtDNA) was computed in DnaSP to identify the potential signal of demographic expansion of the major clades. A goodness-of-fit test was used to investigate whether the observed data fitted a model of recent expansion (Slatkin & Hudson, 1991). In addition, the raggedness indices (r) were calculated in DnaSP to determine whether the concatenated dataset deviated significantly from the population expansion model (Harpending, 1994). A Bayesian skyline plot (BSP) was implemented in BEAST to estimate the effective population size over time (Heled & Drummond, 2008 (Muscarella et al., 2014;R Core Team, 2017). Moreover, the regularization multiplier value was set up with a range of 2-4, which significantly reduced the difference between calibration AUC (area under curve, AUC), evaluation AUC, and omission rate without affecting the evaluation AUC (Radosavljevic, Anderson, & Araújo, 2014). SDM analysis was run in Maxent, and the data were randomly split into two subsets in which the 80% occurrence data was used as training data and the remaining 20% were used as testing data.
The model used the default convergence threshold, 5,000 maximum iterations, and 10 replications. A receiver operating characteristic curve close to 1 indicated better model performance (Phillips et al., 2006).

| Phylogeographic analysis and population divergence time
The phylogenetic tree showed two major evolutionary clades (A  Table 1 and Figure 1). In addition, median-joining network analysis also indicated the divergence of two major clades, which was consistent with the Bayesian phylogenetic tree (Figure 2b). In conclusion, phylogenetic topologies indicated strong geographic variations between these two clades.
Divergence time analysis also resulted in two major clades that represented two derived populations in the species tree with the same topological structure, corresponding to Bayesian tree ( Figure 2c). As such, the species tree supported these two major

| Population genetic structure
A total of 53 individuals generated 37 haplotypes, which considered alignment gaps and invariable sites based on concatenated mtDNA sequences. Most C. kamensis populations had a higher UH ratio among populations except for the LS population (Table 1). Twenty-six unique haplotypes were found in 12 populations, indicating high levels of genetic variation among C. kamensis populations. The average number of nucleotide differences (K) between geographic populations indicated that population DR had a higher K value (K = 68.571) and showed fine DNA heterogeneity (Table 1). The nucleotide diversities (π) of all 13 populations were quite low (π = 0.001-0.016), and the DR population was the highest (π = 0.016) compared to other populations (Table 1). To reduce the small sample bias, adjacent population were merged with few samples into a bigger number population depending on phylogeographic structure, including two major lineages A and B (Table 3 and Figure 2

| Demographic history
Neutral tests indicated that Tajima's D values were negative and no significant difference (p > 0.05) in lineages A and B were found. They exhibited an observed multimodal mismatch frequency distribution (Figure 3a), suggesting past population did not go through sudden expansion, but instead accepted the constant size model. Fu's Fs were detected in both A and B lineages, and were positive, with neither significant nor unimodal distribution. Those suggest that the past population showed no explicit signals of population expansion or equilibrium (Figure 3a). Moreover, the R2 value was not significant for two derived populations, which also supported the constant expansion model (Table 3)

| Species distribution modeling
The SDM analysis showed an optimal prediction of the C.   (Figure 4b). After the LGM, the suitable habitat was restored to a certain extent based on the potential distribution of MH for C. kamensis. Compared to the MH, the current distribution of C. kamensis was increased, especially in the east and north of the QTP increased (Figure 4c and d).

| Phylogeographic structure: QTP uplift and isolation
The phylogenetic results showed that the haplotypes of C. kamensis from different geographic populations split into two distinct clades (Clades A and B, PP = 1) inferred from the BI and BEAST.
Moreover, the Bayesian tree and median-joining network resulted in the same topological structure (Figure 2a and b) (Wang et al., 2008b), which also blocked the gene exchange between C. kamensis populations.
Furthermore, the rapid uplift of the Nyenchentanglha Mountains also occurred during the Middle Pliocene (about 3.7 Mya) (Wu, Liu, Hu, & Ye, 2003) and also formed a major geographic barrier F I G U R E 4 Species distribution model of suitable habitats for C. kamensis estimated via Maxent. The estimates are inferred from the Last Interglacial (LIG), the Last Glacial Maximum (LGM), the Mid-Holocene (MH), and the current climate (Current). The green area represents the suitable habitat for C. kamensis that impeded the genetic exchange between Clades A and B (Zhu, Zhenhan, Zhao, Jianping, & Wang, 2012).
Prior to the Qing-Zang tectonic movement (3.4 Mya), this region had a typical subtropical climate and a large number of mammalian fossils can be found in the QTP, for example, of Hipparion chilongensis, Coelodonta thibetana, and Panthera blytheae, which confirm climatic assumptions (Deng & Ding, 2015;Li et al., 1979). At relatively low altitude, animals can freely roam the plateau and lowland, and at that time, numerous animals obligately inhabited a warm and humid environment (Deng, 2016). Theoretically, there are wide areas on the plateau that was suitable for the survival of C. kamensis at that time, rather than around the QTP, like the current distribution (Feng et al., 1986;Luo et al., 2000). Three viewpoints can explain the current phylogeographical structure of C. kamensis. First, mountains and rivers were formed by the Qing-Zang tectonic movements, which created geographic barriers that reduced the gene flow between isolated populations and promoted allopatric divergence, which often restricted animal dispersal, particularly that of small mammals with weak diffusion capacity such as rodents . For instance, Fan et al. (2012) suggested that the major clades for the differentiation of Apodemus draco that occurred during the Late Pliocene to Middle Pleistocene correspond to the Qing-Zang tectonic movements and subsequent geological events. Second, the Kunlun-Huanghe tectonic movement in the Middle Pleistocene (between 1.1 and 0.6 Mya) caused the average altitude of the QTP to exceed 3,000 m (a.s.l.) after the Qing-Zang tectonic movement, and a large area of the plateau entered the frozen circle Shi, Zheng, Li, & Ye, 1995).
The subsequent Gonghe tectonic movement raised the QTP as a whole to the current altitude during the Late Pleistocene (about 0.15 Mya) . As a result of the blocking of the northward flow of the Indian Ocean warm current and the strengthening of the winter monsoon, caused drier and colder climate, especially for the northern plateau Lei et al., 2014;Wu et al., 2001). This change of climate caused a dramatic ecological shift in the QTP, where forests were replaced by grasslands, glaciers started to develop, and both deserts and permafrost areas formed (Wu et al., 2001). In addition, the uplift of the QTP resulted in shaping the paleolake in the Qaidam Basin during the Late Pliocene (about 2.5 Mya), and a growing winter monsoon caused the shrinkage and disappearance of ancient lakes and an evolution of desertification (Wang, Huang, & Liu, 2002;Wang et al., 2009;Zeng, Feng, & Cao, 2003). Moreover, the uplift of the QTP also promoted aridification and the development of vast deserts in the Gonghe basin (Liu, Wang, Wang, Hideaki, & Abbott, 2012). Consequently, mountains uplift and aridification promoted the divergence between the Qilian Mountains population and the Bayan har Mountains population.
Third, the uplift of the QTP induced changes of plant communities, which can affect the distribution of C. kamensis (Jia et al., 2012;Liu, Wang, Wang, Hideaki, & Abbott, 2006). These environmental factors caused the favorable habitat of C. kamensis to shrink to margin regions or intermountain basins with relatively low altitude in the QTP. C. kamensis show a strong habitat preference, which makes them highly sensitive to environmental change. Consequently, it was easy to understand that the uplift of the QTP caused isolation among populations. A lack of gene exchange between populations let each population evolve in its own direction. In summary, all of these rational factors can lead to the current phylogeographical pattern of C. kamensis.

| Genetic structure and divergence of the C. kamensis population
All total of 53 C. kamensis specimens were collected from 13 geographic populations in the QTP, all of which were part of two recovered populations and were subjected to population genetic analysis.
It is worth noting that there was no shared haplotype among populations, and each population had higher UH ratio, and individuals in the population had higher nucleotide differences (K value, see Table 1).
Moreover, the recovered populations showed similar results based on the phylogeographical patterns of C. kamensis (Table 3). Here, 13 geographic populations and two derived populations showed low genetic diversity (Tables 1 and 3). AMOVA indicated highly restricted gene flow and fragmented divergence throughout the distributional range of the species ( C. kamensis preferred to live in valleys or riversides with leguminous or chenopodiaceous plants (Luo et al., 2000), and SDM analysis further validated its habitat preferences. Thus, geographic barriers restrict the gene exchange between populations, such as Nanorana pleskei (Dicroglossidae) (Zhou et al., 2017). The reason for the high UH ratio and nucleotide differences in C. kamensis populations was that its population density is comparatively small and it has been listed as near threatened in a recent survey (Jiang et al., 2016). The C. kamensis populations showed low genetic diversity (Tables 1 and   3), which was related to the harsh environment on the QTP.  suggested that organisms living at high-altitude areas with high selective pressure develop a specific genotype for survival. In contrast, suitable environments with relaxed selective pressure can accommodate more genotypes and the population possesses higher genetic diversity , such as Phrynocephalus vlangalii.
Thus, the genetic structure of the C. kamensis population is strongly affected by the harsh environment of the QTP.

| Demographic history: Glacial oscillations and population colonization
Neutrality tests for the derived populations indicated that Fu's FS, Tajima's D, and R2 were not significantly negative, and mismatch distributions were multimodal rather than following a Poisson distribution, suggesting that two recovered populations did not go through a sudden expansion in history (Table 3; Figure 3). BSP analysis indicated that the population size fluctuating slightly at the Late Pleistocene (about 0.10 Mya) and showed an increase in population since then (Figure 3). The reason why the C. kamensis population had not experienced rapid expansion lies in the complicated landform configuration and the harsh environment after the uplift of the QTP.
Both extrinsic and intrinsic factors interacted to generate those results. Intrinsic factors include that the population density of C. kamensis is low and that they are not the dominant species on the QTP (Feng et al., 1986;Jiang et al., 2016). This hamster prefers to live alone and has a lower reproductive capacity in comparison with gregarious rodents such as Merions meridianus (Rodentia, Gerbillinae) (Luo et al., 2000). Furthermore, its habitat preference and weaker locomotion were also potential intrinsic factors. Extrinsic factors were that the uplift of the QTP resulted in a decrease in suitable habitat of C. kamensis. The species was restricted to microhabitats and dispersed in valleys and riversides of the QTP (Figure 4). At this point, the Quaternary glacial period exerted little effect on the population size of C. kamensis ( Figure 3) because valleys or riversides are generally at lower altitudes and can temporarily help animals through the harsh environment (Fan et al., 2012;Qu et al., 2010).
The last glacial period only caused a fluctuation in the C. kamensis population size during the Late Pleistocene (Zhu et al., 2012), which was confirmed by BSP analyses (Figure 3). Areas of low altitude (e.g., valleys, riversides, and intermontane basins) are not covered by ice sheets and can thus provide significant suitable microhabitats for organisms during the glacial period (Zheng & Rutter, 1998). A variety of vertebrates in the QTP benefit from glacial refugia, and the populations began to expand after the end of the glacial period, such as Montifringilla adamsi and Pyrgilauda blanfordi (Qu et al., 2010), A. draco (Fan et al., 2012), and N. parkeri (Liu et al., 2015). Therefore, C. kamensis survived the Late Pleistocene glacial period.
The results of SDM analyses indicate that the LGM induced a significant decrease in the C. kamensis suitable habitat in comparison with the LIG (Figure 4). The LGM was larger than the previous glacial period so that suitable habitats of C. kamensis may have been insufficient for their survival (Shi et al., 1998). Gupta, Sharma, & Shah (1992) and Zheng and Rutter (1998) suggested that the whole plateau is covered by a huge ice sheet during the glacial period, and force most species to retreat to favorable habitat at lower elevations on the edge of the plateau during glacial maxima, from which they recolonized the interior during interglacials (Fan et al., 2012;Lei et al., 2014;Tang et al., 2010). Several recent phylogeographic studies support this hypothesis, SDM analysis showed that suitable habitats of N. pleskei were also significantly reduced during the LGM (Liu et al., 2015;Wang et al., 2017). Later, suitable habitats for C. kamensis survival were restored since the Mid-Holocene, and the suitable habitat was able to expand again (Figure 4). The reason was that the warm and humid climate on the QTP was suitable for animals and plants to thrive during the Mid-Holocene, suitable habitats may have increased since the LGM, and consequently, organisms living at the QTP could expand (Shi et al., 1998). In addition, the precipitation in the east of the plateau has decreased from the southeast to the northwest, since a region with high precipitation is located in eastern Tibet and the western Sichuan plateau, while a region with low precipitation is located in the Qaidam Basin (Hu & Liang, 2013).
Theoretically, within a certain range of precipitation, the diversity of rodents correlates positively with precipitation (Abramsky & Rosenzweig, 1984;Ye, Ma, & Feng, 1998). A case study suggested that small mammals from the semiarid environment are more likely to choose a habitat with higher precipitation (Milstead et al., 2007).
This explained why the suitable habitat of C. kamensis increased in the east of the QTP under the current climatic conditions.

| Review on the C. kamensis classification
The classification of Tibetan hamster has been controversial over the past century. Historically, Satunin named the female hamster from the southeastern part of the Q inghai-Tibet Plateau for the first time as an Urocricetus kamensis (Type locality: the Moktschjun River in the upper reaches of the Lantsang River) in 1902 and named the specimen with a pale pelage from Qilian Mountain as the C. kozlovi (Satunin, 1902). Then, Bonhote named a specimen from Lhasa (Tibet autonomous region, China) as a C. lala in 1905 (Bonhote, 1905), and Thomas named the specimen with shorter tail and pale pelage from the Ladakh region adjacent to the northern Qinghai-Tibet Plateau as the C. alticola (Thomas, 1917 Corbet (1978) admitted that C. kamensis and C. alticola are valid species, and the other three species (C. kozlovi, C. lama and C. tibetanus) are synonyms of C. kamensis (Corbet, 1978 tunately, a small sample size can cause errors in its classification (Mayr, Linsley, & Usinger, 1953), such as age variation, seasonal variation, and ecological variation. Wang et al. (1973) pointed out C. alticola is likely to a subspecies of C. kamensis, and then Zheng (1979) confirmed that this conclusion is correct in the scientific investigation of animals and plants in the Ali area of Tibet (China). Feng et al. (1986) and Luo et al. (2000) also suggested that C. alticola is a subspecies of C. kamensis rather than a valid species. In conclusion, whereas values between 2% and 11% had a high probability of being indicative of conspecific populations or valid species. The average K2P distance of sister species was 9.55% (Bradley & Baker, 2001), suggesting that two major clades (A and B) are considered arbitrary as a species. It was worth noting that the K2P distance between hamsters was in the range of 16.7%-26.6% that was much higher than two clades (6.7%). The K2P distance between two major clades of the M. meridianus complex is 8.6% and does not meet interspecific divergence level, for instance, suggesting that smaller degrees of divergence tended to characterize major within-species phylogroups . Our results indicated that the K2P distance between C. kozlovi and C. migratorius was 0.219 and the K2P distance from C. kamensis was 0.005, in addition, suggesting that C. kozlovi is not a subspecies of C. migratorius but a C. kamensis, which is not consistent with Lebedev and Potapova (2008)'s inferences. In summary, the K2P values did not support two clades to be regarded as a separate species and suggested the intraspecific value.
Theoretically, it appears that phenotypical changes are regarded as responding to environmental and climatic changes (Mayr et al., 1953), such as the posterior part of the skull, coat color, and teeth morphology, whereas these have long constituted the basis of rodent classification and should be validated by molecular systematic method (Chevret & Dobigny, 2005;Ge et al., 2012).
Ecological specialization triggers an uncoupling of molecular and phenotypic evolution, and the departure from a phenotypic drift pattern (Qu et al., 2013;Renaud, Chevret, & Michaux, 2007), such as climate, vegetation, and diet. Thus phenotypic plasticity actually represents a fundamental component of evolutionary change, and adaptive plasticity may involve both physiological homeostasis and morphological response (Thompson, 1991). These were closely related to environmental variation. Therefore, previous taxonomic conclusions only with morphological methods did not necessarily reflect authentic phylogenetic relationships. Currently, the taxonomic approaches intend to use DNA which can be as an auxiliary criterion for identifying a species or taxon (Arctander, Arctander, Minelli, Thomas, & Vogler, 2003). Morphologically, the differences between the four Tibetan hamsters are mainly concentrated on tail length and pelage variation (Feng et al., 1986;Luo et al., 2000). These morphological characteristics may be related to habitat surrounding, such as Rattus species (Yom-Tov, Yom-Tov, & Moller, 1999). Here, the divergence time between clades A and B did not support the evolution of them into a new species (Figure 2c) because for a vertebrate to evolve into a new species to it takes at least 2 million years on average and this timeframe is too short to evolve a new species (Avise, Walker, & Johns, 1998).
In addition, the genetic distance failed to the interspecific level based on the complete CYTB data. Based on our results, taken together, four hamsters were synonyms of C. kamensis. Further work is needed to assess morphological variation and subspecific differentiation on C. kamensis.

ACK N OWLED G M ENTS
The work was supported by the National Key Research and Development Program of China (2017YFC0504802) and the National Natural Science Foundation of China (Nos. 30870294, 31372179). We were grateful to Guangjie Luo, Dongdong Jing and Quan Zhou of School of Life Sciences, Lanzhou University, who provided helps in the process of sampling. In addition, we thank to Professor Shaoying Liu (Sichuan Academy of Forestry, Chengdu) for his contribution to specimen in this study.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.

AUTH O R CO NTR I B UTI O N
JC L conceived the study; L D was responsible for obtaining and analyzing data, and wrote the manuscript with the help of JC L.