Phylogeographic Patterns of Microtus fortis (Arvicolinae: Rodentia) in China Based on Mitochondrial DNA Sequences

Jun Gao1,2, Liang-Liang Yue3, Xianhuan Jiang1, Liju Ni4, Muhammad Aqeel Ashraf5, Yuxun Zhou,1 Kai Li1,* and Junhua Xiao1,*

1College of Chemistry, Chemical Engineering, and Biotechnology, Donghua University, Shanghai, 201620, China

2Institutes of Animal Husbandry and Veterinary Sciences Shanghai Academy of Agricultural Sciences, Shanghai, 201106, China

3National Plateau Wetland Research Center, Southwest Forestry University, Kunming, 650224, China

4Shanghai Laboratory Animals Research Center, Shanghai, 201203, China

5Faculty of Science and Natural Resources, University Malaysia Sabah 88400 Kota Kinabalu Sabah, Malaysia

Jun Gao and Liang-Liang Yue contributed equally to this work.

ABSTRACT

To examine the phylogeographic relationships of Microtus fortis in China, we investigated 84 individuals collected from five populations. The mitochondrial cytochrome b gene (cyt b) and control region (CR) were sequenced and 49 haplotypes were observed. No shared haplotype was found among different geographic populations. High Fst values among the populations suggested that fragmentation of habitat has resulted in genetically distinct populations. The trees, inferred from maximum likelihood and Bayesian phylogenetic analysis, highly supported all the M. fortis individuals clustering into one monophyletic lineage. Three main clades are recovered within M. fortis: (1) North group; (2) South group; and (3) GX group. The North Group distributed on the north side of Qinling Mountains-Huaihe river line as well as the South Group was on the south. It suggests this geographic barrier played an important role in differentiation of M. fortis in China. Furthermore, the samples all from Southwest China in the GX group may be an example of ‘refuge within refugia’ in glacial period. According to our molecular clock analysis, the main clades of M. fortis divergence and separated time at around 0.77±0.64 million years ago (Mya) located in the Penultimate Glaciation. Divergences within the three clades taken place during the interglacial period between the Penultimate Glaciation and the Last Glaciation. Bayesian skyline plot indicates the effective population size of M. fortis had been experiencing a downward trend in the past decades, which may due to the habitats loss and environmental degradation.

Article Information

Received 02 November 2016

Revised 10 December 2016

Accepted 13 December 2016

Available online 12 June 2017

Authors’ Contribution

JG, KL and JX conceived and designed the study. JG, LN and YZ collected the samples. XJ and JG extracted DNA and performed mtDNA sequencing. JG, L-LY and KL analyzed the data. JG, KL and MAA wrote the article.

Key words

Microtus fortis, cytochrome b gene, D-loop on mitochondrial DNA, Phylogeographic patterns.

DOI: http://dx.doi.org/10.17582/journal.pjz/2017.49.4.1185.1195

* Corresponding author: likai@dhu.edu.cn;

xiaojunhua@dhu.edu.cn

0030-9923/2017/0004-1185 $ 9.00/0

INTRODUCTION

The genus Microtus comprises more than 60 species distributing throughout the Palearctic and Holarctic (Chaline et al., 1999). The karyotypes of contemporary Microtus species range from 17 (2n) to 64 (2n) (Modi, 1987). Their morphological differences are unremarkable (Triant and DeWoody, 2006), which lead to complicated species groups in this genus (Haring et al., 2011). In the past decade, phylogenetic and phylogeographic analysis have been extensively used to resolve relationship within this genus based on mtDNA and nucDNA markers (Brunhoff et al., 2003; Conroy and Cook, 2000; Bannikova et al., 2010; Jaarola et al., 2004). However, most analysis were carried out on West Palearctic species, while only a few studies focused on East Palearctic Microtus species (Haring et al., 2011).

Microtus fortis (Buchner 1889) is mainly distributed in China and part of Russia, Mongolia and North Korea, neighboring to North China. The vole inhabits in the areas with low elevation and everglade was its preferred habitat (Hu et al., 2006). According to morphological characteristics description, the M. fortis has been classified into five subspecies (Huang et al., 1995; Luo et al., 2000) (1) M. f. pelliceus (Thomas 1911), (2) M. f. dolicocephalus (Mori 1930) (3) M. f. fortis (Buchner 1889), (4) M. f. calamorum (Thomas 1911), and (5) M. f. fujianensis (Hong, 1981). The first three subspecies are distributed in northeastern China on the north side of Qinling Mountains-Huaihe river line, whereas the latter two species on the south side. However, the phylogenetic and phylogeographical relationship of these subspecies remain unclear as well as the subspecies classification of M. f. dolicocephalus, M. f. fujianensis, and M. f. pelliceus is still in debate (Ma, 1986; Tan, 1992).

In this study, 84 individuals were trapped in six different geographical location covering the main habit of above five subspecies. Based on the cytochrome b gene (cyt b) and control region sequence (CR), we examined subspecies taxonomy, genetic diversity, phylogeo-graphical pattern and investigate the population demographic history of M. fortis in China. It also reveals the divergence history of M. fortis and the impact of Quaternary climate change, geographic barriers, and human activities on the M. fortis geographic populations.

MATERIALS AND METHODS

Sampling and DNA extraction

Eighty four individuals from 6 provinces of China (Table I, Fig. 1) were trapped using sunflower seeds baited clamps. The sample tissue (tail and legs) were immediately preserved in absolute ethanol, and finally stored at -20°C. Genome DNA was extracted from the tail tissue following the standard phenol-chloroform method (Sambrook et al., 1989).

PCR and sequencing

Complete cyt b gene (1143bp) and CR sequence (910-916bp) were amplified by two pairs of specific primers: CYTB1 (5’ - CCTGAACACCCGCTAACAAT - 3’) and CYTB2 (5’ - TGGAGGGGTAGTCCTTCCTT - 3’); CR1 (5’ - AAGGAAGGACTACCCCTCCA - 3’) and CR2 (5’ -ATAAGGCCAGGACCAAACCT - 3’). PCR reactions were carried out in 25 µL volume and cycling conditions consisted of 95°C for 5 min, 35 cycles of 95°C for 1 min, 58°C for 1 min, 72°C for 1 min, and a final extension of 72°C or 10 min. PCR products were sequenced in both directions performed by Biosune Inc. (Shanghai, China) using ABI 3730 systems.

Phylogenetic reconstruction and molecular divergence analysis

All the sequences ends were trimmed, manually checked and assembled into contigs using the SeqmanTM in the Lasergene v7.1 (DNAStar Inc., Madison, WI). The sequences of the two loci were combined (cyt b+CR) and aligned by Clustal W in MEGA5.2 (Tamura et al., 2011). Additionally, a total of 12 cyt b sequences of M. fortis and 8 outgroup Microtus species (Table II) were downloaded from GenBank for phylogenetic reconstruction. The program DnaSP v5 (Librado and Rozas, 2009) was used to extract haplotypes from combined sequence alignments. Molecular substitution saturation status were tested using the program DAMBE version 5.6.7 (Xia, 2013) to confirm that the status of the data sets were not over mutated and suitable for phylogenetic analysis.

We used the program MrModeltest version 2.3 (Nylander, 2004) to infer the most suitable evolutionary model for the dataset. The best model suggested by the Akaike information criterion (AIC) was used to reconstruct phylogenetic tree in Bayesian inference (BI) and Maximum Likelihood (ML) analysis.

Table I.- Sampling locations of Microtus fortis in this article.

Group	Location	Number (Year)	LAT - LONG
HLJ	Fuyuan County in Heilongjiang Province	3 (2012)	47.5N 133.5E
NX	Linwu City in NingXia Province	20 (2012)	38.2N 106.2E
HN	Datonghu area in Hunan province	19 (2007)	29.1N 112.5E
FJ	Jianyang City in Fujian Province	20 (2012)	27.1N 117.2E
GX	Linchuan County in Guangxi Province	20 (2012)	25.2N 110.1E
JL	Jilin City in Jilin Province	2 (2013)	43.5N 125.6E

Table II.- The other cyt b sequences information used in phylogenetic analysis.

Sample sequence	Location	Number	Accession Number
M. fortis (YJ-2/3/13/15/78)	Yuanjiang city in Hunan province	5	EU870635/EU870632/EU870637/ EU870636/EU870633
M. fortis (SY-135/136)	Shenyang city in Liaoning Province	2	EU870639/EU870638
M. fortis-682	JiLin Province in China	1	EU126809
M. fortis	Mongolia: Selenge Aymak, upper Ero River 49.083°N, 107.283°E	1	FJ986308
M. fortis	Russia: Chita Region, 50.410 °N,118.110°E	1	FJ986307
M. fortis	Russia: Buryatia, Dzhidinskiy District, Dodo-Ichetuy 50.583°N, 105.500° E	1	FJ986306
M. fortis	Korea kyonggi Prov. Dong Du Chun	1	AF163894
M. middendorffii	Asia	1	FJ986314/FJ986315/AF163898
M. monglicus	Asia	1	FJ986304/FJ986305
M. maximowiczii	Asia	1	FJ986303/ FJ986311
M. oeconomus	Asia	1	AY220018/ AY220028
M. montebelli	Japan, Asia	1	AF163900
M. Kikuchii	Taiwan, Asia	1	AF163896/ AF348082/NC003041
Chionomys nivalis	Syria	1	AY513849
Chionomys roberti	Turkey	1	AY513850

*The sampling locations of twelve downloaded ingroup sequences (Bold) labeled in the Figure 1.

Bayesian inferences (BI) analysis were performed using Mrbayes 3.1.2 (Huelsenbeck et al., 2001). Maximum Likelihood (ML) analysis were performed in MEGA5.2 (Tamura et al., 2011) and statistical confidence was assessed through the bootstrap analysis with 500 replications.

Molecular clock and time calibrations

We used BEAST v1.7.5 (Drummond and Rambaut, 2007) to simulate relaxed molecular clock progress to dating the divergence time for each clade. The analysis started with a randomly generated starting tree, uncorrelated lognormal relaxed molecular clock, and the program’s default prior distributions of model parameters. The MCMC analysis consisted of 50 million generations, sampled every 100 generations. The analysis was assessed using Tracer v1.5 (Rambaut and Drummond 2007). Posterior probabilities (PP) ≥0.95 are considered to be strongly or significantly supported (Ronquist and Huelsenbeck, 2003). Two divergence dates were used as the calibration point (subgenus Alexandromys/ Pallasiimus divergence time =1.19±0.19 Mya and Alexandromys basal radiation time=0.84±0.14 Mya, according to estimation by Bannikova et al. (2010). The molecular clock calibrations were set as lognormal distributions in analyses under a relaxed molecular clock model.

Population polymorphic and dynamic analysis

The haplotype diversity (Hd), nucleotide diversity (π) and the pairwise FST among populations were calculated in DnaSP v5 (Librado and Rozas, 2009). Number of phylo-groups was inferred from the genetic distances and geographic distances by the program SAMOVA (Dupanloup et al., 2002). Molecular variance test (AMVOA) in Arlequin version 3.5 (Excoffier and Lischer, 2010). Fu’s test (Fu, 1997), Ramos-Onsins and Rozas’s test (Ramos-Onsins and Rozas, 2002) and mismatch distributions (Harpending, 1994) are calculated also using the software DnaSP v5 (Librado and Rozas, 2009). In mismatch and neutrality tests, we combined the HLJ and JL populations into a unique population, because the two populations contain only five individuals in total. A mantel test was simulated to infer the isolation by distance hypothesis. This test was run using the program Arlequin employing pairwise FST matrix calculated from its own and the geographic distances matrix from the GenAlex 6 (Peakall and Smouse, 2006) for the populations. Bayesian skyline plot (BSP) (Drummond et al., 2005) was simulated to show the dynamics of the maternal effective population size through time using the program BEAST v1.7.5 (Drummond and Rambaut, 2007). In the BSP progress, a tip date encoded by the sample time were used for time calibration. The MCMC program started by a randomly generated tree and simulated for 50 million generations with a burn-in periods of 12000. The effective sample sizes (ESS) defined by the program authors for the posterior and prior distributions were used to evaluate the MCMC runs.

RESULTS

Phylogenetic analysis

Two fragments of mtDNA were obtained from all 84 samples. The total lengths of the cyt b gene and the CR fragment alignments are 1143bp and 910-916bp, respectively. All the sequences of M. fortis were deposited in NCBI under accession number of KJ081871-KJ081954 for cyt b gene and KJ207290-KJ207373 for CR fragment. 49 haplotypes (gaps were considered as substitution) were extracted from the all 84 individuals’ combined fragments (cyt b + CR) and no shared haplotypes were found among different populations. The overall Hd and π of the populations were 0.969 and 0.021, respectively.

The saturation tests showed that the observed impulse indicator saturation (IIs) are significantly lower, indicated that the data set was not mutation saturated and suitable for phylogenetic analysis. The best model suggested by the AIC was the GTR+I+G model (-lnL= 7220.31). A monophyletic relationship within the haplotypes of M. fortis was inferred by both the Bayesian inference and the ML methods. Furthermore, a pattern of three phylo-groups was found according to the phylogenetic trees with posterior probabilities of 1 in Bayesian tree (Fig. 2), and with bootstrap value over 65 in ML tree (Fig. 3).

Table III.- Molecular variance (AMOVA) of genetic variation for the 3 groups (5 populations) of Microtus fortis.

Source of variation	Variance components (%)	Percentage of variation (%)
Among groups	15.73	51.68
Among populations within groups	10.42	34.24
Within populations	4.28	14.06
Total	30.45

Average F-statistics over all loci, Fixation Indices: FSC=0.709; FST=0.859; FCT=0.517.

Population structure

SAMOVA suggested that the FCT value first got to the high plateau when K=3, indicated that the populations could be clustered into three phylo-groups according to their pairwise genetic divergences (Fig. 4). The first group (North group) comprises the individuals of M. fortis from North China include northeast area (Heilongjiang, Jilin, and Liaoning), north central area (Ningxia) and sequences of Russia, Mongolia and Korea samples. The second group (South group) consists of the samples from south china (Hunan and Fujian). The third group (GX group) formed by the samples of M. fortis from Guangxi in Southwest China. The majority of haplotypes of M. fortis are clustered consistent with their geographical location, except for several haplotypes (Hap 31/34/40/43 in BI tree, Hap33/34/41 in ML tree ) of FJ population forming cluster to the Hunan population in South group of M. fortis (Figs. 2, 3).

The lowest pairwise FST value (FST=0.357<0.5) was observed between HLJ+JL population with NX population and the highest pairwise FST value (FST=0.827) was observed between the NX population with the HN population. The other FST value among these 5 populations all above the 0.5 suggested these geography populations are strongly differentiated from each other.

In molecular variance analysis (AMOVA) (Table III), the most proportion of variances (51.69%) were found among groups, followed by the variances among populations (within groups) and then within the populations. A Mantel test indicates that subdivision within M. fortis dose not fit an isolation-by-distance model with a none significant P value (P=0.63>0.05).

Molecular divergence time

Divergence times of M. fortis based on the Bayesian relaxed molecular dating estimation (Fig. 5) indicates M. fortis differentiated from subgenus Alexandromys at about 0.93 Mya (95% CI=0.52-1.35) and the basal radiation time of M. fortis at about 0.77 Mya (95% CI=0.42-1.11). The north group separated from the south group of M. fortis at about 0.64 Mya (95% CI=0.29-0.99). The GX clade was at the basal position of the phylogenetic tree and it is possible to represents an ancient haplotype of M. fortis.

Table IV.- Neutrality test and Mismatch distribution analyses of M. fortis group.

	Clade/Group (sample size)	Neutrality test				Mismatch distribution analyses
	Clade/Group (sample size)	Fu’s Fs	PFS	R2	PR2	SSD	PSSD	rg	Prg
3 clade/group suggested by phylogenetic tree	GX group (20)	2.539	0.870	0.084	0.083	0.034	0.040	0.069	0.040
	South group (39)	1.152	0.704	0.162	0.954	0.026	0.020	0.034	0.010
	North group (25)	1.497	0.748	0.134	0.560	0.056	0.240	0.130	0.080
5 geograp7hic sampling population	HN population(19) in South group	-0.500	0.410	0.103	0.152	0.047	0.541	0.148	0.339
	FJ population(20) in South group	-0.133	0.458	0.184	0.948	0.072	0.016	0.047	0.134
	NX population(20) in north group	-0.573	0.408	0.085	0.032	0.017	0.653	0.053	0.695
	HLJ+JL population (5) in north group	0.240	0.344	0.249	0.693	0.105	0.301	0.180	0.754
	GX population(20) in GX group	-1.816	0.221	0.134	0.570	0.009	0.704	0.013	0.877

Population historical demography of M. fortis

The mismatch distribution analysis (Table IV) showed no smoothness of the distributions of pairwise genetic differences. Moreover, Fu’s and R2 test suggested a neutral status of the M. fortis populations. No sudden expansion population dynamics can be inferred among all the range. Bayesian skyline plot (BSP) analysis in cyt b, CR and combined data all suggested that the effective population size of M. fortis has been declining since 30-40 years ago (Fig. 6).

DISCUSSION

Phylogeographic structure and taxonomic implication

The phylogenetic trees (Figs. 2, 3) confirmed that the M. fortis belonged to the subgenus Alexandromys within the “Asian lineage” of genus Microtus. The phylogenetic trees also confirmed that M. fortis belongs to the subgenus Alexandromys within the “Asian lineage” of Microtus, which was extensively described (Conroy and Cook, 2000; Bannikova et al., 2010; Jaarola et al., 2004). Three distinct monophyletic clades including North group, South group and GX group of M. fortis were inferred from our phylogenetic trees in China with highly supports. Sequences downloaded from GenBank are also clustered to its clade and also coincided with their geographical distributions. It suggested that the differentiation among the populations could relate to geographical isolations.

The subspecies classification of M. fortis has long been controversial in China. In the former researches, five subspecies of M. fortis distributed in China, as listed, M. f. pelliceus (in northeast) M. f. dolicocepalus (in JL and Liaoning province of northeast) M. f. fortis (in Ningxia province), M. f. calamorum (in Hunan province), M. f. fujianensis (in Fujian province). The GX population was first reported by Liang and Huang (1997) and in our original idea it might be generated by diffusion or migration from Hunan population (M. f. calamorum subspecies). Surprisingly, the GX group is distinctly diverged from other populations, suggested GX group did not belong to the one of five exist known subspecies. In the North group, the NX clade is a monophyletic clade, supporting the former morphology study. HLJ, JL and other north east sequences are polyphyletic. However, due to limited number of northeast samples, it is not easy clarify the relationship among the northeast individuals. The situation is more or less complex in the south group. The FJ population is paraphyletic population, because several haplotypes (Figs. 2, 3) from FJ population are clustered to the HN sub-clade and well supported by statistical tests, the FJ population could has been mixed in the past by individuals from the HN populations.

Divergence time and historical demography of M. fortis

There is a wide variation of mutation rate in mammals. The overall rate of the cytochrome b gene in Microtus has been estimated to be 3-7 times higher (Conroy and Cook, 2000; Brunhoff et al., 2003) than the widely used mtDNA standard rate of 2% sequence divergence per million years (Avise, 1998). It was also reported the divergence rates from more than 30% for most recent splits down to 12-14% per Myr for basal-most nodes (Bannikova et al., 2010). Thus, it is unsuitable to estimate the divergence time of M. fortis using a strict molecular clock.

With a Bayesian method under a relax molecular clock model, it is inferred M. fortis basal radiation time at around 0.77 Mya, around the starting of the glacial period MIS 12-16. The separate time of South-North group is about 0.64Mya within the glacial stage of MIS 12-16. The three main clades radiated from 0.5 to 0.37Mya, located at the interglacial stage between the MIS 12-16 and Penultimate Glaciation (Ehlers and Gibbard, 2007; Schäfer et al., 2002; Zheng et al., 2002; Lehmkuhl and Owen, 2005; Owen and Benn, 2005) (also see the visible graph edited by Yue and Sun, 2014). Although there is no evidence to support the unique glacier cap in China, it’s still safe to infer that cold/dry climate and wide spread mountain glacier blocked gene flow for the species between populations and triggered speciation progress showing in the phylogenetic tree. After the glacier period, a warm and wet climate appears, rivers and Lake Flood again, wetland developed. The suitable habits might largely contribute to the divergence within the M. fortis groups. Although, there is no direct palynology evidence for the development of the wetland, it is still suggested for forests develop by palynology researches either in mainland China (Liu, 1988; Wang, 1992; Axelrod et al., 1996; Manchester et al., 2009) or the world wide (Tiffney and Manchester, 2001; Tiffney, 2008). At least, a warm and wet condition can be inferred, which is suitable for this species.

Population structure and phylogeography of M. fortis

Significant population structure was detected among the populations of the species. However, isolation by distance test was failed in the Mantel Test. Individuals from North China was closely related and clustered into one monophyletic clade, even our sample locations in the north China were far away from each other. The lowest FST was also observed between the NX population with northeast population (HLJ+JL). Furthermore, because no evidence of sudden expansion was find within the northern group, the population of NX and northeast China were more likely experienced in situ processes during the last glaciation rather than bottleneck or long distance migration.

The in situ patterns also occurred in the GX population. The GX clade was at the basal position of the phylogenetic tree within M. fortis groups. Therefore, it seemed that the mtDNA of GX group individuals did not contribute to the postglacial recolonization of populations North. Thus, it was possibly isolated in Guangxi province in Southwest China. This might be an example of ‘refuge within refugia’ which described for some species of the Iberian Peninsula or Calabria (Gómez and Lunt, 2007; Grill et al., 2009). It was supposed that there might be some important refuges in the southwest region (e.g., west Yunnan mountain gorge region, southern mountainous region of Guangdong, Guangxi, Hunan and Jiangxi), where were less affected by the Quaternary glaciers in China (Wang and Liu, 1994).

The North and South groups were located on either side of Qinling Mountains-Huaihe river line. This east-west trending boundary line (32~34.5°N, Fig. 1) is the traditional climate divide between arid North china and humid South china. This important geography line approximates coincided with the 0 January isotherm and the 800 mm isohyet in China. Hence, this geographic barrier may lead to fragmentation of M. fortis populations. Similar phenomenon was also found in Holarctic phylogeography of the root vole (Microtus oeconomus), of which central Asian and the North European lineages are separated by the Ural Mountains (Brunhoff et al., 2003). In China, the Mekong-Salween Divide also blocked gene flow for several plant species (Li et al., 2011; Qiu et al., 2011).

It is considered that natural and human activities are two main factors strongly formed the distributions of organisms. Bayesian skyline plot (BSP) analysis in our study suggests the effective population size of M. fortis declined since 30-40 years before (Fig. 6) which coincided with our sightseeing mentioned above. In the past 50 years, China experienced a dramatic development, which triggered huge change of landscape. Earth management was more likely to fill the lakes to produce new farmland and manufactory industry land. The habitats of M. fortis are mainly in coastal rivers and lakes, Carex meadows, marshes. Since the 1950 to 1997, Yangtze River wetlands have been reclaimed of 7.85×104 square hectometer, the equivalent area of land area of 12.39%. Nearly 1000 natural lakes were disappeared due to the reclamation in China in this period. The wetlands also drastically reduced in recent decades (Lei and Zhang, 2005). Wetland ecological function decreased significantly reducing biodiversity, deterioration of the ecological environment. Unfortunately, the vole species was considered as harmful species for farm activity and water conservancy facilities. As the results, the reduction of this species was considerable.

CONCLUSION

Climate change and human activities are the factors affecting the distribution of organisms. Environmental changes profoundly shaped the current distributions and genetic structures of many plant and animal species in Northern Hemisphere. There was a large land area in China with highly diverse animal and plant resources which provides rich resources for studying molecular evolution. Our study has shed light on the rodent M. fortis phylogeography pattern. Three main groups of M. fortis were recognized in distribution by our mtDNA sequence data. It provided molecular genetic evidence will helpful to solve the subspecies classification of M. fortis. Fragmented distribution and high level of genetic variation among populations suggest the Quaternary climate change, strong geographic barriers and human activities were the important influence factors. Our findings also suggest the some glacial refuge exist in southwest China. The phylogeography study of species in these areas is conducive to reveal the inter- and intraspecific divergence history in China.

ACKNOWLEDGMENTS

The authors are grateful to Huanan Huang, Heng Zhan, Jinxiang Wang, Yuqin Yang, and Jian Huang for their help with sample collection. This work was supported by the Key Project of Science and Technology Commission of Shanghai Municipality (No. 13140900300) and the National Science Foundation of China (No. 31171199). The Fundamental Research Funds for the Central Universities (2232013A3-06), and the DHU Distinguished Young Professor Program (B201308).

Statement of conflict of interest

Authors have declared no conflict of interest.

REFERENCES

Avise, J.C., 1998. The history and purview of phylogeography: a personal reflection. Mol. Ecol., 7: 371-379. https://doi.org/10.1046/j.1365-294x.1998.00391.x

Axelrod, D.I., Al-Shehbaz, I. and Raven, P.H., 1996. History of the modern flora of China. In: Floristic characteristics and diversity of East Asian plants (eds. A. Zhang and S. Wu). Springer, New York, pp. 43-55.

Bannikova, A.A., Lebedev, V.S., Lissovsky, A.A., Matrosova, V., Abramson, N.I., Obolenskaya, E.V. and Tesakov, A.S., 2010. Molecular phylogeny and evolution of the Asian lineage of vole genus Microtus (Rodentia: Arvicolinae) inferred from mitochondrial cytochrome b sequence. Biol. J. Linn. Soc., 99: 595-613. https://doi.org/10.1111/j.1095-8312.2009.01378.x

Brunhoff, C., Galbreath, K.E., Fedorov, V.B., Cook, J.A. and Jaarola, M., 2003. Holarctic phylogeography of the root vole (I): implications for late Quaternary biogeography of high latitudes. Mol. Ecol., 12: 957-968. https://doi.org/10.1046/j.1365-294X.2003.01796.x

Chaline, J., Brunet-Lecomte, P., Montuire, S., Viriot, L. and Courant, F., 1999. Anatomy of the arvicoline radiation (Rodentia): palaeogeographica history and evolutionary data. Ann. Zool. Fenn., 36: 239-267.

Conroy, C.J. and Cook, J.A., 2000. Molecular systematics of a holarctic rodent (Microtus: Muridae). J. Mammal., 81: 344-359. https://doi.org/10.1093/jmammal/81.2.344

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

Drummond, A.J., Rambaut, A., Shapiro, B. and Pybus, O.G., 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol. Biol. Evol., 22: 1185-1192. https://doi.org/10.1093/molbev/msi103

Dupanloup, I., Schneider, S. and Excoffier, L., 2002. A simulated annealing approach to define the genetic structure of populations. Mol. Ecol., 11: 2571–2581. https://doi.org/10.1046/j.1365-294X.2002.01650.x