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.
The names of the sampling locations, their ID's as they appear in Fig 1., the geographical coordinates, sample size per locality (N), number of haplotypes, the names of the haplotypes as they appear in Figures 1 and 2, and the numbers of males and females at each site.
No. of Haplotypes
Haplotype names (numbers)
No. of Males (M) and Females (F)
S6, S15(3), S24, S25, S29(2), S31, S39
8 M, 2 F
S1, S5. S23. S25 (2)
2 M, 3 F
S15(4), S25, S26
4 M, 2 F
1 M, 2 F
S15, S16, S39, S40
1 M, 2 F
S15(4), S25, S26
4 M, 2 F
S14, S15(2), S26
3 M, 1 F
S2(2), S17(5), S18
4 M, 4 F
10 M, 4 F
S8, S11(2), S21, S32, S33(8), S36
10 M, 4 F
2 M, 3 F
P2(3), P4(3), P5(1), P6(1)
2 M, 6 F
P2(2), P6(2), P7, P8, P9, P10
1 M, 7 F
The localities where samples from Bilgin et al.  were used in the current study are shown in bold. From Kemah, only one individual was sampled in Bilgin et al. .
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 FS in DnaSP. Based on simulation models, among various statistics that determine signatures of past population growth, R2 and FSare 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%.
Percent differentiation and the number of fixed differences (in parentheses) between clades S, P, and the Indonesian sample.
The diagonal comparisons indicate within-group differentiation.
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).
Percent distribution of variation when individuals are grouped as in clades S and P.
Amoung Populations within Groups
Percent distribution of variation when individuals are grouped by geographical barriers.
Amoung Populations within Groups
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 FS. 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 FS, 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.
Unpaired t-test results for climate differences (summer precipitation) between clades S and P.
Mean summer precipitation (mm)
tcrit (df = 115)
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.
The results of global and locus-by locus AMOVAs, Fst and Rst
Within Pop. Variance
Among Pop. Variance
Within Pop. Variance
Among Pop. Variance
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 Fsstatistics, 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.
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.
Institute of Environmental Sciences, Boğaziçi University
Department of Ecology Evolution and Environmental Biology, Columbia University
Center for Environmental Research and Conservation, Columbia University
Zübeyde Hanım Sağlık Yüksekokulu, Niğde Üniversitesi
Department of Anthropology, New York University
Stanislaus, One University Circle, California State University
Flannery TF: Mammals of the South–West Pacific and Moluccan Islands. Australia, Reed Books 1995.
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. Ankara, DMİ Yayınları 2000.
Demirsoy A: Genel ve Türkiye Zoocoğrafyası. Ankara, Meteksan AŞ 1996.
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.View ArticlePubMed
Graham RW, Lundelius EL Jr., Graham MA, Schroeder EK, Toomey RS III, Anderson E, Barnosky AD, Burns JA, Churcher CS, Grayson DK, Guthrie RD, Harington CR, Jefferson GT, Martin LD, McDonald HG, Morlan RE, Semken HA Jr., Webb SD, Werdelin L, Wilson MC: Spatial Response of Mammals to Late Quaternary Environmental Fluctuations.Science 1996, 272 (5268) : 1601–1606.View ArticlePubMed
Koopman KF: Chiroptera: Systematics. Handb. Zool., 8 Mammalia, Pt. 60. New York, Walter de Gruyter 1994, 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). Johns Hopkins University Press 2005, 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.View ArticlePubMed
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.View Article
Bilgin R, Karataş A, Çoraman E, Pandurski I, Papadatou E, Morales JC: Molecular Taxonomy and Phylogeography ofMiniopterus schreibersii(Kuhl, 1817) (Chiroptera: Vespertilionidae), in the Eurasian Transition.Biol J Linn Soc 2006, 87: 577–582.View Article
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.View ArticlePubMed
Duran FS: Büyük Atlas. Istanbul, Kanaat Yayinlari 1997.
Hewitt GM: Post–glacial re–colonization of European biota.Biol J Linn Soc 1999, 68: 87–112.View Article
Dubey S, Zaitsev M, Cosson JF, Abdukadier A, Vogel P: Pliocene and Pleistocene diversification and multiple refugia in a Eurasian shrew (Crocidura suaveolensgroup).Mol Phylogenet Evol 2006, 38: 635–647.View ArticlePubMed
Jaarola M, Searle JB: Phylogeography of field voles (Microtus agrestis) in Eurasia inferred from mitochondrial DNA sequences.Mol Ecol 2002, 11 (12) : 2613–2621.View ArticlePubMed
Anthony ELP: Age determination in bats.Ecological and behavioral methods for the study of bats(Edited by: Kunz TH). Washington, D.C., Smithsonian Institution Press 1988, 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.PubMed
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.View ArticlePubMed
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.View Article
Nei M, Tajima F: DNA polymorphism detectable by restriction endonucleases.Genetics 1981, 97 (1) : 145–163.PubMed
Graur D, Li WH: Fundamentals of molecular evolution.2nd Edition Sunderland, Mass., Sinauer Associates 2000, xiv, 481 p..
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.View ArticlePubMed
Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony (and Other Methods). Version 4.06b. Massachusetts, Sinauer Associates 2001.
Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution.Bioinformatics 1998, 14: 817–818.View ArticlePubMed
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.View ArticlePubMed
Kitching IJ, Forey PL, Humphries CJ, Williams D: Cladistics: The Theory and Practice of Parsimony Analysis. New York, Oxford University Press 2000.
Xia X, Xie Z: DAMBE: software package for data analysis in molecular biology and evolution.J Hered 2001, 92 (4) : 371–373.View ArticlePubMed
Schneider S, Roessli D, Excoffier L: Arlequin: a Software Package for Population Genetics. Geneva, Genetics and Biometry Laboratory, Department of Anthropology. University of Geneva 2000.
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.PubMed
Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies.Mol Ecol 2000, 9 (10) : 1657–1659.View ArticlePubMed
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.View ArticlePubMed
Ramos–Onsins SE, Rozas J: Statistical properties of new neutrality tests against population growth.Mol Biol Evol 2002, 19 (12) : 2092–2100.PubMed
Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection.Genetics 1997, 147 (2) : 915–925.PubMed
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.View ArticlePubMed
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.
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.
Leberg PL: Estimating allelic richness: effects of sample size and bottlenecks.Mol Ecol 2002, 11 (11) : 2445–2449.View ArticlePubMed
Peakall R, Smouse PE: GenAlEx 6: Genetic Analysis in Excel. Population genetic software for teaching and research. Canberra, Australia, The Australian National University 2005.
Kalinowski ST: HP–Rare: a computer program for performing rarefaction on measures of allelic diversity.Mol Ecol Notes 2005, 5: 187–189.View Article
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.PubMed
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.View Article
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.PubMed
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.View Article
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.PubMed
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.View ArticlePubMed
Hewitt GM: Post–glacial re–colonization of European biota.Biol J Linn Soc 1999, 68: 87–112.View Article
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.View Article
Bilgin R, Karataş A, Çoraman E, Morales JC: The mitochondrial and nuclear genetic structure ofMyotis capaccinii(Chiroptera: Vespertilionidae) in the Eurasian transition, and its taxonomic implications.Zool Scr 2008, 37 (3) : 253–262.View Article
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.View ArticlePubMed
Emmons LH, Leite YLR, Kock D, Costa LP: A Review of the Named Forms ofPhyllomys(Rodentia: Echimyidae) with the Description of a New Species from Coastal Brazil.Am Mus Novit 2002, 3380: 1–40.View Article
Wilkinson GS, Fleming TH: Migration and evolution of lesser long–nosed batsLeptonycteris curasoae, inferred from mitochondrial DNA.Mol Ecol 1996, 5 (3) : 329–339.PubMed
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.View Article
Ruedi M, Castella V: Genetic consequences of the ice ages on nurseries of the batMyotis myotis: a mitochondrial and nuclear survey.Mol Ecol 2003, 12 (6) : 1527–1540.View ArticlePubMed
Kerth G, Petit E: Colonization and dispersal in a social species, the Bechstein's bat (Myotis bechsteinii).Mol Ecol 2005, 14 (13) : 3943–3950.View ArticlePubMed
Templeton AR: Nested clade analyses of phylogeographic data: testing hypotheses about gene flow and population history.Mol Ecol 1998, 7 (4) : 381–397.View ArticlePubMed
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.View ArticlePubMed
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.View ArticlePubMed
Walker MJC: Climatic changes in Europe during the last glacial/interglacial transition.Quat Int 1995, 28: 63–76.View Article
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.View Article
Lewis SE: Roost Fidelity of Bats: A Review.J Mammal 1995, 76: 481–496.View Article
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.View Article
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.