Pleistocene glacial cycle effects on the phylogeography of the Chinese endemic bat species, Myotis davidii
© You et al; licensee BioMed Central Ltd. 2010
Received: 30 December 2009
Accepted: 10 July 2010
Published: 10 July 2010
Global climatic oscillations, glaciation cycles and the unique geographic topology of China have profoundly influenced species population distributions. In most species, contemporary distributions of populations cannot be fully understood, except in a historical context. Complex patterns of Pleistocene glaciations, as well as other physiographic changes have influenced the distribution of bat species in China. Until this study, there had been no phylogeographical research on Myotis davidii, an endemic Chinese bat. We used a combination of nuclear and mitochondrial DNA markers to investigate genetic diversity, population structure, and the demographic history of M. davidii. In particular, we compared patterns of genetic variation to glacial oscillations, topography, and environmental variation during the Pleistocene in an effort to explain current distributions in light of these historical processes.
M. davidii comprises three lineages (MEP, SWP and SH) based on the results of molecular variance analysis (AMOVA) and phylogenetic analyses. The results of a STRUCTURE analysis reveal multi-hierarchical population structure in M. davidii. Nuclear and mitochondrial genetic markers reveal different levels of gene flow among populations. In the case of mtDNA, populations adhere to an isolation-by-distance model, whereas the individual assignment test reveals considerable gene flow between populations. MDIV analysis indicate that the split of the MEP and SWP/SH lineages, and from the SWP and SH lineages were at 201 ka BP and 158 ka BP, respectively. The results of a mismatch distribution analysis and neutrality tests indicate a population expansion event at 79.17 ka BP and 69.12 ka BP in MEP and SWP, respectively.
The complex demographic history, discontinuous extant distribution of haplotypes, and multiple-hierarchy population structure of M. davidii appear associated with climatic oscillations, topography and eco-environmental variation of China. Additionally, the three regions are genetically differentiated from one another in the entire sample set. The degree of genetic differentiation, based on the analysis of mtDNA and nDNA, suggests a male-mediated gene flow among populations. Refuges were in the MEP, SH and the lower elevations of SWP regions. This study also provides insights for conservation management units (MEP, SWP and SH).
Global climatic oscillations and associated glaciation cycles have profoundly influenced species population distributions . Pleistocene climate oscillations have contributed to the genetic structure of current species  and have played an important role in initiating intra-specific divergences of North American animal taxa . However, climate condition in China during glacial periods appears to have been different from those of the rest of the world and bring about a series of eco-environmental effects, primarily due to the continuing uplift of the Qinghai-Tibet Plateau . The uplift has heavily influenced the climatic and environment changes, and formed a complicated topography (a set of three ladders with the highest point to the west and the lowest to the east). The uplift has also had a decisive effect on the formation of the East-Asia monsoon, which increased the climatic differences between the glacial and interglacial period and influenced the extent of the glacial advance in China . Additionally, climatic change induces changes in vegetation distribution, which also influences the distribution of animal populations through food chain and habitat modifications . Thus, the distribution and change of forest vegetation may have influenced the migration, expansion and even extinction of animal populations which differed considerably from other regions of the world at the same latitude .
David's myotis (Myotis davidii) is endemic to China, occurring in a small province of middle and northern China . This species is considered of least concern by IUCN . The measurements of the bat specimen have on average 31.7 mm forearm length , 237.21 mm wingspan and weight 4.31 g (unpublished data). The wing morphology of this bat species is characterized by low wing loading and low aspect ratio, indicating that M. davidii has typically good manoeuvrability . They prey on flying insects and inhabit rock or karst cavities . Generally, cavities inhabited by this species are located in forests with high levels of insect diversity (unpublished results). However, there are no reports about migration, hibernation or ecological preferences. Previous research  on animal fossils of China has shown that the survival of bats depended mainly on forest cover. Furthermore, bat numbers or activity changes can be related to climate change . Thus, bats are a good model for examining historical processes and distributional patterns following climatic oscillations.
In China, studies on the geographic distribution pattern of vertebrates, such as bats , birds  and amphibians , have described evolutionary histories following the global Pleistocene climatic oscillations and reflected differences with Europe in terms of evolution history in the context of complex regional scenarios . Colonization events, dispersal patterns and migratory behaviours play a key role in determining a species' geographical distribution range and demographic history [15, 16]. In addition, both historical events and ecological factors shape extant genetic diversity and population structure of animals . Usually, analysis of genetic markers can be useful to successfully differentiate fine-scale structures in natural populations [16, 18].
Different markers have their own unique genealogy, which may indicate concordance or divergence from the species' history. Mitochondrial DNA only reveals a species' maternal demographic history. However, microsatellite-derived data is influenced by sex-biased dispersal because of the character of biparental inheritance. Thus, the combination of different markers may actually provide complementary information to test the population demographic history  and extant population structure . We therefore chose to include both nuclear DNA (nDNA) and mitochondrial DNA (mtDNA) analyses in our study of the phylogeography of M. davidii, for which genetic diversity, demographic history, and population structure remain unknown. Our main aims were to identify genetic diversity, describe macrogeographical population genetic patterns, and investigate the population demographic history within a context of climate oscillations, complex topology and eco-environmental variation since the mid-late Pleistocene. We inferred the historical causation of the extant population structure to better understand its distribution and population status. Additionally, we wanted to identify possible conservation management units and ultimately use information from the study to provide some advice on the further protection for M. davidii.
Genetic variability within the studied Myotis davidii populations based on mtDNA and nDNA data.
GenBank Accession No.
N31° 33' E118° 05'
0.950 ± 0.100
0.030 ± 0.019
0.50 ± 0.12
0.68 ± 0.13
N29° 47' E118° 10'
0.952 ± 0.096
0.022 ± 0.014
0.49 ± 0.22
0.63 ± 0.23
N31° 22' E119° 48'
0.949 ± 0.044
0.050 ± 0.027
0.59 ± 0.25
0.60 ± 0.23
N30° 06' E120° 02'
0.600 ± 0.131
0.028 ± 0.016
0.54 ± 0.37
0.69 ± 0.17
N26° 36' E114° 12'
0.956 ± 0.054
0.027 ± 0.015
0.52 ± 0.11
0.71 ± 0.17
N30° 58' E108° 08'
0.53 ± 0.01
0.56 ± 0.03
N29° 16' E107° 50'
0.43 ± 0.26
0.58 ± 0.18
N28° 17' E109° 39'
0.810 ± 0.172
0.015 ± 0.009
0.58 ± 0.22
0.71 ± 0.29
N27° 59' E107° 11'
0.600 ± 0.175
0.007 ± 0.005
0.70 ± 0.42
0.74 ± 0.10
N28° 23' E106° 25'
0.997 ± 0.127
0.032 ± 0.020
0. 63 ± 0.22
0.73 ± 0.24
N22° 36' E100° 43'
0.500 ± 0.250
0.004 ± 0.003
0.69 ± 0.34
0.75 ± 0.13
N22° 46' E100° 05'
0.57 ± 0.29
0.72 ± 0.20
N25° 03' E102° 42'
0.55 ± 0.26
0.71 ± 0.12
N26° 28' E100° 50'
0.61 ± 0.32
0.70 ± 0.19
N25° 24' E110° 41'
0.730 ± 0.220
0.007 ± 0.004
0.61 ± 0.36
0.69 ± 0.29
N24° 46' E113° 34'
0.930 ± 0.165
0.012 ± 0.008
0.57 ± 0.19
0.59 ± 0.26
N24° 44' E113° 31'
0.453 ± 0.245
0.025 ± 0.019
0.44 ± 0.21
0.52 ± 0.27
mtDNA and nDNA genetic diversity
Examination of mitochondrial HVI sequences showed length variation (579-909 bp). The R1 repeat (consisting of three to seven 81 bp repeats) within mtDNA HVI region which responsible for the length variation of HVI region, were excluded from the analysis. After exclusion of the tandem repeats, the 340 bp sequence of the HVI region revealed 53 different haplotypes, defined by variation at 159 polymorphic sites. Thirty-eight of these polymorphisms were the result of indels, and 72 polymorphic sites were parsimony-informative. Haplotype diversity (h) averaged over all populations 0.975 ± 0.002. Haplotype diversity (h) was uniformly high and was the highest in the MEP region (Table 1). Most haplotypes were unique (36 haplotypes, 28.57% of the individuals). Other haplotypes were shared within populations (17 haplotypes, 69.84% of the individuals), among populations (5 haplotypes, 38.89% of the individuals) or among regions (2 haplotypes, 21.43% of the individuals; AH1 shared with YN1/YN2/YN3 and JS shared with YN1/YN4) (Additional file 1). Nucleotide diversity (π) averaged over all populations 0.062 ± 0.007. Patterns of variation in nucleotide diversity across regions were consistent with that of haplotype diversity in that the greatest diversity was found in the MEP and the lowest in the SWP and the SH lineages (Table 1).
Characteristics of eight microsatellite loci for Myotis davidii.
Primer sequence (5'-3')
Fragment size (bp)
HW P value
Phylogenetic and demographic analysis
Demographic history analysis of Myotis davidii
Mismatch distribution analysis
(CI = 95%)
An AMOVA revealed significant genetic variance for all three hierarchical levels examined (among regions, among populations within regions, and within populations) (Additional file 2). The grouping pattern was MEP (composed of AH1, AH2, JS, ZJ, JX, CQ2), SWP (composed of CQ1, GZ1, GZ2, HN, YN1, YN2, YN3, YN4), and SH (composed of GD1, GD2, GX) gave the highest ΦCT value (0.771, P < 0.01). Most of the genetic variance (64.82%) was explained by differences among the three regions (Additional file 2). The subdivisions depicted by the AMOVA are consistent with the macrogeographical structure of China (Fig. 1). In addition, the analysis of AMOVA highlighted a low but significant microsatellite genetic variation (14.45%) among the three regions and a high proportion of genetic of variation (68.98%) within populations (Additional file 2). In mtDNA, the pairwise genetic differences among populations varied from 0.007 to 0.988 (P < 0.01) (Additional file 3), and only 3.13% of pairwise genetic differences were not significant in the whole population. In nDNA, the pairwise genetic difference among populations varied from 0.004 to 0.276 (P < 0.01) (Additional file 3), and 82.14% pairwise genetic differences in the SWP lineage were not significant. Thus the pairwise genetic differences among populations within regions were smaller than among regions in both markers (Additional file 3). These results suggest that the genetic differentiation maximizes among the three regions.
A Mantel test indicated that genetic subdivision within M. davidii fits an isolation-by-distance model (R2 = 0.296, P = 0.031) with respect to mtDNA. However, nuclear markers failed to support isolation-by-distance model (R2 = 0.030, P = 0.895).
Individual assignment test indicated that a greater proportion of individuals of M. davidii (63.49%) were assigned to the populations from which they were sampled. Within the MEP, SWP and SH regions, 44.44%, 23.53% and 19.05% of the individuals, respectively, were assigned to other populations. Thus, results of the individual assignment test also pointed toward higher gene flow among populations within regions, especially within the MEP region. When using individual assignment tests according to the three mtDNA clades, one individual (0.79%) in the SWP (HN) group were assigned to the SH (GD2) group. Two individuals (1.58%) in the MEP (ZJ) group were assigned to the SWP (YN3/YN4) group. Three individuals (2.38%) in the SWP group were assigned to the MEP group, where YN3 and GZ2 originated from ZJ/JS and JX, respectively. Thus, 4.75% of individuals were assigned to other lineages, which indicated a weak gene flow among the three regions.
Genetic structure of populations
Compared with similar study on bat mtDNA , M. davidii showed higher mean nucleotide diversity (mean π = 0.062 ± 0.007). The presence of high h and high π indicated that a sub-lineage formed over long periods of evolutionary time . In addition, a high h and high π indicated high levels of divergence within three regions, which was attributed to secondary contact between previously differentiated allopatric lineages. Like other bat species [18, 19], M. davidii was not detected the continuity gene flow among populations based on mtDNA. In addition, isolation by distance was an important factor for the formation of genetic structure within maternal populations. The pattern of subdivision into subpopulations can be explained by the influence of environmental characteristics. In China, the west-east axis was a direction of the prominent changes in many environmental variables, such as vegetation types, temperature, precipitation, and topography, all of which were due to the continuing uplift of the Qinghai-Tibet Plateau. The dependence of genetic differentiation on a gradient of environmental variables is in agreement with the theoretical model of Doebeli & Dieckmann , showing that processes of evolutionary diversification may lead to sharp geographical differentiation along environmental gradients. The environmental gradients existed for an extended period of time which reflected in mtDNA.
However, the population genetic structure in microsatellite loci was less pronounced than in the case of mtDNA analysis. Bayesian clustering analysis (Fig. 3) and individual assignment test showed frequent migration within regions. And the patterns of migration followed a stepping-stone colonization model. Large and non-significant Fis values also showed frequent transfer in males among populations and led to an unobvious population structure. Patterns of genetic differentiation in two molecular markers were attributed to male dispersal and female philopatry [22, 23]. Male-biased gene flow implied low introgression of mtDNA haplotypes from neighboring populations, and therefore greater structuring in mtDNA as compared with nuclear markers.
Our comparisons showed that the deepest phylogenetic split (SWP/SH vs. MEP, Fig. 2) among the haplotypes corresponded exactly to nDNA-based structure when genotypes were in two groups (K = 2, Fig. 3). The deepest phylogenetic split was similar to the pattern obtained from Rhinolophus ferrumequinum in China . SWP and SH lineages split during stage III of the glacial epoch in the Yun-Gui plateau (154-136 ka BP) , and separated into two geographically defined lineages (Fig. 2), which was in agreement with the results of the nDNA analysis (K = 3, Fig. 3). The population genetic structure was consistent with topographic and geographic characteristics of China. Thus, ecological environmental changes, landscape structure differentiation, and the climate discrepancy are primary determinants of population boundaries and rate of movement among regions [25, 26].
When these regions were considered separately we found clear differences in the phylogeographical signal obtained from the two sets of markers. Multi-locus genotyping data suggest that M. davidii has multiple levels of population structure (Fig. 3), which were not shown by mtDNA. Compared with mtDNA, the concordance between the two markers seemed to break down within the SWP clade. These similarities and discrepancies between the data sets together clarify the history of this species within the SWP region. In the SWP, four sampling locations in Yunnan, formed an independent lineage (Fig. 2), which coincided with the elevation difference of other populations (GZ, CQ and HN; Fig. 1). This result is in agreement with previous studies on birds  and amphibians . However, nDNA analysis indicated that the separation of the YN1-YN2 vs. other populations in the SWP group was maintained following the introduction of a fourth and fifth microsatellite-based cluster, which indicated much more frequent contact between YN1 and YN2 populations than between these two and any other population. Additional discordance at K = 4 indicated postglacial colonization across the whole range of SWP region with evidence of a cline in membership between these clusters. Discordance at K = 6 indicated that CQ1 population might be isolated by geographical barrier, although CQ1 and other locations (HN, GZ1 and GZ2) were all in the same lineage of the phylogenetic tree (Fig. 3). These clusters in YN1-YN2 and CQ1 may represent a discrete, discontinuous jump across a relatively small geographical distance. The gene exchange among M. davidii populations in SWP also may be impeded by the limit of landscape.
In the MEP region, no obvious population structure was found at the higher level of K = 5-7, which indicated a continual population contact with one another. Extensive and overlapping ranges within the MEP region of the Bayesian clustering analysis (Fig. 3) indicated the absence of geographical features in the landscape that would not constitute efficient barriers for bat dispersal. SH also showed no clear population division and indicated a much more stable population structure than that of MEP.
Population demographic reconstructions
Since the Middle and Late Pleistocene, violent glacial-interglacial cycles (at least 24 Dansgaard-Oeschger cycles and several Heinrich events between 115 ka BP and 14 ka BP; on average, each fluctuating glacial cycles persisted for about 1500 years ), complicated topology and eco-environmental changes due to the Qinghai-Tibet Plateau uplift  had important influences on the demographic history of many species [12–15]. With the complexity climatic conditions due to glacial oscillation and Qinghai-Tibet Plateau uplift, Chinese species have experienced a unique evolutionary history .
Based on the divergence time estimated among mtDNA lineages, the demographic history of M. davidii could be traced back to before the Riss glaciation (210-135 ka BP) during the Pleistocene , which was similar to the glacial epoch in middle, southern and southwestern China (223-189 ka BP) . Both markers thus support the scenario that the SWP and SH groups have a similar history and a common origin, while a parallel evolving history is observed for the SWP/SH and MEP lineages.
In the MEP lineage, a population expansion event of M. davidii was occured around 79.17 ka BP (Table 3). This expansion thus took place during the end of the last interglacial period (130-75 ka BP), one of the warmer periods . The expansion time was younger than of the greater horseshoe bat (Rhinolophus ferrumequinum) for Japan (127-191 ka BP)  and was nearly the same as amphibians of the same region (80 ka BP) . The different expansion time with the greater horseshoe bat may be due to the influence of glacial oscillations at different latitudes at the same time . Bats are sensitive to environmental changes , such as the large-scale marine transgressions in the MEP region , which may influence bats' migratory behaviour and geographic distribution . Thus, a stepwise population expansion, inferred by individual assignment test of nDNA, has taken place in association with eco-environmental changes.
In the SWP region, a population expansion event of M. davidii occurred around 69.12 ka BP (Table 3), and originated from relatively high elevation areas in SWP to the MEP during 72-60 ka BP (early stage of the last glacial) . The expansion time was a little earlier than that of the greater horseshoe bat (R. ferrumequinum) for Europe (40-60 ka BP) and during the same period for Russia (14-81 ka BP) . Two haplotypes (21.43% individuals) were shared between YN populations in the SWP and AH1-JS populations in the MEP region. This conclusion is in agreement with previous studies on bats , birds  and amphibians . The result of nDNA analysis (Fig. 3) and an individual assignment test also indicated that some individuals of M. davidii shared the same nDNA alleles in the SWP and MEP regions, and no cline across these populations. The high elevation areas of the SWP region were the first to be affected by the sudden cold weather caused by the lowering of the snow line during the last glacial period (74-11.5 ka BP) . Therefore, some species or individuals in high altitudes may have immigrated to refuges, such as Yunnan, or to other regions in order to avoid such cold climatic conditions.
Based on the analysis of microsatellite data, a possible regional population contact was interpreted, which showed a similar allele site among populations in the SH group (Fig. 3). Although the phylogenetic tree was interrupted between SH and SWP, the relationship between SWP and SH groups was much closer than MEP based on the mtDNA results (Fig. 2), which may be an indication of SWP and SH origin from a common ancestor.
Refuges for Myotis davidii
At least two refuges have been reported in east China and in the lower elevations of the southwestern plateau [14, 35]. In our study, the existence of three mitochondrial lineages, with high nucleotide and allelic diversity directed towards three relict refuges for M. davidii in China, where an initial population expansion took place. There is no clear indication of where the refuge of M. davidii populations were located, but the possible areas are the low-elevation plateau, such as Yunnan area and Sichuan Basin (CQ2 is at its margin), where there were large-scale relict refuges for many species . In our study, the results of nDNA analysis imply a relatively simpler genetic structure for M. davidii in YN1-YN2 and CQ1, which is attributed to relict refuge during the last glaciations. The eastern coast and south hills all were glacial refuges which were indicated by the high nucleotide and allelic diversity of both markers. This trend also has been demonstrated for other species .
Genetically divergent populations are increasingly being recognized as appropriate units for conservation. The results presented here potentially indicate three conservation management units (MEP, SWP and SH). For the different management units, we suggest that protection require several different measures. As for the MEP region, in considering frequent contact among populations, there may be no barrier in population migration due to the lack of mountains. Therefore, it is important to strengthen the protection of these habitats. In the SWP and SH regions, the gene exchange among M. davidii populations may be impeded by the limit of landscape. Therefore, it is important to protect migration corridors. In addition, CQ2 is located in the Three Gorges Area at the entrance of the Sichuan basin. Small mammal fossils in the Three Gorges Area indicate that it provids a corridor for the exchange of fauna occurring in temperate and more tropical areas to the south . The position of CQ2 in the phylogenetic tree also points to a corridor of genetic exchange between the MEP and SWP/SH regions. In addition, the finding that Fst and Φst values between CQ2 and other populations within the MEP region exceed values among all the other populations suggests that populations of M. davidii are especially vulnerable to small population effects and fragmentation. A corridor of genetic exchange, such as CQ2, and prey or colony habitats should be protected to maintain population contact and its population persistence.
M. davidii is considered to be of least concern by IUCN . Also, before 2009 there were no systematic studies made of the distribution limit, status, and phylogeography of Chinese endemic bat species. Our research has shown that the distribution range of this species is more extensive than indicated by previous records. M. davidii also were found to be scarce and were not found in most recorded locations during the last nine years. Therefore, the public, environmental protection agencies, and various stakeholders should spark efforts for its future protection.
High haplotype and nucleotide diversity indicate that M. davidii has had a long evolutionary history. Based on the analysis of MDIV, the differentiation times among the three regions were during the Riss glaciation (210-135 ka BP), and the population differentiations correspond to a series of geological events and glacial cycles. The populations in the MEP and SWP experienced expansion events at about 79.17 ka BP and 69.12 ka BP, respectively. The phylogeographic study showed deep genetic differentiations among populations of M. davidii. Glacial cycles, ecological barriers and topological differences are associated with the extant distribution patterns and phylogenetic clades of M. davidii. The populations of M. davidii formed multihierarchical population structures and also conserved high genetic diversity in each region. The degree of genetic differentiation, based on the analysis of mtDNA and nDNA of M. davidii, suggests a male-mediated gene flow among populations or regions.
Sampling and Mitochondrial DNA amplification
We initiated our investigation of M. davidii in 2001. For 9 years, 126 individuals of M. davidii were collected from 17 localities (Fig. 1). All tissues were collected with a non-lethal method: this involved taking a wing punch biopsy for DNA analysis . Protocols for capturing and handling live bats followed the Regulations of Wildlife Conservation of the People's Republic of China (No. 24) , which were approved by the Wildlife Conservation Association of China.
Total genomic DNA was isolated from tissue using the Qiagen DNAeasy Tissue Kit. Hypervariable region I (HVI) of mitochondrial control region was amplified by polymerase chain reaction using primer pair P-F . In order to guarantee accurate sequencing, all amplification and sequencing was repeated at least once. All PCR reactions were performed in 25.0 μL volumes, containing approximately 10-30 ng of genomic DNA, 2.5 mM MgCl2, 75 mM Tris-HCl (pH9.0), 20 mM (NH4)2SO4, 0.01% (v/v) Tween-20, 0.2 mM dNTP mix, 10 μM of each primer and 1U of Bioline Taq polymerase. The following reaction conditions were used: 3 min at 94°C; 40 cycles of 94°C for 1 min, 55°C for 1 min, 72°C for 1 min; and 72°C for 10 min. DNA sequencing was performed on an ABI 3730 automated DNA sequencers (Applied Biosystems). DNA sequences were edited and aligned using BioEdit 220.127.116.11 . All haplotypes were submitted to GenBank [GenBank: GU013475-GU013527].
Subsequent to screening 16 microsatellite loci originally isolated from either Myotis myotis  or M. daubentonii , 8 polymorphic loci were selected for examining genetic variation in M. davidii (Table 2). All forward primers were labelled with 5'-Fluorescent bases (Table 2), and PCR reactions were conducted in 10 μL volumes containing 10 ng of template DNA, 1.5 mM MgCl2, 75 mM Tris-HCl (pH 9.0), 20 mM (NH4)2SO4, 0.01% (v/v) Tween-20, 0.2 mM dNTP mix, 1 μM of each primer and 0.5 U of Bioline Taq polymerase. Touchdown PCR was used with all primers, with an initial denaturing step of 3 min at 94°C. The annealing temperature of the reaction is decreased 0.5°C each cycle from 60°C to a touchdown at 50°C, at which temperature 20 cycles are carried out. An additional 25 cycles were then performed with 30 s denaturing at 91°C and 30 s annealing at 50°C. No extension steps were included , apart from a 5 min period at 72°C following the final annealing step. Genotyping was performed on an ABI 3730 automated DNA sequencer (Applied Biosystems, Inc.).
Estimates of mitochondrial haplotype variability within different regions and populations included the number of haplotypes (A), haplotype diversity (h) and nucleotide diversity (π). All estimates were made using ARLEQUIN 3.1 . Values for polymorphic sites and parsimony informative sites were estimated in DnaSP, version 4.0 .
Microsatellite genotype data were scored using the Genemarker software 1.75 (SoftGenetics Inc.). Genetic diversity was assessed for each population by calculating expected heterozygosity (HE), observed heterozygosity (HO) and average allelic richness (RS) determined with FSTAT 2.9.3 . Multi-locus Fis was calculated for each population and adjusted for multiple tests using a Bonferroni's correction . Deviation from the Hardy-Weinberg equilibrium (HWE) was tested by permutation using FSTAT 2.9.3 . Tests for linkage disequilibrium between loci for each population were performed with GENEPOP 3.4 .
Phylogenetic analyses of unique haplotypes included both maximum parsimony (MP) and maximum likelihood (ML) algorithms performed in PAUP 4.0b10 . A GTR + G + I model was used, as determined by Modeltest 3.7  (base frequencies: A, 0.3870; C, 0.2174; G, 0.1122; and T, 0.2834; transition/transversion ratio = 4.9960; proportion of invariable sites Pinvar = 0.0892; gamma distribution shape parameter = 2.4244). MP analyses were conducted using a heuristic search option with 100 random sequence additions and tree-bisection-reconnection (TBR) branch swapping. Robustness of the MP trees was assessed by 1000 bootstrap replicates. ML analyses used parameters estimated from trees obtained from MP analyses. ML analyses used heuristic searches with starting trees obtained by NJ followed by TBR branch-swapping. ML nonparametric bootstrap analyses used 100 heuristic searches with starting trees obtained with NJ based on p distances followed by TBR and nearest-neighbor interchange branch-swapping, saving all optimal trees. Bayesian inference was performed with MrBayes 3.1  using default parameters. Two independent parallel runs of four incrementally heated Metropolis-coupled MCMCs (Monte Carlo Markov Chains) were performed, with trees sampled every 100 generations for 1000000 generations to the average standard deviation of split frequency below 0.01. The first 10% of the trees were discarded as 'burnin', and posterior probabilities were estimated for the remaining saved trees. Trees were rooted with sequences of M. daubentoni, M. lucifugus, M. adversus, M. altarium, M. bombinus, M. chinensis and M. pequinius [GenBank: EU447268, U95342, U95341, U944765, GU944764, GU936479, GU944763].
The population demographic expansion was tested by mismatch distribution analysis and neutrality tests. Mismatch distribution analyses were implemented to detect historical demographic expansions. Significant difference from a model of sudden expansion was assessed using the sum of squared deviations (SSD) and the Harpending raggedness index (r) with parametric bootstrapping (10000 replicates). Generally, a smaller and non-significant value (PSSD and Pr > 0.05) indicated population expansion. We estimated the time since expansion (t) with the equation τ = 2 μkt , where τ (tau) is a parameter of the time to expansion in units of mutations, μ is the mutation rate per nucleotide and k is the number of nucleotides of the sequence. Tajima's D and Fu's Fs (based on the infinite-site model) were sensitive to the population expansion. Significant negative values of Tajima's D were interpreted as non-neutral conditions. Negative values of Tajima's D and Fu' Fs significantly deviated from "0", inferred as a signal of demographic expansion. All analyses were conducted with ARLEQUIN 3.1 .
With a coalescent-based approach, the program MDIV was used to estimate the timing of population divergence . MDIV was implemented using a 'finite model' (HKY), with 5000000 Markov chain iterations and a 25% burn-in. Mmax and Tmax were set to 10 and 5, respectively. The divergent times (t) were estimated using the formula t = T pop *(θ/2 μk), where T pop is the time of population divergence, θ is the population mutation rate, μ is mutation rate per nucleotide, and k is the number of nucleotides of the sequence. Likelihood values for T pop and θ were calculated and the value with the highest posterior probability accepted as the best estimate . The program was run multiple times with different random seeds to obtain consistent distribution results.
We used a mutation rate of 20% per million years (Myr) in our study, which was previously applied to bats of the genus Nyctalus . The generation time was estimated to be 2 years.
Population structure analysis
To test for geographical genetic structure, analyses of molecular variance (AMOVA) with 10000 permutations were assessed according to the degree of differentiation between regions (ΦCT), between populations within regions (ΦSC) and between all populations (ΦST) by using ARLEQUIN 3.1 .
To compare patterns of structure obtained from mtDNA haplotype data with those from multi-locus markers, we used STRUCTURE 2.2  to reconstruct the hierarchical relationship among populations, as well as to distinguish historical processes. We applied a mixed model that allows for allele frequency correlation across a set of K genetic clusters. We performed 10 replicate  runs of structure for each value of K. We applied the admixture model with a burn-in of 30000 and a run length of 106. Summary outputs for each value of K were then displayed graphically using DISTRUCT . Assignment tests can be used to identify the source population of individuals, and individuals of unknown origin can be pre-assigned to populations according log likelihood scores . Therefore, we used assignment tests to identify the origin population according to likelihoods from genotypic data , which was implemented with GENALEX 6.0 .
A Mantel test (10000 permutations) to assess the statistical significance of the correlation between genetic and geographical distances was performed with IBDWS 3.14 , to understand the influence of geographical barriers on population differentiation.
We thank Jiangfeng Du for participate in its design and coordination and helped to draft the manuscript. We thank Professor Kaiya Zhou, Professor Bao Liu and Master Xinmin Tian for provide technical help. We thank Professor Fumin Lei and Doctor Gang Song for chorography providing. We thank Shi Li, Limin Shi, Xinhua Wang and Xu Zhu for field support. We thank Zhenzhen Zhang for molecular data providing. This work was supported by the National Natural Science Foundation of China (30770361, 30900132); the National Grand Fundamental Research 973 Program of China (2009CB426305); the Program of Introducing Talents of Discipline to Universities (B07017).
- Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405 (6789): 907-913. 10.1038/35016000.View ArticlePubMedGoogle Scholar
- Avise JC: Phylogeography: the history and formation of species. 2000, Cambridge: Harvard University PressGoogle Scholar
- Lovette IJ: Glacial cycles and the tempo of avian speciation. Trends in Ecology & Evolution. 2005, 20 (2): 57-59. 10.1016/j.tree.2004.11.011.View ArticleGoogle Scholar
- Zhang D, Fengquan L, Jianmin B: Eco-environmental effects of the Qinghai-Tibet Plateau uplift during the Quaternary in China. Environmental Geology. 2000, 39 (12): 1352-1358. 10.1007/s002540000174.View ArticleGoogle Scholar
- Tong G, Chen Y, Wu X, Li Z, Yang Z, Wang S, Cao J: Pleistocene environmental megaevolution as indicated by the sporopollen floras in China. Journal of Geomechanics. 1999, 5 (4): 11-21.Google Scholar
- Smith AT, Johnston CH, Jones G, Rossiter S: Myotis davidii. IUCN. 2009, [http://www.iucnredlist.org]Google Scholar
- Allen GM: The Mammals of China and Mongolia. 1938, New York: American Museum of Natural History PressView ArticleGoogle Scholar
- Norberg UM, Rayner JMV: Ecological Morphology and Flight in Bats (Mammalia; Chiroptera): Wing Adaptations, Flight Performance, Foraging Strategy and Echolocation. Philosophical Transactions of the Royal Society of London Series B, Biological Sciences. 1987, 316 (1179): 335-427. 10.1098/rstb.1987.0030.View ArticleGoogle Scholar
- Liu MY, Xie YH, Ji DM: The list of China vertebrate. 2000, China: Liaoning University PressGoogle Scholar
- Wu XZ, Wang YF, Zheng LP: Fauna in Yanger cave in Hunan Province and paleo-environment in southern margin of Three Gorges of the Yangtze River during the late Pleistocene. Journal of Chongqing Normal University. 2008, 3: 76-83.Google Scholar
- Jones G, Jacobs DS, Kunz TH, Willig MR, Racey PA: Carpe noctem: the importance of bats as bioindicators. Endangered Species Research. 2009, 8: 93-115. 10.3354/esr00182.View ArticleGoogle Scholar
- Xu L, He CF, Jiang TL, Shi LM, Sun KP, Berquist SW, Feng J: Phylogeography and population genetic structure of the great leaf-nose bat (Hipposideros armiger) in China. Journal of Heredity. 2010, doi:10.1093/jhered/esq039Google Scholar
- Song G, Qu Y, Yin Z, Li S, Liu N, Lei F: Phylogeography of the Alcippe morrisonia (Aves: Timaliidae): long population history beyond late Pleistocene glaciations. BMC Evolutionary Biology. 2009, 9 (1): 143-10.1186/1471-2148-9-143.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang H, Yan J, Zhang G, Zhou K: Phylogeography and demographic history of Chinese black-spotted frog populations (Pelophylax nigromaculata): Evidence for Independent refugia expansion and secondary contact. BMC Evolutionary Biology. 2008, 8 (1): 21-10.1186/1471-2148-8-21.PubMed CentralView ArticlePubMedGoogle Scholar
- Flanders J, Jones G, Benda P, Dieta C, Zhang SY, LI G, Sharifi M, Rossiter SJ: Phylogeography of the greater horseshoe bat, Rhinolophus ferrumequinum: contrasting results from mitochondrial and microsatellite data. Molecular Ecology. 2009, 18: 306-318. 10.1111/j.1365-294X.2008.04021.x.View ArticlePubMedGoogle Scholar
- Newman RA, Squire T: Microsatellite variation and fine-scale population structure in the wood frog (Rana sylvatica). Molecular Ecology. 2001, 10 (5): 1087-1100. 10.1046/j.1365-294X.2001.01255.x.View ArticlePubMedGoogle Scholar
- Hewitt GM, Butlin RK: Causes and consequences of population structure. Behavioural Ecology. 1997, 350-372.Google Scholar
- Castella V, Ruedi M, Excoffier L: Contrasted patterns of mitochondrial and nuclear structure among nursery colonies of the bat Myotis myotis. Journal of Evolutionary Biology. 2001, 14 (5): 708-720. 10.1046/j.1420-9101.2001.00331.x.View ArticleGoogle Scholar
- Ruedi M, Walter S, Fischer MC, Scaravelli D, Excoffier L, Heckel G: Italy as a major Ice Age refuge area for the bat Myotis myotis (Chiroptera: Vespertilionidae) in Europe. Molecular Ecology. 2008, 17 (7): 1801-1814. 10.1111/j.1365-294X.2008.03702.x.View ArticlePubMedGoogle Scholar
- Grant WAS, Bowen BW: Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. Journal of Heredity. 1998, 89 (5): 415-426. 10.1093/jhered/89.5.415.View ArticleGoogle Scholar
- Doebeli M, Dieckmann U: Speciation along environmental gradients. Nature. 2003, 421: 259-264. 10.1038/nature01274.View ArticlePubMedGoogle Scholar
- Storz JF, Bhat HR, Kunz TH: Genetic consequences of polygyny and social structure in an Indian fruit bat, Cynopterus sphinx. I. Inbreeding, outbreeding, and population subdivision. Evolution. 2001, 55 (6): 1215-1223.View ArticlePubMedGoogle Scholar
- Kerth G, Mayer F, Petit E: Extreme sex-biased dispersal in the communally breeding, nonmigratory Bechstein's bat (Myotis bechsteinii). Molecular Ecology. 2002, 11 (8): 1491-1498. 10.1046/j.1365-294X.2002.01528.x.View ArticlePubMedGoogle Scholar
- Yi C, Cui Z, Xiong H: Numerical periods of quaternary glaciations in China. Quaternary Sciences. 2005, 25 (5): 609-619.Google Scholar
- Ge X, Ren S, Ma L, Wu G, Liu Y, Yuan S: Multi-stage uplifts of the Qinghai-Tibet plateau and their environmental effects. Earth Science Frontiers. 2006, 13 (6): 118-130.Google Scholar
- Manel S, Schwartz MK, Luikart G, Taberlet P: Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution. 2003, 18 (4): 189-197. 10.1016/S0169-5347(03)00008-9.View ArticleGoogle Scholar
- Bond G, Showers W, Cheseby M, Lotti R, Almasi P, Demenocal P, Priore P, Cullen H, Hajdas I, Bonani G: A pervasive millennial-scale cycle in North Atlantic Holocene and glacial climates. Science. 1997, 278 (5341): 1257-1266. 10.1126/science.278.5341.1257.View ArticleGoogle Scholar
- Gillespie A, Molnar P: Asynchronous maximum advances of mountain and continental glaciers. Reviews of Geophysics. 1995, 33 (3): 311-364. 10.1029/95RG00995.View ArticleGoogle Scholar
- Liu D, Shi Y, Wang R, Zhao Q, Jian Z, Cheng X, Wang P, Wang S, Yuan B, Wu X, et al: Table of the Chinese quaternary stratigraphic correlation remarked with climate change. Quaternary Sciences. 2000, 20 (2): 108-128.Google Scholar
- Watson RT, Zinyowera MC, Moss RH: The regional impacts of climate change: an assessment of vulnerability. 1998, UK: Cambridge University PressGoogle Scholar
- Voris HK: Special Paper 2: maps of Pleistocene sea levels in Southeast Asia: Shorelines, river systems and time durations. Journal of Biogeography. 2000, 27 (5): 1153-1167. 10.1046/j.1365-2699.2000.00489.x.View ArticleGoogle Scholar
- Russell AL, Medellin RA, McCracken GF: Genetic variation and migration in the Mexican free-tailed bat (Tadarida brasiliensis mexicana). Molecular Ecology. 2005, 14 (7): 2207-2222. 10.1111/j.1365-294X.2005.02552.x.View ArticlePubMedGoogle Scholar
- Shi Y: A suggestion to improve the chronology of quaternary glaciations in China. Journal of Glaciology and Geocryology. 2002, 24 (6): 687-692.Google Scholar
- Liu JQ, Ni YY, Chu GQ: Main palaeoclimatic events in the Quaternary. Quaternary Science. 2001, 21 (3): 239-248.Google Scholar
- Zhang R, Jin S, Quan G, Li S, Ye Z, Wang F, Zhang M: Distribution of Mammalian Species in China. 1997, China: China Forestry Publishing House PressGoogle Scholar
- World Heritage Centre: Sichuan Giant Panda Sanctuaries-Wolong, Mt Siguniang and Jiajin Mountains. Natural Heritage List. 2006, [http://www.worldheritagesite.org/sites/giantpanda.html]Google Scholar
- Worthington WJ, Barratt E: A non-lethal method of tissue sampling for genetic studies of chiropterans. Bat Research News. 1996, 37 (1): 1-4.Google Scholar
- National People's Congress Standing Committee: Regulation of wildlife conservation of the People's Republic of China (No. 24). 1988Google Scholar
- Wilkinson GS, Chapman AM: Length and sequence variation in evening bat D-loop mtDNA. Genetics. 1991, 128 (3): 607-617.PubMed CentralPubMedGoogle Scholar
- Hall TA: BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium Series. 1999, 95-98.Google Scholar
- Castella V, Ruedi M: Characterization of highly variable microsatellite loci in the bat Myotis myotis (Chiroptera: Vespertilionidae). Molecular Ecology. 2000, 9 (7): 1000-1002. 10.1046/j.1365-294x.2000.00939-6.x.View ArticlePubMedGoogle Scholar
- Ngamprasertwong T, Mackie IJ, Racey PA, Piertney SB: Spatial distribution of mitochondrial and microsatellite DNA variation in Daubenton's bat within Scotland. Molecular Ecology. 2008, 17 (14): 3243-3258. 10.1111/j.1365-294X.2008.03845.x.View ArticlePubMedGoogle Scholar
- Don RH, Cox PT, Wainwright BJ, Baker K, Mattick JS: 'Touchdown' PCR to circumvent spurious priming during gene amplification. Nucleic Acids Research. 1991, 19 (14): 4008-10.1093/nar/19.14.4008.PubMed CentralView ArticlePubMedGoogle Scholar
- Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary bioinformatics. 2005, 1: 47-50.Google Scholar
- Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003, 19 (18): 2496-2497. 10.1093/bioinformatics/btg359.View ArticlePubMedGoogle Scholar
- Goudet J: FSTAT, a program to estimate and test gene diversities and fixation indexes (updated from Goudet 1995). 2001, [http://www2.unil.ch/popgen/softwares/fstat.htm]2.3.9Google Scholar
- Raymond M, Rousset F: GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism. Journal of Heredity. 1995, 86 (3): 248-249.Google Scholar
- Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony. 2001, Sinauer, Sunderland, MA, 4.0 b10Google Scholar
- Posada D, Crandall KA: Modeltest: testing the model of DNA substitution. Bioinformatics. 1998, 817-818. 10.1093/bioinformatics/14.9.817.Google Scholar
- Huelsenbeck JP, Ronquist F, Nielsen R, Bollback JP: Bayesian inference of phylogeny and its impact on evolutionary biology. Science. 2001, 294 (5550): 2310-10.1126/science.1065889.View ArticlePubMedGoogle Scholar
- Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology & Evolution. 1992, 9 (3): 552-569.Google Scholar
- Nielsen R, Wakeley J: Distinguishing migration from isolation a markov chain Monte Carlo approach. Genetics. 2001, 158 (2): 885-896.PubMed CentralPubMedGoogle Scholar
- Petit E, Excoffier L, Mayer F: No evidence of bottleneck in the postglacial recolonization of Europe by the noctule bat (Nyctalus noctula). Evolution. 1999, 53 (4): 1247-1258. 10.2307/2640827.View ArticleGoogle Scholar
- Peakall ROD, Smouse PE: GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes. 2006, 6 (1): 288-295. 10.1111/j.1471-8286.2005.01155.x.View ArticleGoogle Scholar
- Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data. Genetics. 2000, 155 (2): 945-959.PubMed CentralPubMedGoogle Scholar
- Rosenberg NA: DISTRUCT: a program for the graphical display of population structure. Molecular Ecology Notes. 2004, 4 (1): 137-138. 10.1046/j.1471-8286.2003.00566.x.View ArticleGoogle Scholar
- Paetkau D, Calvert W, Stirling I, Strobeck C: Microsatellite analysis of population structure in Canadian polar bears. Molecular Ecology. 1995, 4 (3): 347-354. 10.1111/j.1365-294X.1995.tb00227.x.View ArticlePubMedGoogle Scholar
- Jensen JL, Bohonak AJ, Kelley ST: Isolation by distance, web service. BMC genetics. 2005, 6 (1): 13-10.1186/1471-2156-6-13.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.