Phylogeography of the South China Field Mouse (Apodemus draco) on the Southeastern Tibetan Plateau Reveals High Genetic Diversity and Glacial Refugia

The southeastern margin of the Tibetan Plateau (SEMTP) is a particularly interesting region due to its topographic complexity and unique geologic history, but phylogeographic studies that focus on this region are rare. In this study, we investigated the phylogeography of the South China field mouse, Apodemus draco, in order to assess the role of geologic and climatic events on the Tibetan Plateau in shaping its genetic structure. We sequenced mitochondrial cytochrome b (cyt b) sequences in 103 individuals from 47 sampling sites. In addition, 23 cyt b sequences were collected from GenBank for analyses. Phylogenetic, demographic and landscape genetic methods were conducted. Seventy-six cyt b haplotypes were found and the genetic diversity was extremely high (π = 0.0368; h = 0.989). Five major evolutionary clades, based on geographic locations, were identified. Demographic analyses implied subclade 1A and subclade 1B experienced population expansions at about 0.052-0.013 Mya and 0.014-0.004 Mya, respectively. The divergence time analysis showed that the split between clade 1 and clade 2 occurred 0.26 Mya, which fell into the extensive glacial period (EGP, 0.5-0.17 Mya). The divergence times of other main clades (2.20-0.55 Mya) were congruent with the periods of the Qingzang Movement (3.6-1.7 Mya) and the Kun-Huang Movement (1.2-0.6 Mya), which were known as the most intense uplift events in the Tibetan Plateau. Our study supported the hypothesis that the SEMTP was a large late Pleistocene refugium, and further inferred that the Gongga Mountain Region and Hongya County were glacial refugia for A. draco in clade 1. We hypothesize that the evolutionary history of A. draco in the SEMTP primarily occurred in two stages. First, an initial divergence would have been shaped by uplift events of the Tibetan Plateau. Then, major glaciations in the Pleistocene added complexity to its demographic history and genetic structure.


Introduction
The Tibetan Plateau is a recognized biodiversity hotspot [1]. Recently, a number of phylogeographic or phylogenetic studies have focused on the Tibetan Plateau due to its complex geologic history [2][3][4][5][6][7][8]. Historically, it uplifted several times, and the most intense uplift events were known as the Qingzang Movement (3.6-1.7 million years ago, Mya) and the Kun-Huang Movement (1.2-0.6 Mya) [9][10][11]. These events caused dramatic environmental changes: grasslands replaced forests and glaciers and deserts developed as the climate gradually became drier, colder and windier [12]. Thus, as a result of the uplift events, the plateau began to diversify topographically by the continual developed mountains and watercourse [13]. Moreover, it also caused its elevation to increase by up to 3000 m, which consequently increased the exposure of the endemic fauna to glacial events [5,9]. So, Quaternary glacial cycles have been regarded as an important factor in promoting the diversification of species in this region [2,14]. Although there are different opinions about the range and extent of glaciers on the plateau [15][16][17], ice-sheets during the most extensive glaciation (0.5-0.17 Mya) were undoubtedly several times larger than those found today [9]. Taken together, the current distribution of endemic organisms and their patterns of genetic differentiation should be greatly influenced by these geologic and climatic events.
Several studies have demonstrated that species distributed on the Tibetan Plateau retreated to the east side of the plateau during cold periods [2,7,18,19]. Qu et al. [19] reported that both the twite (Carduelis flavirostris) and the black redstart (Phoenicurus ochruros) inhabited the northeastern margin of the Tibetan Plateau and maintained stable population sizes through glacial cycles. On the contrary, three other avian species (Montifringilla adamsi, Pyrgilauda blanfordi and Eremophila alpestris), mostly distributed on the inner region of the plateau, experienced population expansions after glacial retreat [19]. These results indicate that the east side of the plateau probably constituted a large refugium during the major glaciations in the Pleistocene.
In this study, we focused on the southeastern margin of the Tibetan Plateau (SEMTP), including western Sichuan Province, northwestern Yunnan Province and eastern Xizang Autonomous Region (Tibet). The large rivers and mountains (e.g. Daxueshan Mountains, Shaluli Mountains, Gongga Mountains, Dadu River, Yalong River, and Jinsha River) in this region could be natural barriers for the postglacial colonization of small mammals (Fig. 1). The magnificent Gongga Mountain, with the highest peak of 7,556 m above the sea level, was an important glaciation centre and contemporary glaciers are still well developed [20]. Conversely, some previous studies indicated that Gonga Mountain region (GMR) was a suitable refuge during glacial cycles [21,22]. Many ancient and endemic species inhabited in this region [21,23]. But so far, none of the studies could provide detailed information.
In order to further assess the role of geologic and climatic events on the Tibetan Plateau in shaping the genetic structures of local fauna, we present a phylogeographic investigation of the South China field mouse (Apodemus draco). Small rodents are good candidates for inferring biotic history from contemporary patterns of genetic variation because their short generation time and relatively limited mobility [24]. This mouse is an endemic species to China, and is widely distributed from west to east, including in the Sichuan, Yunnan, Gansu, Qinghai, Shannxi, Hubei, Hunan, Fujian and Taiwan Provinces [23,25]. However, the eastern part of the Tibetan Plateau was the centre of radiation of the genus Apodemus [18,26]. One recent work sequenced the mitochondrial cytochrome b (cyt b) gene in six A. draco and 197 A. ilex with samples from Yunnan-Guizhou Plateau in China, they reported that A. ilex should be a valid species rather than synonym of A. draco [27]. They also indicated that the phylogenetic pattern of A. ilex was closely related to the complex geographical structures in southwest China [27].
The cyt b gene sequence has been extensively and successfully used in phylogeographic studies for rodents [27][28][29][30][31][32][33]. Therefore, the nucleotide sequences of the cyt b gene were applied to investigate the genetic variation and evolutionary history of A. draco. We examined: (1) the relationship between the genetic structure of A. draco and geologic or climatic events on the SEMTP, and (2) whether the SEMTP, especially the GMR, potentially functioned as a refuge during the Quaternary period. Details of the historical demography of A. draco revealed in this study could provide new insights into how this rodent and other small mammals in this region responded to the geologic or climatic events on the SEMTP.  Table 1. Major rivers and mountain are also shown. The distribution of observed major evolutionary clades in phylogenetic analyses is shown in different colors. The two glacial refugia (Gongga Mountain and Hongya County) were labeled in red color. doi:10.1371/journal.pone.0038184.g001

Ethics Statement
We only used wild rodents in this study. All samples were obtained following the regulations for the implementation of China on the protection of terrestrial wild animals (State Council Decree [1992] No. 13). All our observational and field studies and lab work were approved by Wildlife Protection Office, Sichuan Provincial Forestry Departments (China) and by the Ethics Committee of Sichuan University, China.

Samples and DNA Extraction
One Hundred and three South China field mice were collected from 47 sampling sites in Sichuan Province (Table 1). In addition, 23 cyt b sequences from GenBank were retrieved (GenBank Accession Nos. AB096825, AB109397, AM945800-AM945819, AY389007), which were collected from an additional six sampling sites in the Sichuan and Yunnan Provinces (Table 1) [18,34,35]. We combined the very near sampling sites (which located at the same mountain or village) into one population. In total, we used 126 individuals from 20 populations (Table 1). Cyt b sequences from more individuals were in GenBank, but the length of the sequences was significantly shorter (319 bp or 402 bp) than our target sequences (1050 bp), thus we did not include them here.
Detailed information on sampling localities is shown in Fig. 1 and Table 1. Voucher specimens were deposited in the Museum of the Sichuan Academy of Forestry. Total genomic DNA was extracted from muscle or liver tissues preserved in 95% ethanol using standard procedures [36].

Amplification and Sequencing
Polymerase chain reaction (PCR) was used to amplify the sequences of cyt b with the universal primers L14724/H15915 [37]. PCR conditions included an initial denaturation for 5 min at 95uC, 35 cycles of denaturation for 40 s at 94uC, 45 s annealing at 50-60uC, and a 1 min extension at 72uC followed by a final extension at 72uC for 10 min. The amplification was performed in 25 ml reaction volumes with 2.5 ml of 106 EX Taq buffer (Mg 2+ Free; TaKaRa Biotech), 1.5 ml dNTP (2.5 mM each), 2.0 ml MgCl 2 (25 mM), 1 ml of each primer (10 uM), 0.2 ml EX Taq polymerase (5 U/ml, TaKaRa Biotech), and approximately 200 ng total genomic DNA as template. PCR products were sequenced from both directions with the same PCR primers in an ABI PRISM 3730 DNA sequencer (Applied Biosystems, Inc.).

Genetic Diversity and Protein Sequences
We aligned the sequences in software MEGA 4.0 [38]. Nucleotide saturation was tested by plotting transitions and transversions against the JC69 distance using DAMBE [39].
Haplotype diversity (h) and nucleotide diversity (p) were estimated using DnaSP Ver. 5 [40]. Genetic differentiation between populations was evaluated by estimating pairwise values of F ST with the software Arlequin 3.1 [41]. The significance was assessed by 1000 permutations. Populations contain only one sample were excluded in this analysis. Pairwise comparisons were used as genetic distances in the following analyses, such as isolation by distance. McDonald-Kreitman test [42] was used to test the overall selection on cyt b on the website (http://mkt.uab.es) [43]. In this analysis, the numbers of synonymous and nonsynonymous substitutions in cyt b of A. draco and A. latronum (Nos. GU908339-GU908394) were estimated and their ratio was compared for each species with the same ratio for fixed mutations between the species. The cyt b sequences of A. latronum were used here due to A. latronum was the sister taxon to A. draco.

Phylogenetic Analysis
Two different methods of phylogenetic analysis were used to reconstruct phylogenetic relationships among cyt b haplotypes: Bayesian inference (BI) and maximum parsimony (MP). Apodemus latronum (GU908426-GU908428, GU908433, GU908436) was used as outgroup due to it was the sister taxon to A. draco. In order to avoid the possible bias, A. flavicollis (AB032853) and A. sylvaticus (AB033695) were also chosen as outgroups. The best fit model of DNA substitution was obtained using MrModeltest2.2 [44] under the Akaike Information Criterion (AIC). BI was performed with MrBayes 3.1.2 [45]. For BI, we used two different datasets, partitioned and non-partitioned. For partitioned dataset, the nucleotides sequences were partitioned by codon position. In total, three partitions were given; one each for codon positions 1, 2, and 3 cyt b. Posterior distributions were obtained by the Markov Chain Monte Carlo (MCMC) method with one cold chain and three heated chains for 10,000,000 generations and sampled every 1000 generations. The first 25% were discarded as a conservative burn-in and the remaining samples were used to generate a 50% majority rule consensus tree. PAUP* 4.0 b 10 [46] was used to perform MP analysis. A heuristic search strategy was employed with the tree bisection and reconnection (TBR) branch swapping algorithm, random addition of taxa and 1000 replicates per search. Supports for branches in the MP trees were tested by bootstrap analysis with 1000 replicates.
Since the haplotype network could better detect the relationship among haplotypes, we also constructed a median-joining network (MJN) approach to depict relationships among the cyt b haplotypes of A. draco using Network 4.6.0.0 [47].

Landscape Analyses
We used analysis of molecular variance (AMOVA) to estimate the geographical pattern of population subdivision using different grouping options [48]. To test for the correlation between genetic differentiation and geographic distance, we performed a Mantel test [49] and a spatial autocorrelation analysis. In Mantel test, we used the Euclidean distance between sampling localities, based on GPS coordinates, as the geographic distance. Because we did not have the coordinates of the samples from GenBank, we approximated distance using mountains or county locality information listed in the corresponding papers [18,34,35]. We used three different datasets to perform the Mantel test. First one included 13 populations (exclude populations that only contain one sample); second one further excluded all the samples from other studies since we did not have their exact coordinates. Third one only contained the populations in clade 1 (Jiajin, Gongga and Erlang Mountains, and Hongaya county). The significance of the Mantel test was determined by 1,000 permutations using the isolation by distance (IBD) web service v3.16 (http://ibdws.sdsu. edu) [50]. The spatial autocorrelation analysis was performed in Allele in Space (AIS) [51]. The spatial autocorrelation algorithm calculates the statistic Ay, the average genetic distance between individuals in distance class y. Thus, Ay can be interpreted as the  average genetic distance between pairs of individuals that fall within distance classes [51]. In order to ensure that the arbitrary defined distance class size had no effect on result, analyses were performed using 5, 10 and 15 distance classes, respectively. We used a randomisation procedure consisting of 5,000 replicates to identify distance classes where average genetic distances were significantly larger or smaller than random expectations. Because this algorithm uses genetic and spatial data for every individual, we preformed this analysis on different data set, all individuals and individuals without GenBank samples based on the same reason we explained above.

Molecular Dating
Divergence times among the observed mitochondrial clades on phylogenetic tree were estimated by BEAST v1.7.1 [52] using a Bayesian coalescence approach. We also used partitioned (by codon positions) and non-partitioned data. For partitioned data, we had two different partitioned datasets here. The first one had two partitions (pos1+pos2, pos3), and the second one had three (pos1, pos2, pos3). We used an uncorrelated lognormal relaxed clock with the GTR+I+G model of substitution. All analyses were performed using a coalescent-constant size process of diversifica-tion. The prior for mean mutation rate was specified as a lognormal distribution, with a mean of 0.02, log (standard deviation) of 0.56 and offset of 0.0172. This distribution covered the range from 2% to 8% substitutions per site per Myr (million years), with the 95% highest posterior density (HPD) range of 2.4% to 6.01% per Myr. We used this prior because the standard divergence rate for the mammalian cyt b gene is estimated at 2% per Myr [53], and the evolutionary rate of cyt b in genus Apodemus is estimated at 2.4% per Myr [34,54]. However, a number of studies [28,29,55,56] indicate that the cyt b gene of rodents in general evolved several times faster than the standard rate. Moreover, Fan et al. [22] applied a rate of 2.4-8% per Myr for a phylogeographic analysis of another field mouse, A. latronum in their study. Monophyly of clades was not enforced, and runs were initiated on random starting trees. The MCMC chain was run for 50 million generations, sampling every 1,000 generations. The convergence of the chains to stationary was checked using Tracer 1.5 [57] and the first 25% were discarded as burn-in. The final molecular clock tree was produced in TreeAnnotator 1.7.1 [58] using the mean as node heights. Finally, Figtree v1.3.1 [59] was used to visualize the results. Although the divergence times were Figure 2. Fifty percent majority rule consensus tree from Bayesian analysis of Apodemus draco based on all the haplotypes of cyt b data. Numbers represent node supports inferred from Bayesian posterior probability of non-partitioned data and partitioned data, and maximum parsimony bootstrap analyses, respectively. Only values of the main evolutionary clades are shown. The haplotype names can be found in Table 1. doi:10.1371/journal.pone.0038184.g002  Table 1. Evolutionary clades correspond to the five major clades and three subclades in Fig. 2 only rough estimates here, it will provide useful information when we discuss A. draco' s divergence.

Historical Demography
The demographic history was tested by mismatch distributions and neutrality tests based on the cyt b sequences. First, two tests of selective neutrality were calculated for each set (whole sample and each observed phylogenetic clade). Fu's Fs [60] was calculated in Arlequin 3.1 [41], and Ramos-Onsins and Rozas's R 2 statistics [61] was calculated in DnaSP Ver. 5 [40]. Fu's Fs and R 2 were applied due to the fact that they have been suggested as the most powerful tests for detecting sudden population growth or contractions [61,62]. Fu's Fs is more powerful for large population sizes, whereas R 2 works well for small ones [61]. The sum of squared deviations (SSD) and the raggedness index (r) [63] were also estimated in Arlequin 3.1 to determine whether the sequences deviated significantly from a model of population expansion. Significance of all the statistics was calculated using 1,000 coalescent simulations with the software Arlequin 3.1 and DnaSP Ver. 5. Mismatch analyses [64] were carried out using DnaSP Ver. 5. The time of possible population expansions (t) was calculated according to t = 2 ut [65], where t was the model of the mismatch distribution and u was the mutation rate of the DNA sequence. The value u was derived from u = mk, where m was the mutation rate per nucleotide and k was the number of nucleotides. We used an estimated evolutionary rate of 2-8% per Myr for cyt b (rationale is provided in previous section).

Sequence Information
We generated 103 cyt b gene sequences each consisting of 1050 sites (GenBank Accession Nos. HM162766 -HM162833, JQ424902 -JQ424910). With sequences from the additional 23 individuals, we recovered 76 haplotypes from 126 individuals, and the overall h was 0.982, while the p was 0.0381 (Table 2). Saturation analysis did not show any evidence of saturation, so our sequences were used for further analyses. We were confident that the cyt b sequences obtained in our study represented the mitochondrial cytochrome b gene since all cyt b sequences could be translated into amino acid sequences without interruption and closely matched the previously published cyt b sequence for A. draco. The McDonald-Kreitman test did not provide overall evidence for positive selection on the cyt b gene [42,43]. The value of P n /P s (polymorphism) was larger than D n /D s (divergence), and the difference between them was not significant (Table 3).

Phylogeographic Structure
Based on the AIC test, the GTR+I+G model was chosen for the cyt b. For partitioned data, different models of sequence evolution were detected (pos1: SYM+I; pos2: HKY+I; pos3: GTR+I+G). The 50% majority consensus tree based on BI was shown in Fig. 2. MP tree recovered the same tree topology and was not shown. The phylogenetic trees based on partitioned and non-partitioned datasets exhibited the same topology.
The haplotype network also revealed five clades separated by 9 to 73 mutational steps (Fig. 3), which was consistent with the phylogenetic trees. The subclades within clade 1 only separated by four to six mutational steps, which exhibited very close relationships. Moreover, haplotype YX2 was found at samples both from Yuexi and Zhongguo, HY1 was identified in individuals from Hongya and Erlang Mountain, and JM3 was observed in individuals from Jiajin Mountain and Gongga Mountain.

Genetic Variation and Landscape Genetic Analyses
The F ST values for pairwise populations comparisons ranged from 20.026 to 0.985. Of the pairwise values obtained, only 3 of the 78 were nonsignificant ( Table 4). The AMOVA analysis showed a very strong differentiation among different populations or evolutionary clades, and showed that among populations or clades variation accounted for 66% (p,0.001) or 86% (p,0.001) of the overall variation. However, when we focused on the populations in clade 1, the result exhibited that among populations variation only accounted for 26% of the overall variation, but 48% of the overall variation was contributed by variation among samples (Table 5).
Three different datasets were applied to perform the Mantel test. A significant correlation between genetic (F ST ) and geographic distances was detected in first two datasets (13 populations: p = 0.0009, r = 0.4122; 11 populations: p = 0.0006, r = 0.4944), indicating high levels of geographical structure for A.draco. But for the dataset only including 4 populations in clade 1, the result was non-significant (p = 0.3800, r = 0.1831). Spatial autocorrelation analysis (SAA) also found a significant relationship between genetic and geographic distances in different distance classes (5, 10 and 15, only show the results of 10 and 15). In analyses performed using 10 distance classes, both the whole individuals dataset and without GenBank samples dataset were significant smaller than random expectations over the first two shortest distance classes (geographic distance up to ,200 and 150 km, respectively). Then, when the geographic distance increased (start from 300 and 200 km, respectively), both the datasets yielded significant values of Ay that were higher than random expectations. Similar results were obtained when 5 and 15 distance classes were used for analyses (Fig. 4). Different results were observed when we performed the analyses with individuals only in clade 1. First, it yielded significant smaller value than random expectations at the shortest distance class (,20 km). Then, significant higher values than random expectations were observed at around 60-80 km. But this dataset showed significant smaller and higher values again at around 100 km and 180 km, respectively (Fig. 4).

Estimation of Divergence Times
Although we applied three different datasets to estimate the divergence time, we got the very similar results. Therefore, we described and discussed the divergence times based on the three partitions dataset ( Table 6). The divergence time between clade 5 and the remaining clades occurred during the Pliocene/Pleistocene boundary at approximately 2.20 Mya, although the range of the 95% HPD interval spanned the period from the late Pliocene through the middle Pleistocene, between 3.80 and 0.87 Mya ( Table 6). The divergence between clade 4 and its sister clade and clade 3 and its sister clade occurred in the middle Pleistocene at

Historical Demography
We did not include the data sets of whole sample and clade 1 in mismatch distribution because the phylogenetic trees indicate they exhibit a high level of population structure. But the non-significant values of Fu's Fs and Ramos-Onsins and Rozas's R 2 indicated that the whole population maintained a stable population size. The mismatch distribution analysis indicated a ragged and multimodal distribution ( Table 2), suggesting a relatively constant population size for the clades 2, 3, 4 and 5, this result was also supported by non-significant values of Fu's Fs and R 2 tests for clades 2 and 5 ( Table 2). However, Fu's Fs for clade 4, and Fu's Fs and R 2 for clade 3 showed significant values. Thus, considering the multimodal distribution in mismatch distribution of clades 3 and 4, it was hard to conclude whether they experienced a population expansion.
The Subclade 1C showed multimodal distribution and nonsignificant value in R 2 , but a significant Fu's Fs value. This inconsistent pattern maybe caused by single sampling region of this subclade. We did not calculate the possible expansion time for clade1, because historical expansion in subclade 1A and subclade 1B was also supported by unimodal mismatch distributions and significant negative values of all the neutrality tests ( Table 2). The values of t were different in subclades 1A and 1B, indicating that

Discussion
Sakka et al. [35] studied the phylogeography of four Apodemus species (A.draco, A.latronum, A.peninsulae, and A.agrarius) in the Asian Far East. Results from our analyses were different from this study in several ways. First, in their study, the cyt b haplotypes from Yunnan Province did not cluster with any other haplotypes from Sichuan Province, and they formed a single lineage [35]. However, in our study, individuals from Yunnan Province (populations 2 and 3) clustered together with individuals from Daocheng (population 18), Sichuan Province (Fig. 2). Second, haplotypes from Maerkang, Emei Mountain and Baoxing fell into one lineage [35], but phylogenetic tree in our study clearly showed that individuals from Maerkang (population 4) formed clade 2, and samples from Emei Mountian (population 1) and Jiajin Mountain (population 6) fell into subclade 1A and subclade 1C, respectively (Fig. 2). It is possible that their [35] limited numbers of individuals and sampling sites (individuals: 13 in Yunnan and 22 in Sichuan; sampling sites: 4 in Yunnan and 4 in Sichuan) caused these differences. Based on the topographic complexity and unique geologic history of the SEMTP, it is necessary to use more samples in this region to investigate the relationship between geologic events and local species' evolutionary history. Apodemus draco exhibits a high level of population structure in our study. Most individuals from the same or nearby sampling localities clustered together, with the exception of individuals from Gongga Mountain (in both clades 1 and 2) and individuals from Jiajin Mountain (in different subclades of clade 1). The results of multiple landscape analyses further support the geographic structuring of genetic diversity in this species. As shown in Figure 1, some big rivers and mountains would act like barriers, such as Yalong River, Dadu River, and Qionglaishan Mountain and Jiajin Moutain. But other rivers and mountains did not act as barriers, such as Daxiangling Mountain. Hongya County (population 9) and Emei Mountain (population 1) located at right side of Daxiangling Mountain, but their samples were clustered with samples from its left side (populations 6,11,12,15,16). Therefore, the phylogeographical structure of A. draco was very complicated due to the complex geologic history and topography of the SEMTP. However, an unexpected pattern was observed in clade 4. Individuals from the central part of Sichuan Province, Zhongguo (population 5), clustered together with samples from southern Sichuan Province, including Jinyang, Zhaojue, Yuexi, Jiulong and Shimian Counties (populations 7,8,10,13,14) and the Gongga Mountain (population 15), whereas the samples from Hongya County (population 9) and Emei Mountain (population 1), which were more close to the sampling sites in southern Sichuan, did not fall into clade 4 (Fig. 1). Further studies are needed to elucidate the unique pattern found in clade 4.
The estimated divergence times of A. draco's main splits (2.20-0.55 Mya) were congruent with the periods of the Qingzang Movement (3.6-1.7 Mya) and Kun-Huang Movement (1.2-0.6 Mya) [9][10][11]. These two movements were known as the most intense uplift events in the Tibetan Plateau. These events shaped the geologic configuration of this area, making the topography of the plateau very complex [13]. Because small rodents are sensitive to the changing environment, possess limited migration abilities [22,24], therefore, many phylogenetic and phylogeographic studies [3,4,14,18,22] have demonstrated that uplift events and environmental complexity played an important role in the diversification of local species. It is reasonable to suggest that the initial divergence of A. draco would have been shaped by the uplift events of the Tibetan Plateau, and high levels of genetic diversity have beeen maintained as a result of these events. Although the divergence times in our study was rough estimates, this conclusion was supported by Sakka et al. [35]'s study. In their work, the divergence time between the main lineages of A.draco corresponded to periods from 1.5-1.7 Mya to 2.0-2.2 Mya, and they further suggested that the complex land conditions could explain the divergence of A. draco in Sichuan Province.
The Tibetan Plateau underwent four or five glaciations during the Quaternary period [66], and the frequency and range in this plateau was less than that in other regions in Asia [66,67]. According to geologic data, the extensive glacial period (EGP) on the Tibetan Plateau occurred during the middle Pleistocene (about 0.5 Mya), and continued until 0.17 Mya [66,68,69]. The last glacial period (LGP) occurred 0.08-0.01 Mya, in which the last glacial maximum (LGM) was 0.021-0.017 Mya [9,70]. During the EGP, ice coverage could be permanent at high elevations and central regions of the Tibetan Plateau [69,70,71]. Conversely, the frequency of glaciers in the east was less than that in the west, and indeed, no large ice sheets were found in the SEMTP, so many continuously ice-free areas were likely to have existed in the SEMTP because of its lower elevation and the complicated geologic configuration [66,67,70,71]. Consequently, many areas located in the SEMTP could provide suitable refugia, allowing A. draco to maintain stable population sizes.
Indeed, some recent phylogeographic studies have indicated that the large mountain ranges on the east edge of the Tibetan Plateau could provide substantial refugia for the local species during the Quaternary [19,22,72]. However, this did not mean the major glaciations in the Pleistocene did not have any effect on the evolution of A. draco. The estimated divergence time showed that the split between clade 1 and clade 2 occurred at 0.26 Mya, which fell into the EGP (0.5-0.17 Mya). Thus, it was possible that the largest glaciation further isolated these groups, which could promote the differentiation. Similarly, Zhan et al. [72]'s study showed that the major glaciations during the Pleistocene have had a major impact upon the evolution of the blood pheasant (Ithaginis cruentus).
The IBD and SAA results showed that there was no clear relationship between geographic distance and genetic distance in the central region (clade 1), although there were many huge mountains. It seems that there is a certain connection between the pattern in clade 1 and the regional historical demography. Indeed, two of the subclades in clade 1 experienced population expansions.
Subclade 1A was estimated to have experienced population expansions at about 0.052-0.013 Mya, after the EGP (0.5-0.17 Mya), and probably before the LGM (0.021-0.017 Mya). Most individuals of this subclade were collected from Hongya County (population 9). Its elevation is low, and the sampling site in this county (about 1800 m) was lower than other sampling sites (2450-3700 m) in subclade 1A. Thus, Hongya could have provided suitable habitats for species during the glacial period. During the population expansions, many individuals from nearby areas might have migrated into Hongya County. At the same time, some individuals in these regions could have dispersed into nearby areas. This expansion scenario could explain why several individuals from the Jiajin Mountain (population 6) appeared in subclade 1A.
The expansion of the subclade 1B was inferred to have happened at about 0.014-0.004 Mya, right after the LGM (0.021-0.017 Mya). Most individuals of subclade 1B were collected from Gongga Mountain (population 15). Gongga Mountain located at the eastern margin of the Tibetan Plateau, was an important glaciation centre in the eastern Tibetan Plateau. Here, remains of Quaternary glaciers are still widespread. Glacial accumulation and erosion landforms coexisted, which indicated this region has been glaciated many times [20]. The range and intensity of the Quaternary glaciers in this region were very different in different periods. All contemporary glaciers are located at high elevations (above 3000 m). During the early Pleistocene, no glacier developed in the GMR because the relative low elevation at that time. Glaciers developed during the mid-Pleistocene, but most of them occurred at valley heads or high altitudes [20]. The GMR experienced the most extensive glaciation at about 0.277 Mya [20]. As a result, many small areas in this region were suitable for species to inhabit during different periods. Additionally, many ancient and endemic species inhabited in this region [21,23] (e.g. Flora: Taxus chinensis, Kingdonia uniflora, Cordyceps sinensis, Cercidiphyllum japonicum; Fauna: Ailuropoda melanoleuca, Panthera uncia, Pantholops hodysoni, Rhinopithecus roxellanae), leading to the conclusion that the GMR acted as a refuge during the cold periods [21,22].
Another character of Gongga Mountain is its elevation is very large, ranging from less than 1000 m to more than 7000 m, so many types of micro-environments are available [20,21,73]. Apodemus draco has been found at elevations from 800 m to 4000 m, and it occupies different microhabitats, such as broadleaved forests, mixed broadleaf-conifer forests, scrublands, and meadows [23,74]. Therefore, A. draco could have found available microhabitats in high elevations or simply retreated to lower elevations. Consequently, based on the characteristics of A. draco and the geologic history of the GMR, A. draco could have survived in many parts of the GMR even during the EGP and LGM. This conclusion could explain why individuals from GMR fell into different evolutionary clades (Fig. 2).
Two recent studies found the similar pattern in the blood pheasant (I. cruentus) [73], Sichuan field mouse (A. latronum) and Chinese scrub vole (Neodon irene) [22]. But the Sichuan field mouse and Chinese scrub vole had to find suitable places at high elevations in the SEMTP instead of shifting elevation, because they are only found in high elevations [22,23,74]. Therefore, since there were two refugia for A. draco in clade 1 (Hongya county and GMR) in the central region, the discontinuous gene flow among different sampling regions could form the current observed pattern at clade 1.

Conclusion
The present study demonstrated that the genetic differentiation of the South China field mouse's population in the SEMTP not only was influenced by the complex geologic history and unique topography of the SEMTP, but also was influenced by the major glaciations in the Pleistocene. As a consequence, we hypothesize that its evolutionary history included two stages. First, an initial divergence would have been shaped by the uplift events of the Tibetan Plateau, and high genetic diversity was maintained by environmental and topographic complexity. Then, major glaciations in the Pleistocene added further complexity to its demographic history and genetic structure. The observed genetic differentiation pattern in this work is consistent with that of other recent studies focused on the eastern Tibetan Plateau [8,19,22,72]. Our work, together with the studies from Fan et al. [22] and Zhan et al. [72], indicate that past geologic and climatic events on the eastern Tibetan Plateau could have similar effects on the regional historical demography of many animals in the region.