- Research article
- Open Access
Regionally and climatically restricted patterns of distribution of genetic diversity in a migratory bat species, Miniopterus schreibersii(Chiroptera: Vespertilionidae)
BMC Evolutionary Biology volume 8, Article number: 209 (2008)
Various mechanisms such as geographic barriers and glacial episodes have been proposed as determinants of intra-specific and inter-specific differentiation of populations, and the distribution of their genetic diversity. More recently, habitat and climate differences, and corresponding adaptations have been shown to be forces influencing the phylogeographic evolution of some vertebrates. In this study, we examined the contribution of these various factors on the genetic differentiation of the bent-winged bat, Miniopterus schreibersii, in southeastern Europe and Anatolia.
Results and conclusion
Our results showed differentiation in mitochondrial DNA coupled with weaker nuclear differentiation. We found evidence for restriction of lineages to geographical areas for hundreds of generations. The results showed that the most likely ancestral haplotype was restricted to the same geographic area (the Balkans) for at least 6,000 years. We were able to delineate the migration routes during the population expansion process, which followed the coasts and the inland for different nested mitochondrial clades. Hence, we were able to describe a scenario showing how multiple biotic and abiotic events including glacial periods, climate and historical dispersal patterns complemented each other in causing regional and local differentiation within a species.
The geographic transitions between continents can result in species diversification and endemism, forming regions of allopatry, and contact zones for divergent flora and fauna. Sulawesi, for instance, an island located between Continental Asia and Australia, has elements of both faunal assemblages within its boundaries, hosting an elevated number of species and endemics . Mexico, as another example, being located in a transition zone between tropical central America, and temperate North America is considered to be a megadiversity country . Within the temperate zones, southeastern Europe and Anatolia are located within a similar geographical transition centered in between Europe, Asia and Africa (Figure 1a). As such, this region comprises an interesting area for investigation as a zone of allopatry with geographical barriers to gene flow [3, 4], a contact zone for divergent biota from different continents and climatic regimes, and a refugium for the entire western Palearctics [5, 6].
In this study, we examined the nuclear and mitochondrial genetic structure of the bent-winged bat, Miniopterus schreibersii, in southeastern Europe and Anatolia. M. schreibersii is a colonial, cave-dwelling , and polytypic species with one of the largest Old World ranges among mammals. Its global distribution spans the Palearctics, Africa and Australia, with 17 subspecies described [8, 9]. Recent studies showed high levels of nuclear and mitochondrial genetic structure in a congeneric species, Minipterus natalensis in South Africa . Differentiation has also been reflected in taxonomy of this species in Australia , and southeastern Europe and Anatolia , where three and two subspecies have been described, respectively.
Observing the genetic differentiation in this species in southeastern Europe and Anatolia, our goal was to try teasing apart the differential contribution of various factors that can contribute to its genetic differentiation. Topographical features in the region, including the Taurus Mountains, the eastern Anatolian Diagonal Mountain Chain and the Sea of Marmara (Figure 1a) are potential barriers to gene flow that could result in allopatric differentiation of populations. For instance, The Marmara was proposed as a potential barrier preventing the secondary contact of isolated refugial populations in another bat species, the greater horseshoe bat . The region is also interesting due to different climates regimes that co-exist, as differential adaptations to these could promote population divergence. In addition to drastic differences in precipitation between the coastal and inland regions, there also exist regional differences in climate. For instance, the Mediterranean coast, the Black Sea coast and the Marmara region all show slight differences in terms of temperature and precipitation patterns . The glacial history of the region, especially at the Pleistocene, is another candidate for having caused divergence of populations. The Balkans is one of the main glacial refugial areas for the entire Europe [15, 16], and the Caucasus is an important refugium for western Asia and the Middle East in general [17, 18]. Anatolia is located in between these two major refugial areas, and this proximity could have left a signature on the genetic make-up of populations. To address these questions, we used nuclear microsatellite markers, as well as mtDNA markers, comprising the next logical step to the previous study of Bilgin et al. , which had used only mtDNA data to mainly evaluate the taxonomic status of M. schreibersii in this region. In the end, we were able to construct a scenario outlining the historical details of the genetic differentiation of M. schreibersii in the transition between Europe and Anatolia.
121 individuals of M. schreibersii were sampled in 34 locations in Bulgaria, Greece and Turkey (Figure 1a and Table 1), spanning a range of about 1700 km. This comprises an increase from 58 individuals in 14 locations that was used in the previous study of Bilgin et al. , on the same species. The caves where samples from Bilgin et al.  were used in the current study are indicated in Table 1. Fieldwork was undertaken, in the years 2004 and 2005, in August, to avoid disturbing the nursery colonies. Through inspection of the degree of fusion of the phalangeal epiphyses, the adult status of the sampled bats was confirmed, and no juvenile bats were sampled . Tissue for genetics analysis was collected from each of the wing membranes (plagiopatagium) of individual bats by using biopsy-punches (3 mm diameter) as outlined by Worthington Wilmer and Barratt . The 3 mm holes in each wing are known to heal in approximately four weeks.
For mtDNA analysis, the hyper variable region I of the control region (HV1), tRNA-proline and tRNA-threonine genes, and a partial cytochrome b sequence were sequenced. Half of a biopsy punch was used for each individual's DNA extraction, with a DNeasy Extraction Kit, following the manufacturer's protocols (QIAGEN, Valencia, CA). The primers C and E, and the PCR conditions were used as described in Wilkinson and Chapman  for DNA amplification. This was followed by cycle sequencing, both in 5' and 3' directions, using the primers C and E, respectively. This involved 25 cycles in 10 μl reactions, which were composed of 1 μl of PCR product, 5.7 μl of H2O, 0.3 μl of primer (20 μM), 1 μl of florescent dye (ABI Big Dye) and 2 μl of 5× buffer (provided with the florescent dye). The cycle sequencing parameters for each cycle were 10 seconds of denaturation at 96°C, 5 seconds of annealing at 50°C and 4 minutes of extension at 68°C, followed by a final extension of 7 minutes at 72°C. The sequencing reactions were cleaned using ethanol precipitation and run in an ABI 3730 automated sequencer following the manufacturer's protocol (Applied Biosystems, Inc.). The resulting sequences were assembled using Sequencher 4.1 (Gene Codes Corp.), and aligned using Clustal X  prior to further data analysis.
Four microsatellite primer pairs (Mschreib2, Mschreib3, Mschreib4 and Mschreib5) designed for M. schreibersii  and two pairs (D5S1457 and D6S271) designed for humans, which also amplify in baboons , were used to amplify and genotype six nuclear loci. The PCR amplifications were carried out with a QIAGEN Multiplex kit following the manufacturer's protocol. The reaction volume was cut down by a factor of five (10 μl instead of the manufacturer's recommended 50 μl) to decrease the overall cost, while 1 μl of template DNA was used. The optimal annealing temperature for each locus was found through a gradient PCR. The PCR reactions were successful at the annealing temperature of 57°C for all loci, except Mschreib2 and Mschreib4, which yielded amplicons at 61°C. A homozygote individual was sequenced for each locus using the respective forward and reverse primers, so as to confirm the identity of the repeat motifs with the published ones. The PCR products were run in an ABI Automated Sequencer 3730, using a Rox label (ABI, CA) to size the fragments, in a volume composed of 8.9 μl of High Dye (Applied Biosystems, Inc.), 0.1 μl of Rox, and 1 μl of PCR product for each sample. The final sizes of the fragments were computed in GeneMapper (Applied Biosystems, Inc.) and scored manually. The homozygous individuals were amplified and scored three times to avoid allelic dropout and confirm allele sizes; as for the heterozygous individuals, they were amplified and scored until the sizes matched in two separate runs.
Basic descriptive statistics and genetic diversity parameters such as haplotype diversity (h), nucleotide diversity (π) and number of polymorphic sites [24, 25], were calculated in DnaSP 4.20.2 . These calculations were also made for the reciprocally monophyletic clades S and P, see below, separately. We reconstructed the phylogeny of the haplotypes using distance (neighbor-joining, NJ) and maximum parsimony analyses, and the software PAUP* 4.06b . Nodal support for NJ and maximum parsimony trees were evaluated by 1000 bootstrap replicates. The substitution model to be used for building the distance matrix was selected by applying a likelihood ratio test for the goodness of fit of various substitution models to the data by using the program Modeltest 3.7 . For the neighbor joining tree, the model of nucleotide evolution, selected by Modeltest, which fit the data best was that of Hasegawa et al.  [(HKY85) + Γ , shape parameter = 0.07] and was used in the construction of the NJ tree. In the maximum parsimony analysis, TBR (tree bisection and reconnection) was used as the branch-swapping algorithm as it is the most effective routine for recovering an optimum set of cladograms . Gaps were treated as missing data. The sequences were checked for any evidence for saturation, using percent genetic divergence versus transition/transversion plots in DAMBE . When there was evidence for saturation, the parsimony analysis was adjusted using a step matrix by weighting transversions inversely against transitions, based on their frequency. The ratio of transitions to transversions was calculated following Nei and Kumar  in MEGA4 (Tamura et al., 2007). MEGA4 was also used to compute the genetic distances between clades S and P, using maximum likelihood composite distances. Average number of substitutions per site (Dxy) was calculated using DnaSP 4.20.2.
Differentiation in mtDNA was also explored using analysis of molecular variance (AMOVA) , by calculating the Fst analogue Φst. The program ARLEQUIN 2.0  was used to make the AMOVA calculations for mtDNA. The AMOVA calculations included two separate groupings. One of these was based on grouping clades S and P as separate populations. The second set of regional groups was prepared to test for the effect of geographical barriers to gene flow. These included four groups, which comprised sampling sites to the west of the Marmara Sea (sites 1–18 in Figure 1), between the Marmara Sea, the Taurus and the eastern Anatolian Diagonal (sites 13–24, 30), to the south of the Taurus (25–29), and to the east of the eastern Anatolian Diagonal (sites 31–34).
Nested clade analysis (NCA) was used to determine the causes of intraspecific structure, if any. The program TCS  was used to build a statistical parsimony network and subsequently GeoDIS 2.0  was used for making an NCA. We used the inference key of 11/11/2005 from http://darwin.uvigo.es/download/geodisKey_11Nov05.pdf
Evidence for population expansion was explored using the neutrality test statistics R2  and F S  in DnaSP. Based on simulation models, among various statistics that determine signatures of past population growth, R2 and F S are the two that have the greatest statistical power . The significance of these statistics was assessed with 1000 coalescent simulations. In addition, Bayesian Skyline plots  were used to detect changes in effective population sizes in the past. We estimated the divergence times of reciprocally monophyletic clades using a Bayesian coalescent approach implemented in BEAST 1.4.6. . The HKY + Γ, as revealed by Modeltest, were used as the substitution and site heterogeneity models, respectively. Again, based on Modeltest, The Γ prior was set to 0.07. The tree model prior was determined by using a UPGMA tree to construct an initial tree, and the mean rate and covariance priors were set to uniform values. The chain was run for 100,000,000 generations, making sure that the ESS values of each statistic was at least 1000. Convergence was examined in TRACER 1.4 http://beast.bio.ed.ac.uk/Tracer.
Basic tests of Hardy-Weinberg equilibrium and linkage disequilibrium were carried out in FSTAT . Sequential Bonferroni corrections were made to correct for levels of significance for multiple tests . The software Micro-Checker  was used to test for the presence of stutter bands, large allelic drop-out, and null-alleles. The unbiased number of alleles (based on sample size) was calculated in FSTAT, using a method of rarefaction . The expected and observed heterozygosities and their P values (for the entire data set, and separately for clades S and P), the uncorrected number of alleles, allelic diversity (in terms of effective number of alleles, Ne), and frequency of alleles in clades S and P for each locus, are provided in the additional file. All of these computations were made using Genalex 6.1 . In addition, corrected number of private alleles were calculated using a rarefaction method, and the software HP-rare . Global and locus-by-locus AMOVA were made, and unbiased Fst and Rst estimators [of θst  and ρst , respectively] were also calculated using Genalex 6.1. Both Fst and Rst estimators were calculated to account for the possible underestimation of Fst in detecting genetic differentiation when analyzing microsatellites .
For the analysis of mtDNA, 433 base pairs of mtDNA spanning cytochrome b (58 bp), tRNA-threonine (71 bp), tRNA-proline (66 bp) and HV1 (238 bp) were sequenced for 121 individuals. These sequences have been deposited to GenBank with the accession numbers EU332355–EU332392 and EU332393–EU332402. In this fragment, there were 52 polymorphic sites, 43 of which were parsimony informative. A total of 49 haplotypes were found. The haplotype diversity (h) was 0.926 (S.D. = 0.015) and nucleotide diversity (π) was 0.02593 (S.D. = 0.00327).
Phylogenetic trees for M. schreibersii were constructed using neighbor joining and parsimony methods (Figure 2a). A sample of M. schreibersii from Indonesia was included in the analysis as an outgroup. In the neighbor joining tree, two monophyletic clades have been identified, clade S and clade P, with high bootstrap support (83% and 100%, respectively). A heuristic maximum parsimony search also gave high bootstrap support to the nodes delimiting clades S (98%) and P (81%). In this analysis, the transversions were weighted against transitions based on their frequency (1:7). The C.I and R.I were both 0.908 for this tree.
In terms of descriptive statistics, a total of 38 haplotypes were found as belonging to clade S, in 102 individuals. Within clade S, there were a total of 37 polymorphic sites, 24 of which were parsimony-informative. Haplotype diversity was 0.932 (S.D. = 0.014) and nucleotide diversity was 0.00689 (S.D. = 0.00043). Within clade P, there were ten haplotypes, found in 19 bats, comprised of ten polymorphic sites of which six were parsimony-informative. Haplotype diversity was 0.897 (S.D. = 0.056) and nucleotide diversity was 0.0056 (S.D. = 0.00089).
The number of fixed and percent nucleotide differences between the haplotype groups is provided in Table 2, to give an idea of the extent of differentiation of groups. The percent divergence (modified from Dxy) between clades S and P was 7.94%. The maximum composite likelihood distance between the clades was almost identical (0.078, S.E. = 0.015). The differences between the Indonesian sample and clades S and P were 12.19% and 8.29%, respectively. These differences between regions are indicative of differentiation when compared to the percent nucleotide differences within each haplogroup, which ranged between 0 and 0.68%.
The Bayesian estimate of the time of the split of the S and P clades was calculated as 233,000 years BP (95% Highest Posterior Density = 169,000–299,000 years BP). In these calculations, a molecular clock rate based on 20% per million years differentiation of D-loop in the noctule bat , which falls between the estimate for house mouse [10% per million years ) and bison (30% per million years ] was used.
Another mathematical treatment of the differentiation between clades S and P was made with an AMOVA. The pairwise Φst was 0.92 and significantly greater than zero (P < 0.01, 1000 permutations). The percent of variance attributable to differences between clades was 92%, whereas this variance value within clades was 8% (Table 3). In terms of geographical distribution, both haplotypes belonging to clade S and P appeared on either side of any putative geographic barrier. In addition, using the tree approach, within either clade P or S, there didn't seem to be any geographically meaningful clusters. The AMOVA results, when made using the potential geographic barriers to define populations, showed that the percent of variation attributable to among group differences was 27% (Table 4).
Several additional analyses were made to detect the causes of this genetic structure within clades S and P by building a statistical parsimony network for each (Figure 2b and Figure 2c, respectively). Both of these networks showed star-like phylogenies, which is characteristic of expanding populations . For clade S, haplotypes comprising nodes of expansion are colored in black (Figure 2b). A similar pattern of a star-like network was found for clade P, however there was only one potentially ancestral haplotype (P1, Figure 2c).
Given this pattern of star-like networks, indicative of population expansion, the significance of these patterns was checked through the statistics R2 and F S . As the AMOVA indicated a significant genetic break between clade S and clade P, the analyses were run separately for each clade. For clade S, both of these statistics, whose values were 0.0385 and -32.2963, respectively, indicated significant expansions, as the probability of getting values lower than these was < 0.01 (1000 replicates). For clade P, the pattern was a similar one of significant R2 and F S , whose values were 0.0880 (P < 0.05) and -4.8960 (P < 0.01), respectively. Also the frequency of pairwise nucleotide differences of haplotypes was plotted with expectations under a constant-size model and a model of range expansion. The plots for both clade S and P fit a range expansion model better than a constant size model (Figures 3a–d). The Bayesian skyline plots for clades S and P also showed patterns of population expansion, initiating around 15,000 (Figure 4a) and 5,500 (Figure 4b) years ago, respectively.
Nested Clade Analysis (NCA) was applied to the statistical parsimony networks to infer if there were any geographic patterns within clades S and P. During the building of the nested cladogram, when a particular haplotype could not immediately be assigned to a nesting clade due to its position, it was included in the clade that had the smallest sample size, to increase the statistical strength of that nesting clade. In clade P, no significant associations between haplotypes were found and the hypothesis of panmixia could not be rejected. The nested cladograms for clade S showed four different second level nested clades (Figure 1b), with some significant associations. Based on expectations from the coalescent theory, the haplotype S15, within the purple nested clade (2.1), being the haplotype with the highest frequency, and the one that gave rise to the greatest number of haplotypes, was the most likely ancestral haplotype. This haplotype showed a distribution almost exclusively restricted to the Balkans (Figure 1a). The central haplotype (S3) of the light blue nested clade (2.3) budded off from the ancestral nested clade 2.1, and had a distribution along the Mediterranean coast. The red nested clade (2.4), which was directly connected to S3, was exclusively found along the Mediterranean coast. The green nested clade (2.2), which also derived out of the purple nested clade, showed a distribution predominantly restricted to the Black Sea coast.
In terms of climate, the average for annual precipitation in localities for clade P bats was lower than that for clade S bats (Table 5). The difference was statistically significant (P = 7.84*10-33, unpaired t-test, N = 33). With respect to the precipitation values, there was no overlap in the confidence intervals of the localities of clades S and P.
In microsatellite analyses, none of the loci showed any evidence for linkage disequilibrium, after sequential Bonferroni corrections. The loci Mschreib4 and Mschreib5 were found not to be in Hardy-Weinberg equilibrium. A separate analysis of Hardy-Weinberg equilibrium, for each mtDNA clade, indicated that clade S was not in equilibrium and this probably was the cause for the overall disequilibrium when the two clades were analyzed together. The disequilibrium in these two loci was due to the presence of null alleles, and the microsatellite allele and genotype frequencies were adjusted, using the program Micro-Checker, prior to further data analyses.
A pattern of allelic differences was seen in loci Mschreib2 and Mschreib3, with Clade P individuals having relatively smaller sized alleles. There were private alleles in each population, and they had relatively low frequencies (less than 5%), except three alleles, which had frequencies over 10%. These were alleles 194 in locus Mschreib2 (31.3%), 133 in locus Mschreib3 (16.7%) in clade P, and 190 in locus Mschreib5 (13.4%) in clade S. There were also some differences in the frequencies of certain alleles, although they were not private alleles. For instance the allele size range series between 130–140 bp had a high frequency (33.4%) in clade P, however only two alleles in this series (133 and 137 bp) were found in clade S, in one individual each. For these two loci, based on the corrected number of private alleles, these differences were reflected as a higher number of private alleles for clade P, when compared to clade S (Additional File, Table AF7). Although to a smaller extent, the AMOVA of microsatellite loci, compared based on the two mtDNA clades also supported the break (Table 6), with a cumulative Fst value 0.023 that was significantly greater than zero (P = 0.010, 1000 permutations). The Rst value was comparable, but slightly higher (0.037), suggesting that Fst slightly underestimated genetic differentiation when compared to Rst.
Anatolian Suture Zone
Suture zones are often observed after postglacial range expansions where lineages that diverged in allopatry meet again . In Europe, there are many different examples of an east/west division of lineages of animals and plants that diverged in separate refugia . Genetic suture zones were also observed for bats in the Balkans and Anatolia [53–55]. In this study, for M. schreibersii, a similar suture zone, of two divergent mitochondrial clades, was seen within Anatolia. Although to a smaller extent, this differentiation in mtDNA was also reflected in the nuclear genome. The Bayesian estimate of the time to the most recent common ancestor of these two clades was 233,000 years BP (169,000–299,000 years BP), supporting a scenario of secondary contact after differentiation in isolation during the Pleistocene. The geographical distribution of clade S and clade P are in conformation with the previous study of Bilgin et al. , with clade S being found predominantly to the west, and clade P to the east. These results support the idea that the clade S differentiated in the Balkans, and clade P in refugia to the east of the Turkish border, in Armenia, Georgia or Iran during the glacial maxima in the last ice age. After the conditions became more hospitable, the two clades expanded into their current distribution, forming a suture zone in between, due to secondary contact.
Evaluation of whether the observed mtDNA differentiation was reflected in the nuclear genome, using microsatellites, did not provide any strong evidence for the presence of reproductive isolation and cryptic speciation. Although there was some statistically significant differentiation in microsatellites, the actual quantitative levels were low. Considering the nuclear locus that reflects the mitochondrial divide the best, Mschreib2, allele 194 was found only in individuals belonging to clade P with a frequency of 31.3%. Putting into consideration the point that none of these were homozygotes, 62.5% of the individuals had this allele. However, even with this frequency, it is not possible to assign the remaining 37.5% of individuals to either clade based on this locus. Another locus in M. schreibersii, Mschreib3, showed a similar break with all of the eighteen typed individuals except two in clade P having at least one allele smaller than 140 bases. However, still, two individuals in clade S had alleles in this size range, suggesting that a hypothesis of nuclear gene-flow among the two clades cannot be rejected. Hence it might be possible to conclude that although differentiation in refugia might have resulted in frequency differences in certain allele/locus combinations, this did not necessarily result in reproductive isolation.
Considering the details of the geographic distribution of each clade, clade S showed more of a coastal distribution, as opposed to clade P, which had a distribution restricted to inland. This kind of coastal vs. inland genetic differentiation has been observed in other species including mammals, e.g., Atlantic tree rats , lesser long-nosed bat , birds (e.g., swamp sparrow ), and plants (e.g., western white pine ). The significant differences between climate parameters with respect to distribution of these mitochondrial coastal and inland clades suggest some kind of climatic association for each clade. In Anatolia, coastal areas exhibit a more humid climate, in comparison to more steppe-like and drastically drier climates inland . Precipitation can influence vegetation, insect density and composition in a region . The distribution of the mitochondrial clades being closely associated with different rainfall intensities suggests that local adaptations to different habitats and climate regimes might have been a determinant of the differentiation of the species in the region. This kind of climatic association supports results from another congeneric species, M. natalensis in South Africa, which suggested an association between different biomes and population substructure .
In terms of other potential causal factors underlying this genetic break, none of the topographic barriers (the Marmara Sea, the Taurus Mountains or the eastern Anatolia Diagonal) seemed to be prominent impediments to dispersal. Individuals belonging to either clade S or clade P were found on either side of a putative barrier. Although AMOVA analysis showed some significant differences between regions as defined by geographic barriers (Φst = 0.27), these differences were not as obvious as those seen when comparing clades S and P with each other (Φst = 0.92). Also, the geographic differences did not hold up in all of the comparisons of the regions. For instance, significant differences were found between southern Anatolia and western Anatolia. However, no significant differences were found with comparisons of southern Anatolia and eastern Anatolia, even though Taurus Mountains delimits the border for southern Anatolia with both of the other regions. Hence although there is some genetic structure within the clades, the cause of this is probably not geographic barriers. The study of Bilgin et al.  had suggested that the Marmara and the Taurus were not likely to be impediments to dispersal, and the results of this study suggest that, in addition, the eastern Anatolian Diagonal does not seem to limit gene flow. Although some studies indicate that physical barriers can impede gene flow in bats [13, 61, 62], the results of this study also suggest that similar to two other bat species investigated in this region (Rhinolophus euryale , and Myotis capacinnii ), geographic barriers do not seem to limit gene flow for M. schreibersii.
Scenarios of post-glacial expansion
Further analysis gave additional support to the idea of postglacial expansions as the best explanation for the mtDNA differentiation. Statistical parsimony networks for both clade S and clade P showed star-like phylogenies, which are theoretically indicative of population expansions . For both of these clades, this pattern was supported by the R2 and F s statistics, the mismatch distributions of haplotypes and Bayesian skyline plots. Especially in clade S, certain haplotypes (S15, S3, S33 and S17) potentially represent nodes indicative of multiple bursts of population expansion.
The phylogenetic analyses showed no geographically associated structure within clade S or P. However, nested clade analysis  of the statistical parsimony networks showed subtle differentiation that delineated the details of possible historical dispersal routes. The nested cladograms for clade S showed four different second level nested clades (Figure 1b). The most ancestral haplotype, S15, showed a distribution almost exclusively restricted to the Balkans (Figure 1a). The central haplotype (S3) of the (light blue) nested clade 2.3, budding off from the ancestral nested clade 2.1 suggested a founder event from the Balkans to the Mediterranean coast and a subsequent range expansion of the (red) nested clade 2.4 along that coast. The (green) nested clade 2.2, which also derived out of the nested clade 2.1, showed a distribution predominantly restricted to the Black Sea coast, suggesting another migration route parallel to that coast. These results show the restriction of the most ancestral haplotype to the Balkans, followed by expansion and differentiation of different nested clades along the Black Sea and the Mediterranean coasts. This kind of post-glacial coastal expansion has been documented in relatively few taxa, with classical cases including human expansions such as colonization of eastern China  and western North America .
Using a Bayesian skyline plot approach, the timing of the onset of expansion of the entire clades S and P were found to be around 15,500 years ago and 4,500 years ago, respectively. These limits roughly correspond to the onset of deglaciation of continental ice sheets, after the last Heinrich event, around 15,000 years BP , and the stabilization of climate in the Holocene about 4,000 years BP to near today's temperature . This suggests that the final range expansions of clades S and P occurred after the end of Pleistocene, in a time range when the climate started to get warmer. The initiation of the expansion of these clades corresponding to post-Pleistocene is supported by evidence from other species, which show similar post Pleistocene modes of population expansion .
This type of distribution is expected from a founder event matching a leptokurtic model of migration (with a greater number of long distance migrants) where the populations behind the initial migration front do not advance as readily [15, 67]. NCA confirmed this with an explanation of contiguous and long-range expansion along the Black Sea and Mediterranean coasts respectively. This expansion pattern is also similar to that represented by the "grasshopper paradigm" where central and western Europe have been colonized from a Balkans' refuge, when the Iberian and Italian refugial populations were slowed down in the Pyrenees and the Alps . However, in this case, rather than by another refugial population, the expansion of the most ancestral haplotypes seems to have been blocked by the founder population that had colonized the uninhabited regions.
Using the Bayesian skyline approach for the nested clades, the onset of expansion for nested clades 2.2 and 2.3 were calculated to have started approximately at 4,000 and 6,000 years ago, respectively. Assuming that S15 was in the Balkans when the other clades started to expand, these results imply that this single haplotype has been restricted to the Balkans for at least 6,000 years. This roughly corresponds to 1000 generations, assuming an average generation time of six years for M. schreibersii (DEH, Australia). The Balkans has been proposed to be a glacial refugium for many taxa , including bats . Subsequent migrations of the other nested clades from the Balkans also seemed to follow geographic regions, such as the Black Sea and the Mediterranean with slightly different vegetation types and climate, and again suggest regional and climatic associations following migration events. This represents a form of extensive and long-term restriction of haplotypes and lineages to specific areas, and indicates that individuals belonging to the same matrilines have inhabited their respective geographic regions (Balkans, Black Sea coast, Mediterranean coast and inland) for hundreds of generations. Adaptations to local geographic conditions and climatic regimes can confer advantages to individuals, in terms of habituation to distribution and availability of roosts and resources , increased reproductive success, survival rate, reproductive output and recruitment . The results of this study suggest that these adaptations and preferences might be conserved over hundreds of generations, affecting the intraspecific microevolution of a species and the distribution of its resulting diversity.
Through combined analyses of the effects of various factors to the genetic differentiation of M. schreibersii, we were able to outline some of the details of the evolutionary history of this species in southeastern Europe and Anatolia. There was no evidence available that supported topographic barriers as determinants of genetic differentiation in this species. The results indicated mitochondrial genetic differentiation in glacial refugia, and subsequent range expansions coupled with various types of regional climate patterns. There was also evidence for long-term restriction of matrilines to these climatic regions, for hundreds of generations. Consequently, we were able to show how multiple biotic and abiotic events including glacial periods, climatic associations and historical dispersal patterns complemented each other in causing regional and local differentiation within a species. The next step along the lines of this research should include sampling in the winter, to see whether the observed genetic structure is maintained in the hibernation colonies or not.
Flannery TF: Mammals of the South-West Pacific and Moluccan Islands. 1995, Australia , Reed Books
Conservation International: Hotspots Revisited: Earth's Biologically Richest and Most Threatened Terrestrial Ecoregions. 2005
Akgündüz S: Türkiye'de Yağış , Sıcaklık ve Nem Verilerinin Klimatolojik Analizi. 2000, Ankara , DMİ Yayınları
Demirsoy A: Genel ve Türkiye Zoocoğrafyası. 1996, Ankara , Meteksan AŞ
Horacek I, Hanak V, Gaisler J: Bats of Palearctic region: A taxonomic and biogeographic review. Proceedings of the 8th EBRS. 2002, 1: 11-157.
Veith M, Schmidtler JF, Kosuch J, Baran I, Seitz A: Palaeoclimatic changes explain Anatolian mountain frog evolution: a test for alternating vicariance and dispersal events. Mol Ecol. 2003, 12 (1): 185-199. 10.1046/j.1365-294X.2003.01714.x.
Graham RW, Lundelius EL, Graham MA, Schroeder EK, Toomey RS, Anderson E, Barnosky AD, Burns JA, Churcher CS, Grayson DK, Guthrie RD, Harington CR, Jefferson GT, Martin LD, McDonald HG, Morlan RE, Semken HA, Webb SD, Werdelin L, Wilson MC: Spatial Response of Mammals to Late Quaternary Environmental Fluctuations. Science. 1996, 272 (5268): 1601-1606. 10.1126/science.272.5268.1601.
Koopman KF: Chiroptera: Systematics. Handb. Zool., 8 Mammalia, Pt. 60. 1994, New York , Walter de Gruyter, 217-
Simmons NB: Order Chiroptera. Mammal Species of the World: a taxonomic and geographic reference, third edition, vol 1. Edited by: Wilson DE, Reeder DM. 2005, Johns Hopkins University Press, 312-529.
Miller-Butterworth CM, Jacobs DS, Harley EH: Strong population substructure is correlated with morphology and ecology in a migratory bat. Nature. 2003, 424 (6945): 187-191. 10.1038/nature01742.
Cardinal BR, Christidis L: Mitochondrial DNA and morphology reveal three geographically distinct lineages of the large bentwing bat (Miniopterus schreibersii) in Australia. Aust J Zool. 2000, 48: 1-19. 10.1071/ZO99067.
Bilgin R, Karataş A, Çoraman E, Pandurski I, Papadatou E, Morales JC: Molecular Taxonomy and Phylogeography of Miniopterus schreibersii (Kuhl, 1817) (Chiroptera: Vespertilionidae), in the Eurasian Transition. Biol J Linn Soc. 2006, 87: 577-582. 10.1111/j.1095-8312.2006.00591.x.
Rossiter SJ, Benda P, Dietz C, Zhang S, Jones G: Rangewide phylogeography in the greater horseshoe bat inferred from microsatellites: implications for population history, taxonomy and conservation. Mol Ecol. 2007, 16 (22): 4699-4714. 10.1111/j.1365-294X.2007.03546.x.
Duran FS: Büyük Atlas. 1997, Istanbul , Kanaat Yayinlari
Hewitt GM: Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999, 68: 87-112.
Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF: Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998, 7 (4): 453-464. 10.1046/j.1365-294x.1998.00289.x.
Dubey S, Zaitsev M, Cosson JF, Abdukadier A, Vogel P: Pliocene and Pleistocene diversification and multiple refugia in a Eurasian shrew (Crocidura suaveolens group). Mol Phylogenet Evol. 2006, 38: 635-647. 10.1016/j.ympev.2005.11.005.
Jaarola M, Searle JB: Phylogeography of field voles (Microtus agrestis) in Eurasia inferred from mitochondrial DNA sequences. Mol Ecol. 2002, 11 (12): 2613-2621. 10.1046/j.1365-294X.2002.01639.x.
Anthony ELP: Age determination in bats. Ecological and behavioral methods for the study of bats. Edited by: Kunz TH. 1988, Washington, D.C. , Smithsonian Institution Press, pp. 47-58.
Worthington Wilmer J, Barratt EM: A non-lethal method of tissue sampling for genetic studies of chiropterans. Bat Res News. 1996, 37: 1-3.
Wilkinson GS, Chapman AM: Length and sequence variation in evening bat D-loop mtDNA. Genetics. 1991, 128 (3): 607-617.
Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG: The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 1997, 25: 4876-4882. 10.1093/nar/25.24.4876.
Bayes MK, Smith KL, Alberts SC, Altmann J, Bruford MW: Testing the reliability of microsatellite typing from faecal DNA in the savannah baboon. Conserv Genet. 2000, 1: 173-176. 10.1023/A:1026595324974.
Nei M, Tajima F: DNA polymorphism detectable by restriction endonucleases. Genetics. 1981, 97 (1): 145-163.
Graur D, Li WH: Fundamentals of molecular evolution. 2000, Sunderland, Mass. , Sinauer Associates, xiv, 481 p.-2nd
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.
Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony (and Other Methods). Version 4.06b. 2001, Massachusetts , Sinauer Associates
Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.
Hasegawa M, Kishino H, Yano T: Dating the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985, 22 (160:174): 678-687.
Yang Z: Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J Mol Evol. 1994, 39 (3): 306-314. 10.1007/BF00160154.
Kitching IJ, Forey PL, Humphries CJ, Williams D: Cladistics: The Theory and Practice of Parsimony Analysis. 2000, New York , Oxford University Press
Xia X, Xie Z: DAMBE: software package for data analysis in molecular biology and evolution. J Hered. 2001, 92 (4): 371-373. 10.1093/jhered/92.4.371.
Schneider S, Roessli D, Excoffier L: Arlequin: a Software Package for Population Genetics. 2000, Geneva , Genetics and Biometry Laboratory, Department of Anthropology. University of Geneva
Excoffier L, Smouse PE, Quattro JM: Analysis of Molecular Variance Inferred from Metric Distances among DNA Haplotypes - Application to Human Mitochondrial-DNA Restriction Data. Genetics. 1992, 131 (2): 479-491.
Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9 (10): 1657-1659. 10.1046/j.1365-294x.2000.01020.x.
Posada D, Crandall KA, Templeton AR: GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Mol Ecol. 2000, 9 (4): 487-488. 10.1046/j.1365-294x.2000.00887.x.
Ramos-Onsins SE, Rozas J: Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002, 19 (12): 2092-2100.
Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147 (2): 915-925.
Drummond AJ, Rambaut A, Shapiro B, Pybus OG: Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005, 22 (5): 1185-1192-10.1093/molbev/msi103.
Goudet J: FSTAT, Version 1.2.a program for IBM PC compatibles to calculate Weir and Cockerham's estimators of F-statistics. J Hered. 1995, 86: 485-486.
Rice WR: Analyzing tables of statistical tests. Evolution. 1989, 43: 223-225. 10.2307/2409177.
van Oosterhout C, Hutchinson WF, Wills DPM, Shipley PF: Micro-Checker: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004, 4: 535–538-10.1111/j.1471-8286.2004.00684.x.
Leberg PL: Estimating allelic richness: effects of sample size and bottlenecks. Mol Ecol. 2002, 11 (11): 2445-2449. 10.1046/j.1365-294X.2002.01612.x.
Peakall R, Smouse PE: GenAlEx 6: Genetic Analysis in Excel. Population genetic software for teaching and research. 2005, Canberra, Australia , The Australian National University
Kalinowski ST: HP-Rare: a computer program for performing rarefaction on measures of allelic diversity. Mol Ecol Notes. 2005, 5: 187-189. 10.1111/j.1471-8286.2004.00845.x.
Michalakis Y, Excoffier L: A generic estimation of population subdivision using distances between alleles with special reference for microsatellite loci. Genetics. 1996, 142 (3): 1061-1064.
Goodman SJ: R-ST Calc: a collection of computer programs for calculating estimates of genetic differentiation from microsatellite data and determining their significance. Mol Ecol. 1997, 6 (9): 881-885. 10.1111/j.1365-294X.1997.tb00143.x.
Hardy OJ, Charbonnel N, Freville H, Heuertz M: Microsatellite allele sizes: a simple test to assess their significance on genetic differentiation. Genetics. 2003, 163 (4): 1467-1482.
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.
Nachman MW, Boyer SN, Searle JB, Aquadro CF: Mitochondrial DNA variation and the evolution of Robertsonian chromosomal races of house mice, Mus domesticus. Genetics. 1994, 136 (3): 1105-1120.
Shapiro B, Drummond AJ, Rambaut A, Wilson MC, Matheus PE, Sher AV, Pybus OG, Gilbert MT, Barnes I, Binladen J, Willerslev E, Hansen AJ, Baryshnikov GF, Burns JA, Davydov S, Driver JC, Froese DG, Harington CR, Keddie G, Kosintsev P, Kunz ML, Martin LD, Stephenson RO, Storer J, Tedford R, Zimov S, Cooper A: Rise and fall of the Beringian steppe bison. Science. 2004, 306 (5701): 1561-1565. 10.1126/science.1101074.
Hewitt GM: Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999, 68: 87-112.
Bilgin R, A. F, Çoraman E, Karataş A: Phylogeography of the Mediterranean horseshoe bat, Rhinolophus euryale (Chiroptera: Rhinolophidae), in southeastern Europe and Anatolia. Acta Chiropt. 2008, 10: 41-49. 10.3161/150811008X331072.
Bilgin R, Karataş A, Çoraman E, Morales JC: The mitochondrial and nuclear genetic structure of Myotis capaccinii (Chiroptera: Vespertilionidae) in the Eurasian transition, and its taxonomic implications. Zool Scr. 2008, 37 (3): 253-262. 10.1111/j.1463-6409.2008.00326.x.
Kerth G, Petrov B, Conti A, Anastasov D, Weishaar M, Gazaryan S, Jaquiery J, Konig B, Perrin N, Bruyndonckx N: Communally breeding Bechstein's bats have a stable social system that is independent from the postglacial history and location of the populations. Mol Ecol. 2008, 17 (10): 2368-2381. 10.1111/j.1365-294X.2008.03768.x.
Emmons LH, Leite YLR, Kock D, Costa LP: A Review of the Named Forms of Phyllomys (Rodentia: Echimyidae) with the Description of a New Species from Coastal Brazil. Am Mus Novit. 2002, 3380: 1-40. 10.1206/0003-0082(2002)380<0001:AROTNF>2.0.CO;2.
Wilkinson GS, Fleming TH: Migration and evolution of lesser long-nosed bats Leptonycteris curasoae, inferred from mitochondrial DNA. Mol Ecol. 1996, 5 (3): 329-339. 10.1111/j.1365-294X.1996.tb00324.x.
Greenberg R, Cordero PJ, Droege S, Fleischer RC: Morphological adaptation with no mitochondrial DNA differentiation in the coastal plain swamp swallow. The Auk. 1998, 115 (3): 706-712.
White EE: Chloroplast DNA in Pinus monticola 2. Survey of within-species variability and detection of heteroplasmic individuals. Theor Appl Genet. 1990, 79: 251-255.
Wallner WE: Factors Affecting Insect Population Dynamics: Differences Between Outbreak and Non-Outbreak Species. Annu Rev Entomol. 1987, 32: 317-340. 10.1146/annurev.en.32.010187.001533.
Ruedi M, Castella V: Genetic consequences of the ice ages on nurseries of the bat Myotis myotis: a mitochondrial and nuclear survey. Mol Ecol. 2003, 12 (6): 1527-1540. 10.1046/j.1365-294X.2003.01828.x.
Kerth G, Petit E: Colonization and dispersal in a social species, the Bechstein's bat (Myotis bechsteinii). Mol Ecol. 2005, 14 (13): 3943-3950. 10.1111/j.1365-294X.2005.02719.x.
Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history. Mol Ecol. 1998, 7 (4): 381-397. 10.1046/j.1365-294x.1998.00308.x.
Eshleman JA, Malhi RS, Johnson JR, Kaestle FA, Lorenz J, Smith DG: Mitochondrial DNA and prehistoric settlements: native migrations on the western edge of North America. Hum Biol. 2004, 76 (1): 55-75. 10.1353/hub.2004.0019.
Tanaka M, Cabrera VM, Gonzalez AM, Larruga JM, Takeyasu T, Fuku N, Guo LJ, Hirose R, Fujita Y, Kurata M, Shinoda K, Umetsu K, Yamada Y, Oshida Y, Sato Y, Hattori N, Mizuno Y, Arai Y, Hirose N, Ohta S, Ogawa O, Tanaka Y, Kawamori R, Shamoto-Nagai M, Maruyama W, Shimokata H, Suzuki R, Shimodaira H: Mitochondrial genome variation in eastern Asia and the peopling of Japan. Genome Res. 2004, 14 (10A): 1832-1850. 10.1101/gr.2286304.
Walker MJC: Climatic changes in Europe during the last glacial/interglacial transition. Quat Int. 1995, 28: 63-76. 10.1016/1040-6182(95)00030-M.
Ibrahim KM, Nichols JD, Hewitt GM: Spatial patterns of genetic variation generated by different forms of dispersal during range expansion. Heredity. 1996, 77: 282-291. 10.1038/hdy.1996.142.
Lewis SE: Roost Fidelity of Bats: A Review. J Mammal. 1995, 76: 481-496. 10.2307/1382357.
Franklin AB, Anderson DR, Gutierrez RJ, Burnham KP: Climate, Habitat Quality, and Fitness in Northern Spotted Owl Populations in Northwestern California. Ecol Monogr. 2000, 70: 539-590.
We would like to thank Ivan Pandurski, Eleni Papadatou, Ayşegül Karataş, Hasan Karakaya and Ferhat Toprak for their assistance with obtaining the samples. We also would like to thank Elizabeth M. Hemond, Don J. Melnick, Sergios Orestis Kolokotronis and two anonymous reviewers for their comments on earlier versions of the manuscript, and Boğaziçi University Speleological Society (BÜMAK) for logistical support.
RB conceived the study, carried out the fieldwork and molecular genetic studies, performed statistical analyses and drafted the manuscript. AK coordinated the fieldwork and helped in drafting the manuscript. EÇ participated in the design of the study and performed statistical analyses. TD supervised the laboratory components of the study and helped in drafting the manuscript. JCM supervised the entire research. All authors read and approved the final manuscript.
About this article
Cite this article
Bilgin, R., Karataş, A., Çoraman, E. et al. Regionally and climatically restricted patterns of distribution of genetic diversity in a migratory bat species, Miniopterus schreibersii(Chiroptera: Vespertilionidae). BMC Evol Biol 8, 209 (2008) doi:10.1186/1471-2148-8-209
- Genetic Differentiation
- Suture Zone
- Private Allele
- Mediterranean Coast
- Ancestral Haplotype