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

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 obtained in cooperation with the Natural History Museum of Crete (Petros Lymberakis), the Slovenian Museum of Natural History (Boris Kryštufek), and the University of Agricultural Sciences and Veterinary Medicine Cluj Napoca (Attila Sandor). For the study, 314 samples were analysed from the following European countries: Czech Republic (49), Slovakia (26), Hungary (9), Romania (43), Slovenia (44), Italy (2), Croatia (20), Bosnia and Herzegovina (25), Montenegro (6), Serbia (10), Bulgaria (4), Macedonia (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.

Non-spatial analyses of genetic variability

The sizes of the microsatellite alleles were determined using GeneMarker V1.85 (http://www.softgenetics.com). Data binning was performed with the AutoBin software (http://www4.bordeauxaquitaine.inra.fr/biogeco/layout/set/print/Ressources/Logiciels/Autobin). The presence of the null alleles (false homozygotes), stutter bands and large allele dropouts was tested using Micro-Checker (Van Oosterhout et al., 2004). Bayesian statistics for recognising hybrid categories was run in NewHybrids 1.1 (Anderson & Thompson, 2002).

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^® (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²) were computed for each population in DnaSP v.5 (Rozas et al., 2003). Fu’s F_S(for big sample sizes) and R² (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 isolation-by-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.

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 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 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).

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 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 existenceof 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).