Between the Balkans and the Baltic: Phylogeography of a Common Vole Mitochondrial DNA Lineage Limited to Central Europe

The common vole (Microtus arvalis) has been a model species of small mammal for studying end-glacial colonization history. In the present study we expanded the sampling from central and eastern Europe, analyzing contemporary genetic structure to identify the role of a potential ‘northern glacial refugium’, i.e. a refugium at a higher latitude than the traditional Mediterranean refugia. Altogether we analyzed 786 cytochrome b (cytb) sequences (representing mitochondrial DNA; mtDNA) from the whole of Europe, adding 177 new sequences from central and eastern Europe, and we conducted analyses on eight microsatellite loci for 499 individuals (representing nuclear DNA) from central and eastern Europe, adding data on 311 new specimens. Our new data fill gaps in the vicinity of the Carpathian Mountains, the potential northern refugium, such that there is now dense sampling from the Balkans to the Baltic Sea. Here we present evidence that the Eastern mtDNA lineage of the common vole was present in the vicinity of this Carpathian refugium during the Last Glacial Maximum and the Younger Dryas. The Eastern lineage expanded from this refugium to the Baltic and shows low cytb nucleotide diversity in those most northerly parts of the distribution. Analyses of microsatellites revealed a similar pattern but also showed little differentiation between all of the populations sampled in central and eastern Europe.


Introduction
Phylogeography is a discipline which uses molecular markers to infer population movements at the end of the Pleistocene and during the Holocene, to identify where glacial refugia were  [24,[28][29][30]. The brown area shows the range of the species in Eurasia. The solid line marks the hybrid zone between the western (arvalis) and eastern (obscurus) forms; only the data for the western form are used in this study. Sampling localities are colored according to the cytochrome b lineage. (B) Distribution of common vole cytochrome b samples in central Europe belonging to the Eastern mtDNA lineage, both newly sequenced and from GenBank. The blue circles match the location of the Balkan mtDNA lineage, which makes contact with the Eastern lineage to the south. The yellow circles show the occurrence of the Central mtDNA lineage to the north of the area occupied by the Eastern lineage, including where they make contact. The black stars mark the locations where Microtus rossiaemeridionalis was found among the samples sequenced (KX380179-KX380191). Location numbers on the map are listed in S1, S2 and S3 Tables. doi:10.1371/journal.pone.0168621.g001

Phylogeography of a Common Vole Lineage
Specifically, this broad sampling allows us to consider the following questions relating to occupation of northern glacial refugia by the common vole: (i) Is the pattern of mitochondrial and microsatellite variation consistent with the Eastern lineage deriving from the Carpathian Basin, as opposed to a refugium or refugia located elsewhere? (ii) How does the timing of expansion of the Eastern lineage relate to events during the last glaciation, particularly the LGM? (iii) Is the differentiation that is evident in Poland, based on microsatellite variation, replicated at a wider scale in the Eastern lineage and what is the implication for the refugial population and its subsequent expansion? Overall, the level and density of sampling in the area surrounding this northern refugium provides an exceptionally detailed perspective on the relationship between genetic variation in a contemporary population and the refugium, or refugia, from which it originated.

Materials and Methods
All capture and handling procedures of voles in Poland were approved by the Local Ethics Committee for Animal Experiments in Białystok (Decision no: 5/2013). Samples of voles from Croatia, Czech Republic, Hungary, Romania, Serbia and Slovenia were collected according to the Legislation of Nature Protection and Legislation of Animal Welfare valid in these countries and permission to the Slovenian Museum of Natural History given by the Veterinary Administration of the Republic of Slovenia. No additional permits from Ethics Committee were needed. All capture and handling procedures of the common voles from Russia and Ukraine were approved by the Animal Care and Use Committee of the A.N. Severtsov Institute of Ecology and Evolution of the Russian Academy of Sciences (official letter from the deputy director of the Institute, no: 3485/09). No additional permits were required. Permission from the landowners was obtained where necessary. Samples were not collected from protected locations and field studies did not involve endangered, threatened or protected species.

Mitochondrial DNA
We obtained new mtDNA data from 177 samples of the common vole from Poland, Belarus, Russia, Moldova, Ukraine, Czech Republic, Hungary, Slovenia, Croatia, Romania, and Serbia (Fig 1, S1 Table). Samples were collected between 1995 and 2015 and preserved in 96-98% ethanol. Additionally, our collection included samples from central and southern Ukraine and Armenia but all of them appeared to be the southern vole (Microtus rossiaemeridionalis), a sibling species almost identical in morphology to the common vole [35]. We obtained 13 complete cytb sequences from those samples and they were identified using BLAST (https:// blast.ncbi.nlm.nih.gov/Blast.cgi). They were not included further in the analyses performed in this study. The sequences obtained were deposited in GenBank (accession numbers: KX380179-KX380191).
Total genomic DNA was extracted using Syngen Tissue DNA Mini Kit in accordance with the manufacturer's instructions. The amplification and sequencing of the complete cytb sequence (1143 bp) were conducted according to methods described by Wójcik et al. [10], using primers described by Stojak et al. [30], in two overlapping fragments, each 500-600 bp in length. We also included negative PCR controls without template DNA. Additionally, we used 609 sequences obtained from previous studies (Fig 1A; for the whole list of locations see S1, S2 and S4 Tables in [30]). In all, we analyzed cytb sequences from 786 individuals (Fig 1A,  S1 and S2 Tables).
The nucleotide sequences were aligned and manually checked in BioEdit 7.2.5 [36]. As in our previous study [30], the complete alignment was shortened to 1000 bp. The distribution of polymorphic nucleotide and amino acid sites, number of transversions and transitions, synonymous and non-synonymous changes and position of stop codons in cytb sequences obtained from newly-collected individuals were determined in DnaSP 5.1 [37] using the method of Degli Esposti et al. [38].
A median-joining network was constructed for all 786 sequences in Network 4.6 [39], to establish the overall pattern of mitochondrial genetic variation in the species.
We used Bayesian genealogy sampling in BEAST 1.8 [40], with all 786 sequences, to estimate the time to most recent common ancestor (tMRCA) for each lineage. The analysis was carried out according to the procedure described by Stojak et al. [30], using the HKY+Γ substitution model. A strict molecular clock was compared with an uncorrelated lognormal relaxed molecular clock [41]. Two demographic models were tested: skyline and constant population size [42]. For model selection, path sampling and stepping-stone sampling [43] were used to estimate marginal likelihoods (MLEs) for each model, based on four independent MCMC chains, each comprising 1000 steps of 100 000 generations, following 10 million generations burn-in. In each case, sampling was considered sufficient, as the estimates from the two methods (path sampling and stepping-stone sampling) converged on similar values. The MLEs were then used to obtain the Bayes Factor for each comparison to determine the best-fitting model [44]. The previously estimated clock rate of 3.27 x 10 −7 substitutions/site/year -1 was used to calibrate the genealogy and date tMRCAs of each lineage [28,30]. Posterior distributions of tMRCA and other model parameters were obtained from four independent MCMC simulations (each run for 200 million generations, with effective sample size for each parameter over 200), using TRACER 1.5 [45]. Posterior samples from the four chains were combined, using LogCombiner, part of the BEAST package. The Maximum Clade Credibility (MCC) tree was obtained using TreeAnnotator, part of the BEAST package, and visualized with FigTree 1.4 (http://tree-bio.ed.ac.uk/). The analyses were repeated without sequence data to test the effect of the joint prior distributions on the posterior distributions for each parameter of interest.
Nucleotide and haplotype diversity were calculated in DnaSP. To test for recent population expansion in mitochondrial lineages, we applied three tests in Arlequin 3.5.1.2 [46]: Tajima's [47] D and Fu's [48] F S statistics and the mismatch distribution using sum of squared deviations SSD [49], with significance inferred using 1000 bootstrap replicates. Subsequently, we calculated time since expansion for lineages in which all three tests confirmed recent expansion (Tajima's D and Fu's F S significant, SSD not significant; see [16,30,50]). The estimate was calculated according to a method based on the Tau value (Tau = 2ut, where u is the mutation rate and t is the mean generation time [51,52]) using an online tool described by Schenekar and Weiss [53] (mismatch calculator, kindly provided by Stephen Weiss). Tau values were calculated in Arlequin, the substitution rate of 3.27 x 10 -7 substitutions/site/year -1 was applied and the generation time was assumed to be one year [28,30]. To examine the demographic expansion of the Eastern lineage, a Bayesian skyline plot was prepared in TRACER, according to the method described in [30]. In this analysis we used all 318 sequences belonging to this lineage.
Geographic variation in nucleotide diversity (π) for the Eastern lineage was mapped with interpolation using ArcGIS 10. . Populations with only one individual were excluded to avoid cases where π is equal to 0. Altogether we used 293 sequences from 74 populations. Interpolation was carried out using an Inverse Distance Weight model (IDW, power = 1, based on 12 neighbors) which uses a linear-weighted combination set of sample points, assigned as a function of the distance of an input point from the output cell location. This means that the greater the distance is, the less influence the input cell obtained has on the results.

Microsatellite DNA
We used 311 new specimens for microsatellite analysis, obtained from the Carpathian Basin, Russia and the Balkans, together with 188 specimens which were already analyzed in the previous study of Stojak et al. [34]. We initially used the previous panel of 11 microsatellite loci as described by Stojak et al. [34], following the same protocols. However, this was reduced for this study by the exclusion of loci AV12, Moe5 and MSM2, due to excessive missing data or risk in identification of false alleles as a result of non-specific amplification. Negative PCR controls without template DNA were included for each set of samples and each multiplex. We repeated the PCR for 15% of samples to test for error. All genotypes were independently scored three times. Altogether, we conducted the study on 499 individuals from 61 localities (Fig 1B). Sample sizes for populations varied, ranging from 1-16 individuals. A complete list of populations, localities and number of individuals are given in S3 Table. Only populations with five or more individuals were included in all population-level analyses. Observed (H O ) and expected heterozygosity (H E ) were calculated in Arlequin. FSTAT v. 2.9.3 [54] was used to determine deviations from Hardy-Weinberg Equilibrium (HWE), allelic richness (AR) and to test for linkage disequilibrium between all pairs of loci in each population (using 10,000 simulations). Null alleles were identified in Micro-Checker v. 2.2.3 [55]. Pairwise F ST values [56] between populations within each species were estimated in FSTAT, determining significance using 10,000 permutations under sequential Bonferroni corrections [57]. The relationship between pairwise genetic (F ST ) and geographic distances was tested using a Mantel test [58] in IBDWS v. 3.23 [59].
To establish the number of genetic clusters occurring across central Europe we initially used Bayesian analysis in STRUCTURE 2.3.4 [60]. Ten independent runs of 500,000 generations with 100,000 generations of burn-in were performed for each value of K (1-20) under the admixture model and with the assumption that allele frequencies among populations are correlated. The optimal number of clusters were identified using the Evanno method (ΔK; [61]) and mean posterior probability of the data (Ln(K); log probability). Calculations were conducted in STRUCTURE HARVESTER [62]. We assigned populations to a particular cluster according to assignment probability q ! 0.8. Populations with assignment probability of 0.2 < q < 0.8 were classified as 'intermediate' [34,63]. We used Clumpp 1.1.2 [64] and Distruct 1.1 [65] to assemble 10 independent runs of best-fitting values of K and to generate graphical representations.
Additionally, we performed Spatial Principal Components Analysis (sPCA; [66]) in R studio [67], to describe genetic differentiation and identify potential genetic structure across all populations from central Europe. Within that package, a G test was used to check for global (overall) structure and an L test for local (smaller scale) structure. Both tests were conducted using 999 permutations. The algorithm used in sPCA is modified such that the principal axes maximize spatial autocorrelation instead of correlation. Unlike methods such as STRUC-TURE, sPCA does not take into consideration Hardy-Weinberg Equilibrium (HWE) or the extent of linkage disequilibrium [66,68].
Interpolation of allelic richness (AR) and expected heterozygosity (H E ) were conducted for 42 populations to identify potential hotspots of genetic diversity to complement the nucleotide diversity analysis with cytb. Only populations with 5 or more individuals were included. Two Russian populations were excluded from this analysis because they were geographically distant from all the other populations. Analyses were carried out using ArcGIS 10.3.1, adopting the same model as for cytb.

Mitochondrial cytb sequences
Amplification of the complete cytb sequence (1143 bp) was successful for all 177 new samples of the common vole (S1 Table). We found no contamination in negative controls. The distribution of polymorphic nucleotide sites showed that substitutions are mostly in the third codon position, with a higher ratio of transitions (106Ts:13Tv; S4 Table). Non-synonymous changes occurred at the first codon position more often than at the second. Substitutions were mostly in the third codon position, with a higher ratio of transitions and only one (final) stop codon was identified (nucleotide positions 1141-1143). We found a much higher number of synonymous than non-synonymous changes (99: 20). These results are in agreement with a pattern that has been previously described by Irwin et al. [69] in mammals (and other taxa).
The 177 new sequences had 126 polymorphic sites, 85 of which were phylogenetically informative, and 97 haplotypes were identified. All 786 sequences analyzed contained 267 polymorphic sites, 181 of which were phylogenetically informative, and 309 haplotypes were identified. No nuclear mitochondrial DNA segments (numts) [70,71] were found. PCR amplifications always produced only one band and we found no ambiguities in the sequences obtained-there were no unexpected stop codons, deletions, insertions or changes in the reading frame; all sequences were in concordance with the rest of cytb sequences of the common vole from GenBank.
The network analysis allowed us to assign the new sequences to the six previously described mtDNA lineages [24,[28][29][30]72]. Eight new sequences from the Balkan region were assigned to the Balkan lineage, three from north-western Poland to the Central lineage and the remaining 166 to the Eastern lineage (Fig 2).
In the cytb genealogy, illustrated by the MCC tree and obtained from MCMC simulations with a coalescent model, all of the sequences were again included in the six previously identified lineages (S1 Fig). Using all 786 sequences, the marginal likelihood estimates (MLEs) for skyline and constant population size demographic models were -9474.7 and -9684.5, respectively. A coalescent skyline model was therefore chosen over a constant population size model (Bayes Factors BF: 209.8) and an uncorrelated relaxed lognormal molecular clock was chosen over a strict molecular clock (MLEs: -9474.7 and -9509.4 respectively; BF: 34.7). In the Bayesian genealogy, all six lineages were highly supported (posterior probabilities either 0.99 or 1.0). While the Balkan and Italian lineages were entirely distinct from all of the other lineages, the Western-North and Western-South lineages formed a distinct and well-supported clade (PP = 0.95), as did the Central and Eastern lineages (PP = 0.99). The deeper splits in the genealogy were not supported and the relationship between the lineages is therefore uncertain.
Genetic variation in the three lineages from eastern Europe, which are those of interest here, is summarized in Table 2 Table 1. Time to most recent common ancestor (tMRCA) for the common vole population and for the six cytochrome b lineages, with median and 95% highest posterior density (HPD) limits estimated from 720 million post-burn-in generations obtained from four independent MCMC simulations in BEAST.

95% HPD lower limit (kya)
Median (kya) 95% HPD upper limit (kya) The interpolation of nucleotide diversity in the Eastern lineage showed that this is highest in the Carpathian area and decreases northwards with increasing distance from this area (Fig 3).

Microsatellites
We analyzed eight microsatellite loci in 499 individuals of the common vole (S5 Table). There were only 0.33% missing data among all samples. No errors were found in the re-  Table). Only five populations significantly differed from HWE expectations after Bonferroni correction. Tests of linkage disequilibrium over all populations revealed significant linkage at all loci for two populations from Aleksinac (Serbia) and Česky Dub (Czech Republic).
Overall population differentiation was rather low (F ST = 0.040, 95% CI: 0.032-0.048) and this value is nearly identical after removing the Moe6 locus (F ST = 0.041, 95% CI: 0.031-0.050). There was significant differentiation between some pairs of populations (S7 Table). A Mantel test of the relationship between genetic and geographic distances revealed a weak and non-significant isolation by distance pattern for the common vole in central and eastern Europe (r = 0.1725, P = 0.068). This overall pattern of genetic differentiation is consistent with that from earlier studies of the common vole [31,34].
The values of Ln(K) and ΔK, obtained with STRUCTURE were not entirely consistent (S3 Fig). From Ln(K) there was support for K = 9, while from ΔK three values were supported (K = 2, 3 and 9). With K of two or three, there was a clear pattern of group membership within populations and this was reflected in the geographical distribution of the populations (Fig 4  and S4 Fig). In contrast, with K = 9 there was little consistency of group membership within populations (Fig 4).
In the sPCA analysis, only the highest eigenvalue is dominant (S5A and S5B Fig).  plots with spatial representation of values for each population revealed that genetic differentiation between them was rather low and no distinct separation was observed (S6 Fig).
The interpolations of both allelic richness (AR, Fig 5A) and expected heterozygosity (H E, Fig 5B) showed that diversity is highest in the north-eastern part of the Carpathian Basin and this diversity decreases gradually away from this area.

Discussion
The common vole Microtus arvalis is a particularly well-studied phylogeographic model species in Europe. A number of previous studies have revealed multiple mtDNA lineages of this species, based on cytb sequences, and the distribution of these lineages has been related to end-glacial biogeographic history [24,[28][29][30]72]. The use of other genetic markers, such as microsatellites or those on the Y-chromosome, have generally revealed patterns congruent to those found with mtDNA [28,31,73,74].
Here, we substantially enhanced the geographic coverage of cytb sequences of M. arvalis in central and eastern Europe. In our analysis of cytb, we made use of all 786 sequences available (including the samples in this paper and in GenBank) to confirm the six mtDNA lineages previously described in the species (Western-South, Western-North, Italian, Central, Eastern and Balkan; Fig 1A). Our material also included samples from central and southern Ukraine and Armenia but all of them appeared to be from the southern vole (Microtus rossiaemeridionalis), a sibling species almost identical in morphology to the common vole [35], suggesting that the range of the common vole does not extend into these areas.
While European temperate species were traditionally thought to have occupied Mediterranean glacial refugia, the possibility of more northerly refugia has been promoted since the 1980s [6,9,[11][12][13]16], however studies of present-day populations in the vicinity of these refugia and expansion routes from them have not been exceptionally detailed. The common vole lineages which may derive from such northern refugia are the Western-North lineage (with a refugium located somewhere in central France; [28]), the Central lineage (with a refugium located somewhere in the Alpine region; [73]) and the Eastern lineage (with a refugium located in the Carpathians; [30]). However, the status of any of these putative common vole refugia, as a source for contemporary genetic variation, has not been thoroughly examined. The Eastern lineage and the status of its proposed Carpathian refugium are the focus of the present study. The existence of northern refugia is very interesting in a physiological sense, as the survival of a temperate species in these refugia would suggest that it has high resistance to the severe, dry and sub-arctic climatic conditions present at such latitudes during the LGM and Younger Dryas [4,5,18,21,22]. Such populations may have been able to expand rapidly on amelioration of the climate after the end of the cold period. Northern refugia may therefore have made a substantial contribution to the contemporary genetic diversity of high latitudes in Europe [9,13,14].
We analyzed 318 cytb sequences from the Eastern lineage, with excellent coverage of its distribution between the Balkans and Baltic Sea. This area which would have been affected by severe glacial conditions during the LGM [4,5], but nevertheless includes the vicinity of the Carpathian Mountains, proposed as a refugium for a variety of animals and plants, including hornbeam (Carpinus betulus; [13]), crested newt (Triturus cristatus; [75]) and brown bear (Ursus arctos; [76]). The Carpathian region has a well-established fossil history, including remains of common voles and other small mammals, from the time of the LGM [32,33].
To address the first question of our study, we tested several properties of the Eastern lineage, to establish if the pattern of its mitochondrial variation is consistent with its derivation from a refugium located in the vicinity of the Carpathian Mountains. The distribution of the lineage clearly supports this hypothesis, as it is present in a wide area of central Europe, with the Carpathians in the centre of its range. Moreover, we observed a continuous spread of the lineage north of the Carpathians. Three diversity indices also provide evidence in support of the hypothesis (Figs 3 and 5). Relatively high levels of genetic diversity are expected in refugial areas for temperate species, due to the survival of large, demographically stable populations over a long period of past climatic oscillations [77,78]. Here we show, using both mtDNA (nucleotide diversity) and microsatellite markers (allelic richness and expected heterozygosity) that there are hotspots of genetic diversity for the common vole in the Carpathians and their close proximity. The interpolations and plots showed a distinct reduction of genetic diversity with distance from the Carpathian Basin and are also consistent with a northward route of recolonization from there (Figs 3 and 5). Another, previously proposed, refugial area for this lineage is the Balkan Peninsula [24], but this seems unlikely as this area is occupied by another mtDNA lineage [29]. Although there are fossil records of the common vole from the vicinity Phylogeography of a Common Vole Lineage of the Carpathians at the time of the LGM [32,33], ancient bone fragments without associated DNA cannot give an indication of which lineage was present at that time [79][80][81]. A very similar pattern of contemporary DNA variation to that of the common vole is observed in the bank vole where the 'Carpathian' lineage is more widespread in central Europe than the 'Balkan' lineage, which is present only in close proximity to its refugial area [9]. In the case of the bank vole, gene flow and demographic expansion were also observed from the Carpathian region to the south, not only northwards. This pattern was inferred using an isolation by migration (IM) model, which employs a Bayesian coalescent method and relies on the assumption that there is continuing gene flow between the two descendant populations which were derived from the split of the ancestral population. The estimated rate of gene flow into the 'Carpathian' lineage was equal to zero, while the gene flow outward from this lineage was significant [9].
It may be also hypothesized that the Eastern lineage of the common vole derived from refugia located somewhere close to the Ural Mountains or the Russian Plains and spread from such distant areas to central Europe after the LGM. For instance, Deffontaine et al. [82] suggested that the Eastern clade of the bank vole most likely derived from a refugium in the vicinity of the Urals and expanded its range to central Europe. However, the range of the Eastern lineage of the common vole does not include areas to the east of the Vladimir region, where the obscurus chromosome race occurs. Therefore, this hypothesis does not seem to be supported.
The second purpose of our study was to establish the timing of expansion of the Eastern lineage and to determine how it relates to events which occurred during the last glaciation. Coalescent genealogy sampling in BEAST, using the molecular clock rate that was established from common vole ancient DNA [28], inferred a median tMRCA for the Eastern lineage of 22.4 kya (95% HPD limits of 13.5 and 29.8 kya). This indicates that the Eastern lineage originated from a genetic bottleneck around the time of the LGM, when the species was known to be present in the Carpathians [32,33]. Based on the mismatch Tau value, the median time for the onset of expansion of the Eastern lineage was 8.3 kya (95% CI from 4.8 to 13.9 kya). This date is broadly consistent with the time of demographic expansion that was estimated from coalescent genealogy sampling, in the present and earlier studies, at around 10 kya (S2 Fig,  also see Fig 4 in [30]). These various estimates suggest that the Eastern lineage of the common vole was present in the Carpathian Basin during both the LGM and YD periods and re-colonized central and eastern Europe from this refugium at the beginning of the Holocene, after the last glacial retreat.
The final purpose of our study was to determine if the microsatellite differentiation that we found in our previous study of genetic variation in common voles from Poland [34] would be confirmed in a broader area of central Europe and, if so, what information it reveals about the population history of the species in the region. The previous study [34] described an easternwestern subdivision in Poland for both M. arvalis and M. agrestis. In this study, STRUCTURE analysis with samples of M. arvalis suggested that this pattern is also evident on a wider scale in central Europe. The eastern group occurs across eastern Poland, eastern Hungary, Russia and Serbia and the western group is present in western Poland, western Hungary and the Czech Republic (Fig 4 and S4 Fig). Again, as in the previous study, there were mixed (presumably hybrid) populations with intermediate levels of assignation to one of the two genetic groups. These populations are located between the eastern and western groups, in a broad contact zone (see Fig 4 and S4 Fig). However, there is a need to be cautious with the STRUCTURE results, because three possible values of K or genetic clusters (K = 2, 3 or 9) were inferred. The choice of K = 2 is based on the a priori data from [34], and cannot therefore be viewed as definitive. Results from a second spatial genetics program, sPCA, showed very weak differentiation between all populations across central Europe. The sPCA revealed no genetic structure and no division into genetic groups was indicated (see S6 Fig). These results were supported by low overall differentiation (F ST = 0.04) and significant differentiation between only some pairs of populations. There was a similar lack of confirmation of the STRUCTURE results in the previous study within Poland [34], where STRUCTURE revealed two genetic groups in Poland but Discriminant Analysis of Principal Components (DAPC) showed weak differentiation across the country.
This study is the first to examine genetic variation, using microsatellite data, in common vole populations from a wide area of central Europe. Previous studies were confined to western Europe [28,31]. The samples used in these microsatellite analyses were collected from an area where the Eastern mtDNA lineage predominates, and only five populations from north-western Poland were assigned to the Central mtDNA lineage. As in our earlier study [34], we did not find congruence between the mtDNA and microsatellite results, as the few Central lineage populations were not distinctive from the Eastern lineage. This differs from what has been found in Switzerland, where microsatellites also show differentiation in the contact zone of the Central and Western mtDNA lineages of common vole [31,73,74]. The number of loci used in the different studies is similar: eight loci in this study, 11 loci in [34] and 12 loci in the study in Switzerland [31]. Given that relatively few samples that were not from the Eastern lineage were examined in our study, it is unwise to use our data to make strong statements about differentiation between lineages. However, we are in a position to discuss the degree of substructure within the Eastern lineage.
We assume that the weak eastern-western effect may reflect historical or ecological factors that influenced the common vole in this region. Division into two genetic groups may have resulted from environmental factors during the late Pleistocene and early Holocene. During the YD river corridors were very broad and covered by permafrost patches [83] and probably the Vistula river was a geographic barrier hindering gene flow between vole populations during that time. Around 10 kya, during the period of rapid warming at the beginning of the Holocene, deciduous forest expanded and the Polish territory became influenced by oceanic climate [84]. This resulted in flooding over a large area, which could have formed a barrier to re-colonization in the common vole. The middle Vistula river valley, built of sand dunes, was mostly reshaped by wind between 10.7 and 10.5 kya but its stabilization continued until 9.3 kya. In the case of the Hungarian Plain this process took even longer [84]. The K = 2 pattern on the map (Fig 4) seems to correspond quite well to the Vistula and Danube river valleys. The pattern observed for K = 3 has similarities to that obtained for K = 2. In that case the western group still has the same distribution while the eastern group shows a degree of subdivision. According to the conclusions presented above and our previous study [34], we assume that the pattern obtained for K = 2 may best explain the observed diversity in the common vole populations in central Europe. In either case, the pattern has become largely obscured by later gene flow, which explains the lack of structure that was recovered with the sPCA.
Extensive sampling allowed us to perform an in-depth analysis of genetic diversity throughout the Eastern lineage of the common vole. Results of our study suggest that this lineage expanded from the Carpathian refugium and support an assumption that the Carpathians play an important role as a biodiversity hotspot (see in [85] for a recent review). Our study also supports a hypothesis that a part of central and eastern Europe (present-day Poland) does indeed form a phylogeographic suture zone [10]. Other species display a contact zone in Poland, for example, the bank vole, the common shrew (Sorex araneus) and the weasel [10,16,86]. Phylogeographic investigations in the central and eastern parts of Europe are important to examine colonization processes, including the meeting of lineages and the relative impact of southern and northern refugia on the postglacial colonization history of Europe.  Table. List of those collection localities from central and eastern Europe mapped in Fig  1B which provided previously published data for this study but which were not used to generate new cytb or microsatellite data. The full list of previously published cytb sequences that were used in this study are available in S1, S2 and S4 Tables in [30]. The numbers in square brackets are consistent with the References section. The reference not used in the References section is marked by ' Ã ' and its full description is given below the table. Murariu, S. Pavlova, L. Rychlik, S. Wąsik, P. Veligurov, V. Vohralik, and K. Zub for their help in collecting samples. We thank M. Górny for much helpful advice on GIS software, I. Ruczyńska for help in the lab and B. Johnson for help with use of Clumpp and Distruct. We also thank the three anonymous reviewers for their helpful comments, which have improved this manuscript.