Glacial allopatry vs. postglacial parapatry and peripatry: the case of hedgehogs

Although hedgehogs are well-known examples of postglacial recolonisation, the specific processes that shape their population structures have not been examined by detailed sampling and fast-evolving genetic markers in combination with model based clustering methods. This study aims to analyse the impacts of isolation within glacial refugia and of postglacial expansion on the population structure of the Northern White-breasted hedgehog (Erinaceus roumanicus). It also discusses the role of the processes at edges of species distribution in its evolutionary history. The maternally inherited mitochondrial control region and the bi-parentally inherited nuclear microsatellites were used to examine samples within the Central Europe, Balkan Peninsula and adjacent islands. Bayesian coalescent inference and neutrality tests proposed a recent increase in the population size. The most pronounced pattern of population structure involved differentiation of the insular populations in the Mediterranean Sea and the population within the contact zone with E. europaeus in Central Europe. An interspecies hybrid was detected for the first time in Central Europe. A low genetic diversity was observed in Crete, while the highest genetic distances among individuals were found in Romania. The recent population in the post-refugial area related to the Balkan Peninsula shows a complex pattern with pronounced subpopulations located mainly in the Pannonian Basin and at the Adriatic and Pontic coasts. Detailed analyses indicate that parapatry and peripatry may not be the only factors that limit range expansion, but also strong microevolutionary forces that may change the genetic structure of the species. Here we present evidence showing that population differentiation may occur not only during the glacial restriction of the range into the refugia, but also during the interglacial range expansion. Population differentiation at the Balkan Peninsula and adjacent regions could be ascribed to diversification in steppe/forest biomes and complicated geomorphology, including pronounced geographic barriers as Carpathians.

(2009) revealed a positive correlation between body size and temperature and a negative correlation between body size and summer precipitation.
Association between particular species and their habitats give rise to the need of comparative phylogeographic studies. Dependency of Erinaceus hedgehogs on deciduous woodland (Rautio et al., 2014;Jensen, 2004) is in concordance with similar phylogeography observed in Oak and Silver Fir (Hewitt, 2000;Seddon et al., 2001) and implies that the extension of forests was followed by the hedgehogs. Some literature highlights an E. roumanicus preference for lowlands and open landscapes, particularly during the spatial expansion (Bolfíková & Hulva, 2012;Seddon et al., 2002). However, modern studies based on past environmental niche modelling are missing. Moreover, ecological considerations are based rather on knowledge of Erinaceus europaeus, while data on E. roumanicus are scarce. Different position of these two species ranges on often overlooked oceaniccontinental axis, represented by more humid and less seasonal climate on the west and the opposite in the east, should be also considered (Stewart et al., 2010). The southern parts of E. roumanicus range are located mainly in drier biomes in the Mediterranean Basin and the Pontic-Caspian steppe. The northern parts of the range are situated mainly in the temperate deciduous forests. Adaptations to the more dry open and seasonal habitats in E. roumanicus could not be excluded (Bolfíková & Hulva, 2012).
Do the classical models have a potential to further develop phylogeogeographic paradigm using up to date tools of the landscape genetics? Here, we provide new information on the hedgehog's evolutionary history using a combination of fast evolving markers, including nuclear microsatellites and mitochondrial control region. By analysing multiple E. roumanicus samples ranging from the post-refugial area in the Balkan Peninsula to the parapatric contact zone with E. europaeus in central Europe, we aim to: (1) assess the population structure of E. roumanicus within the studied area using the spatial and non-spatial Bayesian clustering methods and discuss observed patterns within landscape genetic framework in the light of geomorphology and habitat quality of the region; (2) reveal the history of island populations and suggest possible way of colonisation of these islands. He hypothesize that insular effects like the founder effect, genetic drift and other factors which decrease genetic variation highly affected these populations; (3) analyse the populations in the secondary contact zone from Central Europe. We hypothesize that respective populations were affected by species interactions, including introgression and reinforcement; (4) analyse the populations from the post-refugial areas at the Balkan peninsula and assess the hypothesis of occurrence of the glacial subrefugia.

Sampling and genotyping
Samples obtained from muscle tissue or digits or ears were collected from road-killed animals or from museum collections. No special permit is needed from the Institutional Review Board. The tissues were fixed in 96% ethanol and stored at −20 • C. The samples were  (17) and Greece (59), which included several islands (the Ionian islands, Euboea, Cyclades, and Crete). Each sample's origin, GPS coordinates, and date of collection along with the name of the person who collected it are listed in Table S1. A map (Fig. 1), indicating the geographic origins of the samples, was created using Geographic Information System (ESRI, 2014). DNA was extracted using a commercial DNA Blood and Tissue Kit (Qiagen, Hilden, Germany). A combination of nine microsatellite markers and mtDNA control region sequences was used. Laboratory procedures were performed according to the protocols by Bolfíková & Hulva (2012), details may be found in Table S3. Each dataset was analysed separately as specified below.
Evaluation of the genetic structuring of the microsatellite dataset was performed with a Bayesian Clustering Analysis in Structure 2.2 (Pritchard, Stephens & Donnely, 2000). Structure determines the most likely assignment of individuals to clusters (K ) based on analysis of likelihood. Structure is powerful to detect clinal variability and high geographic admixture in the dataset (Chen et al., 2007). The parameters of the final run were as follows: 1,000,000 MCMC steps, with burn-in-period 100,000 and 5 iterations for each K . An Admixture Model of Correlated Frequencies was used. The number of tested population was within the K = 1-14 interval. The results of the analyses were combined and evaluated by Structure Harvester (Earl & VonHoldt, 2012) using the Evanno method (Evanno, Regnaut & Goudet, 2005). Graphic visualisation was provided using Distruct 1.1 (Rosenberg, 2004). The inbreeding coefficient (F IS ) for each population was tested by GenePop (Rousset, 2008). Departures from the Hardy-Weinberg equilibrium (HWE), number of alleles (N a ) and allelic richness (A R ) were estimated in FSTAT (Goudet, 1995). Allelic richness was used to estimate an allele numbers, across all loci, corrected for equal sample size in all analysed populations. The population with the smallest number of individuals was N = 24 from Crete. The expected (H E ) and observed (H O ) heterozygosities were computed using Genetix (Belkhir et al., 1999).
Sequences of the mtDNA partial control region were edited using Geneious (Kearse et al., 2012) and aligned using MAFFT (Katoh & Standley, 2013). Particular haplotypes were identified in DnaSP v.5 (Rozas et al., 2003). The haplotype data were submitted to GenBank R (Accession numbers are from KY489901 to KY489953). The relationships among the haplotypes were visualised by the median-joining network (Bandelt, Forster & Röhl, 1999) using Network 4.5 (http://www.fluxus-engineering.com). The demographic parameters were estimated by summary statistics of the genetic variability for each population (recognised by the Bayesian clustering analysis of the haplotype data using the Geneland software; details in spatial analyses below). Tests for haplotype diversity (h; represents the probability that two randomly chosen haplotypes are different), nucleotide diversity (π; the average number of nucleotide substitutions per site between two sequences) and neutrality tests (Tajima's D, Fu's F S and R 2 ) were computed for each population in DnaSP v.5 (Rozas et al., 2003). Fu's F S (for big sample sizes) and R 2 (for small sample sizes) are considered as the most powerful tests (Ramos-Onsins & Rozas, 2002;Ramírez-Soriano et al., 2008). Significantly negative results indicate presence of selective sweep or strong population growth signal in the data. The course of the past population size was estimated with coalescent-based Bayesian skyline plots (Drummond et al., 2005). We inferred the model of sequence evolution in jModelTest (Posada & Crandall, 1998) using Bayesian information criterion (BIC). The Bayesian skyline plot analysis was conducted using the BEAST 1.4.8 program (Drummond & Rambaut, 2007) with a strict molecular clock and a GTR+I+G substitution model. The MCMC procedure was run three times for each species with 30,000,000 iterations, and the genealogy and parameters of the model were stored every 1,000 iterations. The results were combined in LogCombiner, and burn-in was set to 10,000,000 iterations in each run. The Bayesian skyline plot and the convergence of chains were assessed by Tracer v. 1.4 (Rambaut et al., 2014).

Spatial analyses of genetic variability
To analyse the spatial genetic architecture of E. roumanicus in the area of interest, we used two spatial model based clustering approaches, including software Geneland 4.0.4 (Guillot, Mortier & Estoup, 2005) and TESS (Chen et al., 2007;Durand et al., 2009). Geneland encompasses Bayesian clustering of individual genotypes based on separation of random mating subpopulations and produces Dirichlet cells centred on subpopulation territories. Geneland can analyse both haplotype data and codominant data in independent runs. The parameters for the Geneland analysis were set for each marker type as follows: 500,000 MCMC iterations, thinning = 100, five independent runs and K = 1-15. The posterior probability of the subpopulation membership for the most supported run was computed for each pixel of the spatial domain, using a burn-in of 50 iterations. An uncorrelated Frequency Model was chosen. TESS is based on minimizing of Wahlund effect and produces Dirichlet cells based on individuals and is using microsatellite data. We chose BYM model due to lower Deviance Information Criterion (DIC) and higher stability of algorithm. The first run of program was set to Kmax = 2-10 populations, then we plotted Kmax versus DIC values and selected the optimal Kmax (see details in Fig. S1). The final run of TESS was set to Kmax = 5-6, with 50,000 sweeps, a burn-in 30,000, the interaction parameter φ = 0.6 and 100 iteration for each K. The lowest DIC values were performed for Kmax = 6 and we chose as representative the run within 30% of runs with lowest DIC values, displaying five colours.
Combination of these approaches may be useful in exploring spatial population structure (Balkenhol et al., 2015), as Geneland performs well in detecting recent linear barriers to gene flow (Blair et al., 2012), while TESS addresses situations with barriers of complex shapes and permeability for migrants, as well as secondary fusions of subpopulations (Chen et al., 2007).
The Mantel test (Mantel, 1967) was used to characterise the presence of isolationby-distance to evaluate the correlation between congruent similarity/dissimilarity matrices, which addressed the relationship between the geographic and genetic distances of observations across a landscape provided by the Alleles in Space (AIS) software (Miller, 2005). The number of iterations used was 1,000. We used the AIS software for the Genetic Landscape Shape Analysis to compute the genetic distances, which were assigned to the connector centres of each locality using Delauney's triangulation. The interpolation procedure, which determined the presumed genetic distances in a lattice frame containing 15-km units, covered the area of interest, and a graph of the interpolated distances was created. Peaks in this graph represent an area of high genetic distances and should be correlated with areas of restricted gene flow. Both marker types were independently analysed using AIS software.

RESULTS
We did not detect any genotyping errors that resulted from a large allele dropout or stutter bands in 260 samples that were successfully genotyped at nuclear markers. All genotyped individuals are indicated in Table S1. Micro-Checker detected null alleles in several microsatellite loci, and the estimated frequencies differed between loci and populations. The estimated frequencies were lower than 0.2432 for all populations examined. In the Cretan population, two loci, EEU12H (0.2432) and EEU3 (0.2099), showed high null allele frequencies using the methods described by Van Oosterhout et al. (2004). Primers were designed for E. europaeus (Becher & Griffiths, 1997;Henderson et al., 2000) and tested by Bolfíková & Hulva (2012). All genotypes were originally obtained from the tissues, which are less susceptible to the presence of null alleles than non-invasive samples. Further analyses proved that the Cretan population was influenced by the founder effect, genetic drift and inbreeding. Thus, the null allele estimation was probably influenced by the high homozygote frequency in the population.
Sequences produced by Sanger sequencing were of good quality and some individuals with ambiguous sites were re-sequenced or skipped from the analyses. A total of 296 partial control region sequences (404 bp) were obtained and represented 53 haplotypes among the analysed samples (see details in Table S1).

Non-spatial analyses of genetic variability
The nine selected microsatellite loci were polymorphic. The total number of alleles in a single locus varied from 5 to 23 for the entire dataset. One sample from the Slovakia (Trencianske Teplice) was omitted from further analysis because it was suspected of having a hybrid origin. A follow-up analysis with the NewHybrids program (Anderson & Thompson, 2002) using the original dataset from Bolfíková & Hulva (2012) confirmed that the sample was a backcrossed hybrid of E. europaeus and E. roumanicus; the mitochondrial DNA originated from E. roumanicus.
From the total number of tested clusters (K = 1-14), division into the three populations (K = 3) obtained the highest support according to Evanno, Regnaut & Goudet (2005). The populations were differentiated as follows: (i) samples from Crete and from the Hellenic peninsula in southern Greece; (ii) samples from the Balkan area (Serbia, Bulgaria, Macedonia, Greece, Montenegro, Croatia and Bosnia and Herzegovina) plus Slovenia and Romania; and (iii) samples from Central Europe (Czech Republic, Slovak Republic and Hungary) (Fig. 2A). The second highest K value was observed for a division into the seven populations (K = 7). An increasing number of clusters resulted in further differentiation of the Balkan population by separation of the population from Pannonia (at K = 4), differentiation of the Slovenian population (K = 5) and separation of populations from the west and east of the Balkans (from K = 6). However, there is certain degree of admixture among particular clusters. Two different landscape genetic algorithms supported similar spatial pattern with three (Geneland, Fig. 5D) and five (TESS, Fig. 2B) clusters.  The aforementioned differentiation was also tested without the Czech and Slovakian samples, whose genetic architecture could have been influenced by phenomena associated with the recently developed contact zone. This analysis of Structure supported the previous subdivision of our dataset.
The genetic variability parameters are given for the seven populations (Table 1) and three populations (Table S2). H O was significantly lower than H E and GenePop analyses showed a significant lack of heterozygotes (p < 0.01) for each population. Crete had the lowest A R and H E . The highest F IS was detected in Romanian population.
The median-joining network didn't demonstrate reciprocally monophyletic lineages in the studied area (Fig. 3). Haplotypes are shared across populations, but some trends are Figure 3 Median-joining network of the mitochondrial control region haplotypes for 296 Erinaceus roumanicus samples. Haplotypes are represented by circles, which are proportional to the haplotype frequency. The hypothesized haplotypes are represented by red dots. Colour codes are used according to the population assignments of the samples using mtDNA data from the Geneland analysis (see Fig. 5C or Table S1). Lines between nodes represent a nucleotide changes. visible. One unique haplotype shared with the majority of individuals from Crete and a few other haplotypes carried by one or two individuals are representing the population of South Balkan (in blue on Fig. 3). Similar pattern is representing Central Europe (in orange on Fig. 3). The North Balkan population (in pink on Fig. 3) has two locally common haplotypes, but individuals from Central Europe and the Dinaric area are also share them. The Carpathian population (in yellow on Fig. 3) had several haplotypes in common with other populations but also several unique haplotypes. Samples from Dinaric area (in green on Fig. 3) were wide-spread in the network and produced star-like pattern around haplotype that is common in Slovenia, Hungary and Croatia. Black sea coast population (in violet on Fig. 3) has several unique haplotypes which are widespread within the network.
The genetic variability parameters for mitochondrial data were estimated for six populations according to the results of analysis in Geneland (Table 2). All tests of neutrality Table 2 The summary statistics of the partial control region of mitochondrial DNA in six populations of Erinaceus roumanicus, which were recognised by Bayesian analysis using Geneland. Colours in bracelets correspond to colours used in Fig. 5C. Number of individuals (N ), number of haplotypes (N h ), haplotype diversity (H d ), number of polymorphic sites (N p ), nucleotide polymorphism (π ), Tajima's D (D), Fu's Fs (F S ) and R 2 are given for each population and for whole dataset.  yielded negative values indicating recent population growth, only Central European (containing Czech Republic, East Slovakia and North Hungary) population neutrality tests were not significant ( Table 2). The Bayesian skyline plot indicated a long period with a stable population size and its recent increase (Fig. 4).

Spatial analyses of genetic variability
The Mantel test indicated a significant isolation by distance (r = 0.303; p < 0.001) within our dataset. Geneland analysis of microsatellite markers determined that the division of our dataset into three populations (Central Europe, mainland area, and south Greece plus islands) was the most likely to occur (Fig. 5D). Analysis implemented by TESS revealed the most fitting Kmax = 6, which spatially represents five clusters corresponding to population from Crete, Czech Republic, Pannonian area, Eastern and Western part of Balkan Peninsula (Fig. 2B). This pattern is in concordance with results of Structure analysis for K = 5 ( Fig. 2A). Mitochondrial analysis divided our dataset into the following six clusters: (i) South Balkan (southern and central Greece, Peloponnese, Crete and the adjacent islands); (ii) Dinaric area (north Greece, Macedonia, south Bosnia and Hercegovina, south Croatia, Montenegro and south Serbia); (iii) Black sea coast (Thrace and most of Bulgaria and southeast Romania); (iv) Carpathian area (northwest Romania, east Slovak Republic and east Hungary); (v) North Balkan (north Serbia, north Bosnia and Herzegovina, north Croatia, Slovenia, southwest Hungary and north Italy); and (vi) Central Europe (northwest Hungary, west Slovak Republic and Czech Republic) ( Fig. 5C; Table S1). Genetic Landscape Shape Analysis using AIS software detected very low genetic distances in Crete and the Czech Republic based on the mitochondrial data and nuclear data in particular. However, high genetic distances were observed with the microsatellites among samples in regions to the south and east of the Balkan Peninsula, with the highest peaks observed in Romania (Fig. 5B). The mitochondrial data showed flat genetic landscapes (Fig. 5A).

DISCUSSION
Based on the results obtained using both non-spatial and spatial individual based clustering approaches, the population from the southernmost region of the range (including southern Greece and the adjacent islands) and the population from the contact zone with E. europaeus were separated from the rest of the range. This pattern indicates that processes acting at peripatric and parapatric range edges shaped the genetic architecture of the species substantially. Further differentiation in clustering hierarchy provides evidence for occurrence of subpopulations within southeast part of the range. This result could be ascribed to complicated geomorphology and land cover of the region and may imply also subdivision of refugial populations in the past.

Insular populations in the Mediterranean Sea
All analyses indicate strong differentiation of the individuals from Crete, Cyclades and from Peloponnese, Euboea and southern and central Greece. As the largest island of Greece and the fifth largest Mediterranean island, Crete is typified by impoverished native mammalian fauna and the high impact of recently introduced species. However, Crete is also known as a generator of novelty in the numerous lineages that were able to colonise the island. From a macroevolutionary viewpoint, Crete is the place where pygmy forms of deer, elephants and hippopotamuses occurred (Van der Geer et al., 2010). From the microevolutionary viewpoint, the cryptic bat species occur there (Hulva et al., 2007;Hulva et al., 2010). Peloponnesus is also typified by its cryptic variability in vertebrates (Gvoždík Figure 5 Genetic landscape shape interpolation analysis. Genetic landscape shape interpolation analysis based on variation of (A) the control region of mtDNA using 296 Erinaceus roumanicus samples and (B) nine nuclear microsatellite loci using 260 E. roumanicus samples. The X and Y axes reflect the geographic coordinates within the study area, while the Z axis represents the average genetic distances between the analysed samples. The Bayesian landscape genetic analysis using Geneland illustrates the spatial distribution of the populations using (C) the control region of mtDNA and (D) nine nuclear microsatellite loci.  et al., 2010). The isolation of this region was recently enhanced by the 1893 construction of the Corinth Canal, which effectively separated Peloponnesus from the mainland. However, our analyses do not confirm a reciprocal monophyly for the E. r. nesiotes subspecies, which supports the findings by Schaschl, Lymberakis & Suchentruck (2002). Moreover, we ascertain that in addition to Crete, the range of this population includes other islands and reaches continental Greece. The assessment of the geographic origins and dispersal scenarios of this population will require more detailed analysis. However, the presumably absent land bridge between Crete and the mainland during the Pleistocene period (Schule, 1993), the absence of hedgehogs in Crete's fossil record (Lax & Strasser, 1992) and the frequent introductions of hedgehogs onto the islands by humans (Bolfíková et al., 2013) indicate that the role of the anthropogenic factor in the phylogeography of the southern cluster is highly probable.
The first archaeological evidence of hedgehogs in Crete was collected from an Early Minoan tomb (Jarman, 1996), but the author concluded that the hedgehog might have been an intrusion (e.g., entered the tomb much later). Using allozymes and partial cytochrome b sequences, Schaschl, Lymberakis & Suchentruck (2002) proposed that hedgehogs were introduced to Crete by archaic settlers from the mainland. This population was previously described as a subspecies, E.r. nesiotes, by Dorotha Bate in 1906 based on phenotypic data. These facts imply the occurrence of insular syndrome in the Cretan population. The markedly decreased genetic variation values in the insular population imply that factors such as the founder effect and genetic drift played substantial roles. The agreement between the genetic variation patterns and the phenotypic data indicates that insular syndrome, including dwarfism in hedgehogs from the Greek islands, has occurred (Kryštufek et al., 2009;Giagia-Athanasopoulou & Markakis, 1996). Despite the uncertainty surrounding the origins of differentiation for the southern cluster, the aforementioned results indicate a pronounced role for peripatric processes in the microevolution of E. roumanicus. Elevated evolutionary rates in islands are obvious also in phenotypic evolution. The trends in body size changes of insular hedgehogs are negatively correlated with the island's distance from the mainland, however, the insular response was not uniform (Kryštufek et al., 2009).

Secondary contact zone with E. europeaus in Central Europe
The population from Central Europe is genetically different from the rest of the mainland, as suggested by all Bayesian clustering analyses. This cluster is also characterised by its lower genetic variation, even though, in some parameters still comparable to post-refugial populations. This may be the outcome of spatial expansion, which is usually followed by a decrease in genetic variability due to the bottleneck effect. More complex phenomena, such as allele surfing (increase of an allele frequency in the wave front of an expanding population), might alter the genetic architecture of the population's expanding edge and change the genetic variation as well (Klopfstein, Currat & Excoffier, 2006;Excoffier & Ray, 2008;Excoffier, Foll & Petit, 2009). Demographic analyses did not confirmed statistically significant population size increase. Thus, we suspect that the aforementioned phenomena played important roles in shaping the population genetics of Central European hedgehogs. Other factors that caused population differentiation in this case may be associated with parapatry and its consequent species interactions. Ecological interactions and their consequent selections may have influenced the diversity of the studied microsatellite loci by genetic hitchhiking, while hybridisation and introgression in the secondary contact zones may have further shaped the allopatrically evolved lineages and potentially cover the signal of recent population expansion. These potential factors will require further study and genomic tools.
We are the first to report a single hybrid sample from the range margins of E. europaeus in Central Europe (Slovakia-Trenčianske Teplice) using traditional genetic markers.
Interspecies interactions in the sympatry zone may have previously played roles in forming reproductive barriers. We hypothesised that hybridisation occurred when the two species came into contact for the first time, when the isolating mechanisms were not yet formed. Afterwards, during the reinforcement phase, species developed mechanisms for reproductive isolation and the hybrid event frequency decreased for follow-up generations. The history of possible hybridisation events in the Central European population should be tested using approaches that include whole genome data because the introgression modes may be complex. Future investigations should address the contact zone at the Slovenian/Italian border (located to the south of the Alps) and the effect of altitudinal heterogeneity on both species. This region involves the contact zone for both species and the transition zone between the Slovenian mountains and the Padan Plain.
The lower genetic diversity associated with the range margins and the novel selective pressures applied by interactions with related species may be why the Central European population was highly divergent from the other populations in the dataset.

Post-refugial areas at the Balkan Peninsula and adjacent regions
The separation of the Pannonian population represents the most pronounced pattern after the split of parapatric and peripatric populations. The Pannonian basin is a confined area from the view of geographical isolation by the Carpathian arch and vegetation cover typified by grass steppe, i.e., it possess habitat similarity with the part of species range in Pontic-Caspian steppe. Consequently, it represents a pronounced biogeographical region of Europe and a presumed continental interglacial refugium of continental and steppe faunal elements (Stewart et al., 2010;Ricanova et al., 2011). In E. roumanicus, this pattern might mirror geographical isolation, ecological differentiation, and also the quality of ancestral habitat and larger extent of steppe biome during glacial periods.
Further differentiation occurs between populations from the west of the peninsula at the Adriatic coast and the east at Black Sea coast. This differentiation might be ascribed to isolation by distance and barrier role of Carpathian Mountains. For example, based on the results provided by AIS, the largest genetic distances between the samples in our dataset are located in Romania. This pattern implies a restricted level of gene flow within the population, which might be ascribed to the high altitudinal heterogeneity of the Carpathian Mountains. However, this pattern may also mirror occurrence of forest refugia at Adriatic and Pontic regions. The Black Sea shore likely served as a glacial refugium for multiple animal species, such as Bombina bombina (Fijarczyk et al., 2011) and Emys orbicularis (Joger et al., 2007). The presence of trees during the LGM has also been proven (Tzedakis, 2004). Palynological analyses have confirmed occurrence of isolated forests within the refugia in response to climatic and spatial variabilities (Huntley, 1999), for example in the Pindos Mountains in northwestern Greece, where levels of precipitation were higher (Tzedakis, 2004). In the case of the European oak Quercus, two separate Balkan subrefugia (Greece and the west coast of the Black Sea) have been located (Brewer et al., 2002). King & Ferris (1998) discovered unique haplotypes of the European alder (Alnus glutinosa) in the regions around Bulgaria and Greece. Based on these examples, Tzedakis (2004) presumed the existence of temperate forests during the last glacial maximum in the mountainous and coastal areas, whereas two previously identified subrefugia occurred in western Greece and in eastern Bulgaria near the Black Sea.
The population from Slovenia is separated from the Balkan genotypes in Structure analyses. This pattern suggests a gene flow limitation due to the geographical isolation, caused by the Alps, the Dinarids and/or specific climatic conditions. For example, a fully grown forest was present during the LGM in Slovenia (Tzedakis, 2004). Another hypothesis suggests a possible species interaction with E. europaeus and, consequently, a local parapatic evolutionary process (see above).

CONCLUSIONS
The most pronounced pattern within our dataset was represented by the differentiation of populations at the edges of the recent range. This study indicates that peripatry and parapatry might be not only limiting factors to range expansion, but may also provide strong microevolutionary forces that shape the genetic structure of the species. Traditionally, population differentiation in temperate species has been ascribed to allopatric processes during glacial periods. Here, we provide an alternative example, showing that population differentiation may occur not only during the glacial restriction of the range into refugia, but also during the phase of interglacial range expansion. During the interglacial periods, peripatric modes of evolution for many taxa might be more frequent due to the rise in sea level, which cause increased shelf island formation. Parapatric processes are accelerated by the range expansions of lineages isolated in the refugia and by the formation of potential suture zones.
Recent population of E. roumanicus in post-refugial area at Balkan Peninsula and adjacent regions shows differentiation with particular subpopulation located in areas of interglacial steppic refugia and glacial forest refugia. A detailed habitat model and past niche modelling will be necessary to infer, which type of habitat was crucial for survival of E. roumanicus during glacial periods. However, certain degree of admixture among particular clusters is obvious, mirroring complex geomorphology and shape of geographical barriers, as well as altitudinal heterogenetity within Balkan Peninsula. The situation might be further complicated by the fact, that a mosaic of particular biomes often existed in refugia (Huntley & Birks, 1983), increasing the proportion of ecotones, and the biotas were often ''disharmonious'', i.e., particular ecosystems often contain species that inhabit different habitats today (Tyrberg, 1991).