Skip to main content

Contrasting evolutionary histories of the legless lizards slow worms (Anguis) shaped by the topography of the Balkan Peninsula

Abstract

Background

Genetic architecture of a species is a result of historical changes in population size and extent of distribution related to climatic and environmental factors and contemporary processes of dispersal and gene flow. Population-size and range contractions, expansions and shifts have a substantial effect on genetic diversity and intraspecific divergence, which is further shaped by gene-flow limiting barriers. The Balkans, as one of the most important sources of European biodiversity, is a region where many temperate species persisted during the Pleistocene glaciations and where high topographic heterogeneity offers suitable conditions for local adaptations of populations. In this study, we investigated the phylogeographical patterns and demographic histories of four species of semifossorial slow-worm lizards (genus Anguis) present in the Balkan Peninsula, and tested the relationship between genetic diversity and topographic heterogeneity of the inhabited ranges.

Results

We inferred phylogenetic relationships, compared genetic structure and historical demography of slow worms using nucleotide sequence variation of mitochondrial DNA. Four Anguis species with mostly parapatric distributions occur in the Balkan Peninsula. They show different levels of genetic diversity. A signature of population growth was detected in all four species but with various courses in particular populations. We found a strong correlation between genetic diversity of slow-worm populations and topographic ruggedness of the ranges (mountain systems) they inhabit. Areas with more rugged terrain harbour higher genetic diversity.

Conclusions

Phylogeographical pattern of the genus Anguis in the Balkans is concordant with the refugia-within-refugia model previously proposed for both several other taxa in the region and other main European Peninsulas. While slow-worm populations from the southern refugia mostly have restricted distributions and have not dispersed much from their refugial areas, populations from the extra-Mediterranean refugia in northern parts of the Balkans have colonized vast areas of eastern, central, and western Europe. Besides climatic historical events, the heterogeneous topography of the Balkans has also played an important role in shaping genetic diversity of slow worms.

Background

Diversity of European biota has been strongly influenced by global climatic and environmental changes in the Quaternary. Toward the end of the Pleistocene, repeated climatic oscillations led to extinctions of many phylogenetic lineages from vast northern areas during glacial periods followed by re-colonisations during interglacials [13]. Many plant and animal lineages survived cold and dry glacials in relatively stable and hospitable environments. In Europe these were located in three Mediterranean peninsulas: Iberian, Italian, and Balkan. This general biogeographical model has been expanded to a more complex view acknowledging long-term persistence of cold-tolerant species in central and northern Europe during glacials and survival in multiple refugia located within the Mediterranean peninsulas [4, 5]. Demographic stability of populations in southern refugia enabled them to diverge, which has resulted in high diversity in all three main refugial regions. In contrast, northern populations established during re-colonization are generally characterized by lower taxonomic and genetic diversity.

In comparison to the Iberian and Italian peninsulas, the Balkans has remained much less studied in terms of the biogeographical history of the species distributed there, although it is richer both in biodiversity and paleoendemics [68]. The Balkan Peninsula is not isolated by one extended mountain range such as the Pyrenees of the Iberian and the Alps of the Italian Peninsula, and so there are fewer dispersal barriers to the north. This allowed postglacial expansion of populations from the Balkan refugia to central and northern Europe [1, 5, 9]. On the other hand, the Balkan Peninsula is a region with high topographic and climatic heterogeneity, showing a strong contrast between the eastern/western and northern/southern parts. In the east and north, the surface is formed by plains or plateaus and the mountain slopes are generally gentle, while in the west and south the Dinarides and Hellenides rise steeply from the coastal strip [10]. Each of the Balkan mountain chains also has a different tectonic and sedimentary history, and while they all underwent complex folding and faulting in the process of the Alpine orogenesis, the intensity was different [11]. All this geographical variation offers suitable conditions for local adaptations of populations, which could promote divergence and subsequent diversification [12, 13].

We, and others [1417] have been studying the evolutionary history of legless lizards of the genus Anguis (family Anguidae) within its Western Palearctic range. This genus comprises five species, four of which occur in the Balkans [14, 15]. While Anguis cephallonica Werner, 1894 and A. graeca Bedriaga, 1881 are Balkan endemics with rather restricted distribution in the south of the peninsula, ranges of A. fragilis Linnaeus, 1758 and A. colchica (Nordmann, 1840) are at the continental scale and cover vast areas of Europe and western Asia [14, 1719]. Considering the semifossorial lifestyle and high site tenacity [20, 21], one might expect restricted occurrence of slow worms. However, distribution of slow worms in the Balkan Peninsula seems to be more or less continuous with gaps probably only in agricultural regions and extremely high altitudes [19, 22, 23]. Nevertheless, details of the species ranges within the Balkans, contact zones of multiple species, and detailed intraspecific genetic structure in respect to geography and ecology still remain widely unknown.

In this study we collected and analysed data originating from the Balkan slow-worm populations with the aim to i) provide a detailed picture of distribution; ii) infer historical relationships of populations and describe genetic diversity; iii) reconstruct biogeographical histories of the Balkan slow-worm populations during the Quaternary. Finally, we tested iv) whether the genetic diversity observed in the Balkan slow worms is driven by specifics of topography. Dispersal barriers would most likely coincide with the extensive and variously rugged mountain ranges of the Balkan Peninsula, thus we expected the slow-worm genetic diversity to be correlated with topographic variation of this region.

Methods

Sampling

Since the Balkan Peninsula represents an important evolutionary centre of the genus Anguis, we devoted this study to slow-worm populations from this region. Our sampling strategy focused on equally representing the whole Balkan region as well as all four Balkan species. Sampling effort also took into account that these species vary in distribution ranges and inter-specific genetic diversification and may have low population densities in some areas. Tissue samples were obtained mainly from road-killed individuals or alternatively from living animals as oral swabs, blood droplets, or miniature skin biopsies. This sampling procedure did not affect survival of the captured animals. No experimental research was carried out on these animals in this study. All samples were preserved in 96 % ethanol. A portion, 732 base pairs (bp), of the mitochondrial gene for NADH dehydrogenase subunit 2 (ND2) was targeted. Newly produced nucleotide sequences were supplemented to previously published sequences from the Balkans [1417] to complete a total of 231 specimens from 187 localities. Based on the mtDNA identity, we represented all four Balkan Anguis species, namely 110 A. fragilis, 56 A. colchica, 49 A. graeca, 16 A. cephallonica (Fig. 1; Additional file 1: Table S1). To put our Balkan data into a complex phylogeographical context, we compiled an additional dataset supplemented by all known haplotypes, including those from outside the Balkans, previously published by [14, 15]: A. fragilis (f1–f15), A. colchica (c1–c12), A. graeca (g1–g16), A. cephallonica (ce1, ce2), A. veronensis (v1–15); and [16]: A. fragilis (AF01–AF07), A. colchica (AC01, AC02). The resulting dataset contained 271 sequences, excluding outgroup. Following previous works of our team [14, 15], we used the sister genus Pseudopus as outgroup (P. apodus thracius from Albania, Pat1, GenBank No. FJ666589).

Fig. 1
figure 1

Maximum likelihood (ML) phylogenetic tree of Anguis species and their distributions in the Balkans based on a fragment of mtDNA (ND2). Anguis veronensis (in white) occurs outside the area studied here. Numbers at nodes show Bayesian posterior probabilities and ML bootstrap support values. Yellow lines denote contact zones between two species. Numbers correspond to the locality numbers as given in Additional file 1: Table S1

Laboratory procedures

Total genomic DNA was extracted using various commercial kits and following respective manufacturer protocols. We amplified > 1400 bp-long portion of mtDNA comprising the complete ND2 gene, five subsequent transfer RNA (tRNAs) genes and the light-strand replication origin using primers (L4437n, H5934) and protocol following [14]. We sequenced only the first half of the amplicon using the internal reverse primer AND2inR2 [14], which was also used in PCR amplifications in cases of samples with degraded DNA, using the same protocol. Alternatively, the internal reverse primer AND2inRc [14] was used in A. cephallonica for both PCR amplifications (in degraded DNA) and sequencing. The final stretch contained 732 bp-long fragment of ND2 after trimming the low quality ends. The sequencing was performed by Macrogen Inc. (Seoul, South Korea or Amsterdam, Netherlands; http://www.macrogen.com) and new sequences have been deposited in GenBank under accession numbers KX020147–KX020322 (Additional file 1: Table S1).

DNA sequence evaluation, phylogenetic analyses, and haplotype networks

The protein-coding ND2 fragments (732 bp) were aligned manually. No stop codons were detected when the sequences were translated using the vertebrate mitochondrial genetic code in the program DnaSP 5.10 [24]. The same program was used to calculate uncorrected p-distances among the main lineages or haplogroups within each taxon, and to estimate the number of haplotypes (h), haplotype diversity (Hd), number of segregating sites (S), nucleotide diversity (π), and Watterson’s theta (θ W) for each of these lineages or haplogroups.

For phylogenetic analyses we used the all-individuals dataset supplemented by distinct published haplotypes from outside the Balkans to obtain a complex picture of the phylogenetic relationships within the Balkan Peninsula and in the framework of the whole genus. The best-fit codon-partitioning schemes and the best-fit substitution models were selected using PartitionFinder v1.1.1. [25], according to the Bayesian information criterion (BIC), separately for each dataset and methodological approach (i.e. models available in the used software). Phylogenetic trees were inferred using the Bayesian approach (BA) and maximum likelihood (ML) with the software MrBayes 3.2 [26] and RAxML 8.0 [27], respectively. Each codon position treated separately was selected as the best-fit partitioning scheme for both BA and ML with the best-fit substitution models for the BA analysis as follows: HKY + G (1st codon position), HKY + I (2nd codon position), and HKY + G (3rd codon position); and for the ML analysis: GTR + G in each codon position. The ML clade support was assessed by 1,000 bootstrap pseudoreplicates. The MrBayes analysis was set as follows: two separate runs, with four chains for each run, 10 million generations with samples saved every 100th generation. The convergence of the two runs was confirmed by the convergence diagnostics (average standard deviation of split frequencies, potential scale reduction factor). First 20 % of trees were discarded as the burn-in after inspection for stationarity of log-likelihood scores of sampled trees in Tracer 1.6 [28] (all parameters had effective sample size > 200). Majority-rule consensus tree was drawn from the post-burn-in samples and posterior probabilities were calculated as the frequency of samples recovering any particular clade.

Haplotype-network approaches can be more effective for presentation of intraspecific evolution than the tree-based phylogenetic approaches [29]. Therefore, we also constructed haplotype networks for individual species (or main clades in A. colchica) using the 95 % limit of parsimony as implemented in TCS 1.21 [30]. To infer possible connections to a network when cases of highly divergent haplotypes were detected (two haplotypes in A. graeca, and one in A. cephallonica), we also applied a fixed connection limit at a higher number of steps allowing visualization of their likely connections to the networks constructed under the 95 % limit of parsimony.

Demographic analyses

The past population dynamics of the main population groups were inferred using the Bayesian coalescent-based approach of the Bayesian skyline plots (BSPs; [31]) as implemented in BEAST 2.1 [32]. This method computes the effective population size (N e ) through time directly from sampled sequences and does not require a specific a priori assumed demographic model. Main population groups correspond to monophyletic groups. In a single case of several basal haplogroups of A. fragilis, the population group was defined geographically (‘Slovenian’ populations). Preliminary analyses were run using both strict molecular clock and uncorrelated lognormal relaxed molecular clock. Since the parameter of the standard deviation of the uncorrelated lognormal relaxed clock was close to zero, the final analyses were run enforcing the strict molecular clock model. A uniform prior for the substitution rate with the initial value 0.0065 substitution/site/lineage/Myr (as suggested for the used mtDNA marker in anguid lizards; [33]) was set as no internal calibration point was available. Using PartitionFinder v1.1.1. [25], all codon positions treated together as one partition and the HKY substitution model were selected as the best-fit partitioning scheme and the best-fit model, respectively, for each population group. The final BSP analyses were run in duplicates to check for consistency between runs, each for at least 10 million generations (or more according to each dataset until the effective sample size [ESS] > 200 was achieved) and sampled every 1000 generations (or more, accordingly, to save 10,000 samples). Convergence, ESS, stationarity, and the appropriate number of generations to be discarded as burn-in (10 %) were assessed using Tracer 1.6 [28]. The resulting BSPs were also summarized in Tracer 1.6 with the maximum times as the median of the root height parameter.

In addition, the mismatch distributions (MD) were calculated as the distributions of the observed pairwise nucleotide differences and the expected values under a growing- or declining-population model using DnaSP 5.10 [24]. The occurrence of historical demographic changes was assessed by the neutrality-test statistics of the Fu’s F S [34] Tajima’s D [35], and Ramos-Onsins and Rozas’s R 2 [36] calculated in DnaSP 5.10 with the estimation of the statistical significance using 10,000 coalescent simulations.

Genetic diversity and topographic heterogeneity

Since a more complex topography is more likely to limit dispersal and gene flow, we hypothesized that regions with higher topographic heterogeneity (terrain ruggedness) will be inhabited by slow-worm lineages characterized by higher genetic diversity. To test for this relationship we performed regression analyses of nucleotide diversity (π) with the terrain ruggedness index (TRI). TRI is a measure of topographic heterogeneity calculated as a sum change in elevation between a grid cell and its eight neighbour cells in a grid network [37]. Cell TRI values are then averaged across specific areas such as mountains. Values of TRI were derived from digital elevation model based on the data from the NASA Shuttle Radar Topographic Mission (SRTM-3; available at http://srtm.usgs.gov) with a spatial resolution of approximately 3 arc-sec (~100 × 100 m) with a final resample to 30 arc-sec (~1 × 1 km) using GRASS GIS 7.1 [38]. The polygon network was created for the selected topographic units, mountain ranges, which respect distributions of evolutionary lineages or haplogroups (Apuseni Mts., Carpathians, Dinarides, Hellenides, Prealps, Peloponnese, Macedonian-Thracian Massif, Stara Planina Mts.; Additional file 2: Table S2, Additional file 3: Figure S1 and Additional file 4: Table S3). Since higher genetic diversity might be expected in geographically larger areas, we controlled for the effect of the topographic-unit size. In the multiple linear regressions, we regressed nucleotide diversity of individual phylogenetic lineages/haplogroups against ‘extreme’ values of TRI calculated for the topographic units (‘extreme’ values were taken from the highest 25 % of data: 3rd quartile (Q3), and median and modus of the values above Q3), and against topographic-unit sizes (in km2). The ‘extreme’ values of TRI were selected with the aim to preferentially study the influence of steeper terrain, presumably posing stronger barriers to gene flow and resulting thus in higher probability of lineage divergence. Due to the controversy about the biogeographical significance of the Apuseni Mts. as a separate unit within the Carpathians [39, 40], we performed two separate analyses with samples from the Apuseni Mts. included and excluded, respectively, within the group of the Carpathian samples. The GIS analyses were performed using ArcGIS 10.1 (ESRI) and the multiple linear regressions were carried out using STATISTICA version 12 [41].

Results

Phylogeny, species distributions and contact zones

The maximum likelihood and Bayesian phylogenetic analysis provided topologies concordant with previous studies [14, 15, 17]. The southernmost species, A. cephallonica, forms a clade with A. veronensis from the Italian Peninsula, while the other three species (A. fragilis, A. colchica, and A. graeca) form a separate clade, in which the Balkan endemic A. graeca is in a sister position to the eastern widespread species, A. colchica (Fig. 1).

Anguis fragilis is distributed in the northwestern and central parts of the Balkan Peninsula from the Julian Alps and the southeastern Prealps, along the Dinarides to the Macedonian-Thracian Massif, and only marginally in the northern Hellenides (Figs. 1 and 10). Anguis colchica is documented from the Carpathians, the Balkanides along the Stara Planina Mts. in southeastern Serbia and central Bulgaria, and from the Black Sea region (Strandzha Mts.). Anguis graeca is mostly confined to the Hellenides in the southern Balkans where it is distributed from the northern Peloponnese, along the Pindus Mts. and the Albanian Mts. to the southernmost Dinaric region (southern Montenegro) and western Macedonian-Thracian Massif (northeastern Rep. Macedonia). Anguis cephallonica was found in the Peloponnese and Kephallonia Island (not sampled in Zakynthos and Ithaki islands in this study where the species was documented previously [42, 43]).

Our detailed sampling also revealed several areas where haplotypes of different species could be found in distances from ca. 15 to 80 km, indicating the existence of contact zones (Fig. 1). One such contact zone between A. fragilis and A. colchica was detected in eastern and southeastern Serbia and central-western Bulgaria (sites 34–37, 67–68, 73–76 for A. fragilis; sites 109–111, 115–116, 117–119 for A. colchica). Three zones of contact were further detected between A. fragilis and A. graeca in southernmost Montenegro (sites 61 and 139; sympatric occurrence), northwestern Rep. Macedonia (sites 63 and 158), and in the tri-border area of Serbia, Bulgaria and Rep. Macedonia (sites 64, 70–72, 160). Sympatric occurrence of A. graeca and A. cephallonica was confirmed from northern Peloponnese (sites 165, 166, 172; and 174–176).

Genetic diversity and phylogeographical patterns

The dataset built up from the Balkan specimens contained 231 ingroup (Anguis) sequences, which yielded a total of 100 haplotypes. Nucleotide diversity was higher in the two Balkan endemics, A. graeca (π = 1.17 ± 0.11 %) and A. cephallonica (π = 0.81 ± 0.21 %), than in the Balkan populations of the two northerly distributed taxa, A. colchica (Incerta clade; π = 0.66 ± 0.05 %) and A. fragilis (π = 0.34 ± 0.04 %; Table 1).

Table 1 Summary of genetic polymorphism and results of neutrality tests for the Balkan populations of four species of the genus Anguis

Anguis fragilis shows relatively low genetic variation, with 34 haplotypes identified among 110 individuals (intraspecific p-distance ≤ 1.1 %; Additional file 5: Table S4; Fig. 2). The basal radiation was detected in the northwest of the Balkans, in the northern Dinarides and southeastern Prealps (sites 1–7). Haplotypes from this basal radiation do not form a monophylum and may be divided into three Slovenian haplogroups, which we name in accordance to the detected distributions as follows: North Adriatic (sites 1, 2, 7); Carniolan (sites 5, 6); and Alpine-Pannonian (sites 3, 4). In earlier studies, haplotypes belonging to the latter haplogroup were also found outside the Balkans, i.e. in northeastern Italy (haplotype f8 – [15]) and the Pannonian Plain (haplotypes AF04, AF05 – [16]). Another haplogroup from the basal radiation (haplotypes f14, f15 – [15]; and AF07 – [16]) and a single haplotype (f7 [14]) conform to populations from Western Europe (Spain, France). All other A. fragilis haplotypes cluster into one large unit that might be divided into two geographically separated haplogroups: the northern we hereafter name the Illyrian-Central European haplogroup (ICE), and the southern one (the South Balkan haplogroup, SB). Haplotypes from the ICE haplogroup were also detected in Central Europe and southern Great Britain (haplotypes f1–f3, f12, f13 – [14, 15]; AF01–AF03 – [16]). The ICE haplogroup is paraphyletic in respect to the SB haplogroup. Nevertheless, the SB haplogroup is geographically well defined, confined to the Macedonian-Thracian Massif and only slightly penetrating to the northern Hellenides (Figs. 2b and 10). The ICE haplogroup is distributed along the Dinarides and surrounding lowland areas.

Fig. 2
figure 2

Anguis fragilis, (a) maximum likelihood (ML) phylogeny, (b) geographical distributions, and (c) parsimony haplotype network of the main haplogroups in the Balkans. Numbers at nodes in the tree represent Bayesian posterior probabilities (pp) and ML bootstrap support values (pp below 0.50 and bootstrap values below 50 are not shown). Locality numbers (in parentheses) follow sample IDs and correspond to the numbers in Additional file 1: Table S1. White circles in the network represent extralimital haplotypes as detected in previous studies (Gvoždík et al. [14, 15]; Szabó & Vörös, [16]; Thanou et al. [17])

In Anguis colchica, a deep intraspecific divergence (4.3 % p-distance; Additional file 5: Table S4) was found separating two clades of a different geographical origin (Fig. 3a). One clade is widespread and corresponds to the subspecies A. colchica incerta [14], hereafter named the Incerta clade, while the second clade was detected in the Black Sea coastal area, therefore named the Pontic clade. Outside the Balkans, A. colchica forms two additional clades distributed in the Caucasus (A. c. colchica; the Colchica clade) and the southern Caspian region (A. c. orientalis; the Orientalis clade); see also [14]. The Pontic clade is currently only known from the Strandzha Mts. in southeastern Bulgaria. The mtDNA polymorphism of the Pontic clade is relatively high (8 haplotypes within 13 specimens) in respect to its restricted geographical range (Fig. 3b, c). The Incerta clade (16 haplotypes within 43 specimens) is widespread along the Carpathians and the Stara Planina Mts., with relatively high genetic variation and diversified into three main well-supported lineages: (i) Stara-Planina lineage in the region of the Stara Planina Mts. and the northern foothills, reaching the Serbian Carpathians (sites 111, 112); (ii) Banatian lineage detected in the Banat (southwestern Carpathians); and (iii) Carpathian lineage present in most of the Carpathians with a further sub-structure forming at least four haplogroups; Carpathian I–IV (Fig. 3a). The Carpathian I haplogroup seems to be confined to Transylvania (Apuseni Mts. and their vicinity; sites 94–99), while the other three are partially sympatric. The Carpathian lineage contains haplotypes that were also detected outside the Balkans in earlier studies (c1–c6, c12 – [14, 15]; AC01, AC02 – [16]; Fig. 3c).

Fig. 3
figure 3

Anguis colchica, (a) maximum likelihood (ML) phylogeny, (b) geographical distribution of the main haplogroups in the Balkans, and (c) parsimony haplotype networks of the two Balkan clades (Incerta and Pontic). See the legend to Fig. 2 for more details

Of the two Balkan endemics, A. graeca shows a higher nucleotide but comparable haplotype diversity (29 haplotypes detected among 49 individuals) than the less widespread A. cephallonica. The genetic structure of A. graeca is complex, characterized by many haplogroups but without deep divergences (Fig. 4). Only two detected haplotypes (KJ634800, KJ634801) are relatively divergent both from each other and from all other haplotypes. They both originate from the same location in the northern Peloponnese (site 172). Geographical distributions of most haplogroups are restricted to small areas, mainly in the mountains of Albania (Fig. 4b). Only three haplogroups have wider distribution: one in central and southern mainland Greece, northern Peloponnese and Euboea Island (graeca I); the second in western Greece, Corfu Island and southern Albania (graeca V); and the third one in Rep. Macedonia (graeca XI).

Fig. 4
figure 4

Anguis graeca, (a) maximum likelihood (ML) phylogeny, (b) geographical distribution, and (c) parsimony haplotype network of the main haplogroups. See the legend to Fig. 2 for more details

The Peloponnese endemic, A. cephallonica, has a similarly complex phylogeographical structure with 13 haplotypes detected among 16 specimens (Fig. 5). One haplotype (KJ634795) originating from the Mani Peninsula in the south forms a lineage (hereafter the Mani lineage) divergent from all the other haplotypes, which form a well-supported monophylum (hereafter the Widespread lineage; 2.4 % p-distance; Additional file 5: Table S4; Fig. 5a). The Widespread lineage displays an inner diversification with several haplogroups distributed around the Peloponnese and Kephallonia Island with east–west longitudinal structure and higher diversity in the central Peloponnese (Fig. 5b, c).

Fig. 5
figure 5

Anguis cephallonica, (a) ML phylogeny, (b) geographical distributions, and (c) parsimony haplotype network of the main haplogroups. See the legend to Fig. 2 for more details

Historical demography

The Bayesian skyline plots (BSPs; Figs. 6, 7 and 8) gave evidence of population growth in all tested groups, with the exception of the ‘Slovenian’ populations of A. fragilis (Fig. 6c) and the Carpathian populations of A. colchica, although a mild and relatively recent (during the last ca 80 Ky) population growth was detected in the Carpathian lineage (Fig. 7b). A sharp population growth was detected in the Stara-Planina lineage of A. colchica also since ca 80 Kya (Fig. 7c). Comparing the two main clades of A. colchica, population growth started earlier in the Pontic clade (ca 200 Kya; Fig. 7d) than in the Incerta clade (80 Kya; Fig. 7a). Considerable population growth was also detected during the last 150 Ky in the ICE + SB haplogroups of A. fragilis (Fig. 6a), or since ca 50 Kya when only the SB haplogroup was analysed (Fig. 6b). Anguis graeca was analysed as a single population due to its complex genetic variation with many haplogroups. The BSP showed a substantial population growth starting about 700 Kya, the population being stable during the Middle Pleistocene and slightly declining during the last ca 80 Kya (Fig. 8a). In the widespread lineage of A. cephallonica, a sign of population growth was detected about 300 Kya ago and the lineage has been stable since the last 100 Kya (Fig. 8b).

Fig. 6
figure 6

Mismatch distributions (MD) and Bayesian skyline plots (BSP) of Anguis fragilis lineages distributed in the Balkans. a Illyrian-Central European + South Balkan haplogroups together; (b) South Balkan haplogroup separately; (c) Carniolan + Alpine-Pannonian + North Adriatic (altogether = ‘Slovenian’) haplogroups together

Fig. 7
figure 7

Mismatch distributions (MD) and Bayesian skyline plots (BSP) of Anguis colchica lineages distributed in the Balkans. a Incerta clade; (b) Carpathian lineage of the Incerta clade; (c) Stara-Planina lineage of the Incerta clade; (d) Pontic clade

Fig. 8
figure 8

Mismatch distributions (MD) and Bayesian skyline plots (BSP) of Anguis graeca (a) and A. cephallonica (b)

The complementary mismatch distributions (MDs; Figs. 6, 7 and 8) showed a ragged distribution of the observed values of pairwise differences in the predominantly Slovenian A. fragilis (Fig. 6c), the Incerta clade of A. colchica (Fig. 7a) and its Carpathian lineage (Fig. 7b), and to some extent also in A. graeca and the widespread lineage of A. cephallonica (Fig. 8a, b). In the other analysed population groups the observed values mirrored the values expected for a growing- or declining-population model. The neutrality tests showed significant departures from the neutrality in the majority of our groups, except for the predominantly Slovenian haplogroups of A. fragilis, the Incerta clade of A. colchica and its Carpathian lineage (Table 1).

Genetic diversity and topographic heterogeneity

Multiple linear regressions of nucleotide diversity (π) of the lineages/haplogroups plotted against the Q3, median above Q3, and modus above Q3 of the terrain ruggedness index (TRI), and an area size inhabited by these lineages/haplogroups, were statistically significant (Table 2). Partial regression analyses, however, revealed that only TRI values, not the area size, had a significant effect on the nucleotide diversity (Table 2, Fig. 9 and Additional file 2: Table S2, Additional file 3: Figure S1 and Additional file 4: Table S3). Standardized (beta) regression coefficients were highly significant both when samples from the Apuseni Mts. were included among the Carpathian samples as well as when they were treated separately.

Table 2 Results of the multiple linear regressions between nucleotide diversity (π), topographic heterogeneity [estimated as the third quartile (Q3) of the terrain ruggedness index (TRI), and median and modus calculated for data above Q3], and the area size of the topographic units inhabited by particular slow-worm lineages/haplogroups
Fig. 9
figure 9

Linear regressions between nucleotide diversity (π) of the Balkan slow-worm evolutionary lineages/haplogroups and modus above the third quartile of the terrain ruggedness index of particular topographic units (mountain systems). In the first analysis (a) the Apuseni Mts. were treated as a separate unit, while in the second analysis (b) the Apuseni Mts. were considered to be a part of the Carpathians. b - regression coefficient, P - probability value. Legends: 1 – Apuseni Mts., 2 – Stara Planina Mts., 3 – Macedonian-Thracian Massif, 4 – Dinarides, 5 – Carpathians, 6 – Prealps, 7 – Peloponnese, 8 – Hellenides (without Peloponnese)

Discussion

Distribution of slow worms in the Balkans and contact zones

Due to relatively hard-to-interpret morphology and description of several vaguely defined forms and their intermediates in the Balkan Peninsula, the distribution of slow worms remained problematic and conflicting [18, 22, 44]. Recent molecular-phylogenetic studies [14, 15] recognized four species of the genus Anguis within the Balkans and have also painted the first coarse-grained picture of their distribution, but the precise ranges have remained to be revealed. Here, based on extensive sampling and molecular identification, we show detailed distribution of all four species inhabiting the Balkan Peninsula (Figs. 1 and 10). The Balkan slow worms are characterized by mostly parapatric distributions, to large extent corresponding with major geomorphological units of the peninsula. We acknowledge that the distribution patterns revealed here may not fully represent species distributions due to the specific characteristics of the used mtDNA marker (maternal and clonal inheritance, reduced effective population size, sex-specific dispersal, relatively common interspecific introgression). However, the overall phylogenetic patterns we found are vastly concordant to previously published ones based on both mtDNA and nuDNA markers [14, 15].

Fig. 10
figure 10

Pleistocene refugia (R) and proposed dispersal postglacial routes of slow worms in the Balkans. Approximate species distributions given in colour shading correspond to the colour code in Fig. 1. Question marks denote missing distribution data

Among our studied species, Anguis cephallonica occupies the smallest range limited to the Peloponnese Peninsula and the islands of Kephallonia, Ithaki, and Zakynthos [17, 42, 43]. The distributions of the other three species principally follows the main mountain ranges in the Balkans; A. fragilis is distributed in the Dinarides and Macedonian-Thracian Massif, A. colchica in the Carpathians, Stara Planina, and Strandzha Mts., and A. graeca in the Hellenides.

It appears that while the ranges of A. fragilis and A. graeca each meet with ranges of two other species in the Balkans (furthermore, A. fragilis also forms a contact zone with A. veronensis outside the Balkan Peninsula, see [15]), A. colchica and A. cephallonica only come into contact with one other species. Parts of the contact zones presumably originated by crossing natural barriers such as mountain ridges or river valleys. For instance, the range of A. graeca crosses the Vardar River valley and extends from the Hellenides into the Macedonian-Thracian Massif where it forms a contact zone with A. fragilis. On the other hand, A. fragilis inhabiting predominantly the Dinarides, Macedonian-Thracian Massif, and their vicinity seems to have extended its range to the south, across the northern borderline of the Hellenides, where it forms a contact zone with A. graeca (Fig. 1). Historical demographic model indicates that an expansion of the SB haplogroup of A. fragilis could probably have happened during the Holocene (Fig. 6b). Anguis graeca and A. cephallonica form a contact zone and partial sympatry in the northern Peloponnese where they both might have come into contact repeatedly as climatic oscillations and resulting sea-level changes led to repeated connection and disconnection of the peninsula and the mainland during the Pleistocene [17, 45]. It seems that ranges of A. graeca and A. colchica do not come into recent contact because A. fragilis populations are embedded between them.

Multiple refugia and colonization routes

All four species of slow worms show high levels of intraspecific genetic differentiation in the Balkans and are sub-structured into several divergent lineages or haplogroups. This genetic structure was shaped by local restrictions of ranges into multiple Pleistocene refugia located in the Peloponnese (A. cephallonica), Hellenides (A. graeca), southern Carpathians (A. colchica), and northwestern Dinarides (A. fragilis) (Fig. 10). Existence of several smaller and isolated refugia that harboured slow-worm populations during the Pleistocene climatic oscillations within the Balkans is in concordance with the refugia-within-refugia model originally proposed for the Iberian Peninsula [46], and also suggested for the Italian Peninsula based on the phylogeography of A. veronensis [15]. This pattern might have more general applicability in the Balkans where multiple refugia were corroborated in both animals e.g. [8, 4749] and plants e.g. [50, 51]. They are located either in the Mediterranean region (e.g. Adriatic coast, Peloponnese; [52]) or in non-Mediterranean parts of the peninsula (Carpathians, and the Prealps region between the Dinarides and Alps; [4, 5, 53].

The biogeographical histories of slow worms from southern and northern Balkan refugia differ. The ICE haplogroup of A. fragilis and several haplogroups of the Carpathian lineage of A. colchica colonized broad areas of temperate Europe from their northern extra-Mediterranean refugia. On the contrary, A. cephallonica, A. graeca, the Pontic and Stara-Planina lineages of A. colchica, and the South Balkan haplogroup of A. fragilis did not disperse much from their southern Mediterranean refugia and their distribution has remained more localized south of the Danube River (Fig. 10).

In the case of Anguis fragilis our results indicate the existence of at least three separate Pleistocene refugia. The South Balkan haplogroup predominantly occurs in the Macedonian-Thracian Massif, where a refugium was presumably located. Outside this mountain range the SB haplogroup only dispersed to the northernmost Hellenides, probably recently, as a common and widespread haplotype was detected there. Populations of the ICE haplogroup colonized vast parts of the western Balkans, but also central and northwestern Europe from a refugium presumably located in the Dinarides. This happened relatively rapidly, which is indicated by (i) a star-like pattern of the haplotype network and low genetic variation of the ICE haplogroup and (ii) the broad area presumably colonized from a single source population [54]. The situation could be vividly illustrated using f1 haplotype (Fig. 2a, c): it is found not only throughout the central and western Balkans, but also in central Europe and as far as the British Isles spread over an area of approximate length of 2000 km [14]. The pattern of the haplotype network and current distribution of A. fragilis suggests not only quick expansion to the north, but also a gradual north-to-south/west-to-east expansion during the Pleistocene, which is very rare in terrestrial animals (Fig. 10; [5557]).

We detected relatively high haplotype diversity of A. fragilis in the northern Adriatic region (mainly in Slovenia; Fig. 2b, c). Also the BSP analysis demonstrated population stability for these ‘Slovenian’ haplogroups indicating a long-term survival of slow-worm populations in this region. Such persistence in refugia at foothills of the Alps has been described in several temperate amphibian and reptile species e.g. [5865]. This region was also probably important in shaping genetic diversity of A. veronensis, the species whose main part of the distribution range is located in the Apennine Peninsula [15]. However, the Prealpine slow-worm populations also contributed to the colonization of the Pannonian Basin as indicated by the phylogeographical pattern when extralimital samples were included (Fig. 2a, c; haplotype AF05; [16]).

The Carpathians formed an important extra-Mediterranean refugium of many temperate and cold-adapted species e.g. [40, 56, 66, 67]. This was mainly possible because most of the mountain range remained ice-free during the last glacial maximum [68]. In some taxa, distinct phylogenetic lineages have been detected with distribution restricted to the Carpathians, which indicates their long-term in situ survival (e.g. the newt Lissotriton vulgaris, [49, 69]; the toad Bombina variegata, [70]). These populations also contributed to the postglacial colonization of Europe. In the Carpathians or their close vicinity we discovered haplotypes of three geographically well-separated lineages of A. colchica (Stara-Planina, Banatian, and Carpathian lineages within the Incerta clade; Fig. 3). While the Stara-Planina lineage (which is currently also present in the Serbian Carpathians) presumably survived in a refugium outside the Carpathians, the Carpathian and Banatian lineages are together comprised of several haplogroups that could be traced to multiple microrefugia within the Carpathians. Close affinity of these haplotypes (or even identity in some cases, e.g. haplotypes c1, c6) to those from central and north-eastern Europe [14, 15] suggests that these areas were historically colonized from the Carpathian refugia. A very similar colonization pattern of the northern and eastern Europe from the Romanian Carpathians has been described in a rodent Clethrionomys glareolus [71].

Despite its limited distribution in the Balkans, the Pontic clade of A. colchica shows relatively high mtDNA polymorphism. Close phylogenetic relationships of the southeast Bulgarian and Anatolian populations (own unpublished data) indicate that the Pontic lineage might have colonized the Black Sea region of the Balkans during the Pleistocene when the peninsula was accessible from northern Anatolia via terrestrial route [72, 73].

The Peloponnese, inhabited by endemic A. cephallonica, and the region west of the Pindus Mts. (with high haplotype diversity of A. graeca) have favourable geography with deep long valleys providing stable climatic conditions. Consequently it is known for high endemism of numerous plants, invertebrates, and vertebrates [17, 7476]. Multiple refugia in the region have already been proposed [52]. Further in the north, most of Albania and northwestern Greece are surrounded by mountain ranges characterized by steep slopes and deep valleys which could have had a strong isolating effect on A. graeca during the Plio-Pleistocene and allowed divergence of its lineages. In contrast, the overall flat Skadar region enabled colonization of southern parts of the region of present-day Montenegro and forming a narrow zone of sympatric occurrence with A. fragilis. The existence of several distinct haplogroups in A. graeca indicates that this species has a longer and complex evolutionary history. Overall high intraspecific genetic diversity with up to 3.6 % in p-distances (Additional file 5: Table S4) suggests older diversification events probably associated with multiple refugia, e.g. in central and southern Albania, northwestern Greece, and northern Peloponnese where the most divergent haplotypes were found.

Correlation of genetic diversity and topographic heterogeneity

Phylogeographical analysis of all Balkan slow-worm species showed different patterns of intraspecific divergences and genetic diversity for each studied species, presumably mirroring their different, contrasting, evolutionary histories. Specifically, lineages with more pronounced genetic structure inhabit landscapes with higher terrain ruggedness, i.e. higher altitudinal differences, more numerous and deeper valleys, and steeper slopes. Our regression analysis indeed confirms this pattern with high significance – lineages with higher nucleotide diversity inhabit mountain systems characterized by higher elevational differences, i.e. rugged terrain (Table 2, Fig. 9 and Additional file 2: Table S2, Additional file 3: Figure S1 and Additional file 4: Table S3).

The general pattern described as southern richness and northern purity [3] is typical for many taxa on a broad continental scale and can also be observed in slow worms: the species with highest genetic diversity are A. graeca and A. cephallonica inhabiting the very south of the genus range in the Balkans. A detailed view reveals that even within the relatively small ranges of these species, the highest diversity can be found in smaller and more southerly located areas, corresponding to local microrefugia (or refugia within refugia; [46]). The situation is however different for the two northerly occurring species, A. colchica and A. fragilis, in which the populations with highest diversity occur in more northerly-located areas in the Balkans. More pronounced altitudinal differences, steep exposed slopes, and generally more heterogeneous landscapes create numerous effective barriers preventing dispersal of small legless lizards, such as slow worms, in which the dispersal ability is also limited by semifossorial lifestyle [77]. Such combination of life history and habitat characteristics provides suitable predispositions for isolation and subsequent divergence of populations. On the other hand, lowlands, plains, low-hill regions and slightly rolling landscapes offer fewer barriers to dispersal and gene flow, and thus divergence occurs less often. Our observations on correlation of slow-worm genetic diversity with topographic ruggedness are fully in concordance with the fact that 33 (63 %) of the 52 identified Mediterranean refugia are situated in submontane and montane areas [52].

Conclusions

Our study uncovered mitochondrial DNA variation and distribution of four Anguis species and hidden diversity of their populations in the Balkans. These species have mostly parapatric distributions that correspond with major mountain ranges. We showed that biogeography of the genus in the Balkans is concordant with the refugia-within-refugia model previously proposed for other main European Peninsulas. The role of Mediterranean as well as extra-Mediterranean refugia was detected in the evolutionary history of slow worms with varying ages and degrees of post-glacial recolonization. Beside climatic historical events, we consider the complex topography of the Balkans as one of the most important factors in shaping recent genetic diversity of slow worms. Topographic heterogeneity seems to be a good predictor of both genetic and species diversity, in general. The pattern observed on slow-worm refugia in the Balkan Peninsula thus illustrates and highlights the fact that many global biodiversity hotspots and endemism centres are located in montane regions [7882]. As it has been suggested in other taxa [13, 8386], complex mountain topography offers conditions that could facilitate genetic isolation and divergence and result thus in a high rate of speciation.

Availability of data and materials

New sequences have been deposited in GenBank (accession numbers KX020147–KX020322) and other input data are provided in Additional file 1: Table S1, Additional file 2: Table S2, Additional file 3: Figure S1, Additional file 4: Table S3 and Additional file 5: Table S4 of this study.

Ethics

Not applicable.

Consent to publish.

Not applicable.

References

  1. Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF. Comparative phylogeography and postglacial colonization routes in Europe. Mol Ecol. 1998;7:453–64.

    Article  CAS  PubMed  Google Scholar 

  2. Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999;68:87–112.

    Article  Google Scholar 

  3. Hewitt GM. The genetic legacy of the Quarternary ice ages. Nature. 2000;405:907–13.

    Article  CAS  PubMed  Google Scholar 

  4. Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Philos Trans R Soc Lond B Biol Sci. 2010;277:661–71.

    Article  Google Scholar 

  5. Schmitt T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007;4:11.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Gaston KJ, David R. Hotspots across Europe. Biodiversity Letters. 1994;2:108–16.

    Article  Google Scholar 

  7. Džukić G, Kalezić ML. The Biodiversity of Amphibians and Reptiles in the Balkan Peninsula. In: Griffiths HI, editor. Balkan Biodiversity. Kluwer Academic Publishers: Dordrecht; 2004. p. 167–792.

    Google Scholar 

  8. Hewitt GM. Mediterranean Peninsulas: The Evolution of Hotspots. In: Zachos FE, Habel JC, editors. Biodiversity Hotspots. Springer Publishers: Berlin Heidelberg; 2011. p. 123–47.

    Chapter  Google Scholar 

  9. Schmitt T, Varga Z. Extra-Mediterranean refugia: The rule and not the exception? Front Zool. 2012;9:22.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Ager DV. The geology of Europe. London: McGraw-Hill; 1980.

    Google Scholar 

  11. Reed JM, Kryštufek B, Eastwood WJ. The physical geography of the Balkans and nomenclature of place names. In: Griffiths HI, editor. Balkan Biodiversity. Dordrecht: Kluwer Academic Publishers; 2004. p. 9–22.

    Chapter  Google Scholar 

  12. McRae BH. Isolation by resistance. Evolution. 2006;60:1551–61.

    Article  PubMed  Google Scholar 

  13. Guarnizo CE, Cannatella DC. Genetic divergence within frog species is greater in topographically more complex regions. J Zool Sys Evol Res. 2013;51:333–40.

    Google Scholar 

  14. Gvoždík V, Jandzik D, Lymberakis P, Jablonski D, Moravec J. Slow Worm, Anguis fragilis (Reptilia: Anguidae) as a species complex: Genetic structure reveals deep divergences. Mol Phylogenet Evol. 2010;55:460–72.

    Article  PubMed  Google Scholar 

  15. Gvoždík V, Benkovský N, Crottini A, Bellati A, Moravec J, Romano A, et al. An ancient lineage of slow worms, genus Anguis (Squamata: Anguidae), survived in the Italian Peninsula. Mol Phylogenet Evol. 2013;69:1077–92.

    Article  PubMed  Google Scholar 

  16. Szabó K, Vörös J. Distribution and hybridization of Anguis fragilis and A. colchica in Hungary. Amphibia-Reptilia. 2014;35:135–40.

    Article  Google Scholar 

  17. Thanou E, Giokas S, Kornilios P. Phylogeography and genetic structure of the slow worms Anguis cephallonica and Anguis graeca (Squamata: Anguidae) from the southern Balkan Peninsula. Amphibia-Reptilia. 2014;35:263–9.

    Article  Google Scholar 

  18. Dely OG. Anguis fragilis Linnaeus 1758 – Blindschleiche. In: Böhme W, editor. Handbuch der Reptilien und Amphibien Europas. Band 1. Echsen (Sauria) 1. Wiesbaden: AULA-Verlag; 1981. p. 241–58.

    Google Scholar 

  19. Sillero N, Campos J, Bonardi A, Corti C, Creemers R, Crochet P-A, et al. Updated distribution and biogeography of amphibians and reptiles of Europe. Amphibia-Reptilia. 2014;35:1–31.

    Article  Google Scholar 

  20. Stumpel AHP. Biometrical and ecological data from a Netherlands population of Anguis fragilis (Reptilia, Sauria, Anguidae). Amphibia-Reptilia. 1985;6:181–94.

    Article  Google Scholar 

  21. Hubble DS, Hurst DT. Population structure and translocation of the Slow-worm, Anguis fragilis L. Herpetological Bulletin. 2006;97:8–13.

    Google Scholar 

  22. Džukić G. Taxonomic and biogeographic characteristics of the slow-worm (Anguis fragilis Linnaeus, 1758) in Yugoslavia and on the Balcan Peninsula. Scopolia. 1987;12:1–47.

    Google Scholar 

  23. Stojanov A, Tzankov N, Naumov B. Die Amphibien und Reptilien Bulgariens. Frankfurt am Main: Edition Chimaira; 2011.

    Google Scholar 

  24. Librado P, Rozas J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.

    Article  CAS  PubMed  Google Scholar 

  25. Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29:1695–701.

    Article  CAS  PubMed  Google Scholar 

  26. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61:539–42.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Rambaut A, Suchard MA, Xie W, Drummond A. Tracer v 1.6. 2013. http://beast.bio.ed.ac.uk/Tracer. Accessed 1 Feb 2016.

  29. Posada D, Crandall KA. Intraspecific gene genealogies: trees grafting into networks. Trends Ecol Evol. 2001;16:37–45.

    Article  PubMed  Google Scholar 

  30. Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–9.

    Article  CAS  PubMed  Google Scholar 

  31. Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005;22:1185–92.

    Article  CAS  PubMed  Google Scholar 

  32. Bouckaert R, Heled J, Kühner D, Vaughan T, Wu CH, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10, e1003537.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Macey RJ, Schulte JA, Larson A, Tuniyev BS, Orlov N, Papenfuss TJ. Molecular phylogenetics, tRNA evolution, and historical biogeography in Anguid lizards and related taxonomic families. Mol Phylogenet Evol. 1999;12:250–72.

    Article  CAS  PubMed  Google Scholar 

  34. Fu Y-X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:913–25.

    Google Scholar 

  35. Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. Ramos-Onsins SE, Rozas J. Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002;19:2092–100.

    Article  CAS  PubMed  Google Scholar 

  37. Riley SJ, DeGloria SD, Elliot R. A terrain ruggedness index that quantifies topographic heterogeneity. Intermt J Sci. 1999;5:23–7.

    Google Scholar 

  38. GRASS Development Team. Geographic Resources Analysis Support System (GRASS) Software, Version 7.0. 2015. Open Source Geospatial Foundation. http://grass.osgeo.org. Accessed 1 February 2016.

  39. Magri D, Vendramin GG, Comps B, Dupanloup I, Geburek T, Gömöry D, et al. A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences. New Phytol. 2006;171:199–221.

    Article  CAS  PubMed  Google Scholar 

  40. Bálint M, Ujvárosi L, Theissinger K, Lehrian S, Mészáros N, Pauls SU. The Carpathians as a Major Diversity Hotspot in Europe. In: Zachos FE, Habel JC, editors. Biodiversity Hotspots. Berlin Heidelberg: Springer Publishers; 2011. p. 189–205.

    Chapter  Google Scholar 

  41. Hill T, Lewicki P. STATISTICA: Methods and Applications. Tulsa: StatSoft; 2007.

    Google Scholar 

  42. Grillitsch H, Cabela A. Zum systematischen Status der Blindschleichen (Squamata: Anguidae) der Peloponnes und der südlichen Ionischen Inseln (Griechenland). Herpetozoa. 1990;2:131–53.

    Google Scholar 

  43. Mayer W, Grillitsch H, Cabela A. Proteinelektrophoretische Untersuchungen zur Systematik der südgriechischen Blindschleiche (Squamata: Anguidae). Herpetozoa. 1991;4:157–65.

    Google Scholar 

  44. Cabela A, Grillitsch H. Zum systematischen Status der Blindschleiche (Anguis fragilis Linnaeus, 1758) von Nordgriechenland und Albanien (Squamata: Anguidae). Herpetozoa. 1989;2:51–69.

    Google Scholar 

  45. Lymberakis P, Poulakakis N. Three Continents Claiming an Archipelago: The Evolution of Aegean’s Herpetofaunal Diversity. Diversity. 2010;10:233–55.

    Article  Google Scholar 

  46. Gómez A, Lunt DH. Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography in Southern European Refugia. In: Weiss S, Ferrand N, editors. Evolutionary Perspectives on Origins and Conservation of European Biodiversity. Dordrech: Springer Publishers; 2007. p. 155–88.

    Google Scholar 

  47. Ursenbacher S, Schweiger S, Tomović L, Crnobrnja-Isailović J, Fumagalli L, Mayer W. Molecular phylogeography of the nose-horned viper (Vipera ammodytes, Linnaeus (1758)): Evidence for high genetic diversity and multiple refugia in the Balkan peninsula. Mol Phylogenet Evol. 2008;46:1116–28.

    Article  CAS  PubMed  Google Scholar 

  48. Podnar M, Mađarić BB, Mayer W. Non-concordant phylogeographical patterns of three widely codistributed endemic Western Balkans lacertid lizards (Reptilia, Lacertidae) shaped by specific habitat requirements and different responses to Pleistocene climatic oscillations. J Zool Syst Evol Res. 2014;52:119–29.

    Article  Google Scholar 

  49. Pabijan M, Zieliński P, Dudek K, Chloupek M, Sotiropoulos K, Liana M, et al. The dissection of a Pleistocene refugium: phylogeography of the smooth newt, Lissotriton vulgaris, in the Balkans. J Biogeogr. 2015;42:671–83.

    Article  Google Scholar 

  50. Kramp K, Huck S, Niketić M, Tomović G, Schmitt T. Multiple glacial refugia and complex postglacial range shifts of the obligatory woodland plant Polygonatum verticillatum (Convallariaceae). Plant Biol. 2009;11:392–404.

    Article  CAS  PubMed  Google Scholar 

  51. Surina B, Schonswetter P, Schneeweiss GM. Quaternary range dynamics of ecologically divergent species (Edraianthus serpyllifolius and E. tenuifolius, Campanulaceae) within the Balkan refugium. J Biogeogr. 2011;38:1381–93.

    Article  Google Scholar 

  52. Medail F, Diadema K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J Biogeogr. 2009;36:1333–45.

    Article  Google Scholar 

  53. Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23:564–71.

    Article  PubMed  Google Scholar 

  54. Excoffier L, Foll M, Petit RJ. Genetic consequences of range expansions. Annu Rev Ecol Evol Syst. 2009;40:481–501.

    Article  Google Scholar 

  55. Deffontaine V, Libois R, Kotlík P, Sommer R, Nieberding C, Paradis E, et al. Beyond the Mediterranean peninsulas: evidence of central European glacial refugia for a temperate forest mammal species, the bank vole (Clethrionomys glareolus). Mol Ecol. 2005;14:1727–39.

    Article  CAS  PubMed  Google Scholar 

  56. Kotlík P, Deffontaine V, Mascheretti S, Zima J, Michaux JR, Searle JB. A northern glacial refugium for bank voles (Clethrionomys glareolus). Proc Natl Acad Sci U S A. 2006;103:14860–4.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Kryštufek B, Bužan EV, Hutchinson WF, Hanfling B. Phylogeography of the rare Balkan endemic Martino’s vole, Dinaromys bogdanovi, reveals strong differentiation within the western Balkan Peninsula. Mol Ecol. 2007;16:1221–32.

    Article  PubMed  Google Scholar 

  58. Canestrelli D, Cimmaruta R, Nascetti G. Phylogeography and historical demography of the Italian treefrog Hyla intermedia reveals multiple refugia, range expansions and secondary contacts within the Italian peninsula. Mol Ecol. 2007;16:4808–21.

    Article  CAS  PubMed  Google Scholar 

  59. Canestrelli D, Salvi D, Maura M, Bologna MA, Nascetti G. One species, three Pleistocene evolutionary histories: Phylogeography of the Italian crested newt, Triturus carnifex. PLoS One. 2012;7:e41754.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Sotiropoulos K, Eleftherakos K, Džukić G, Kalezić ML, Legakis A, Polymeni RM. Phylogeny and biogeography of the alpine newt Mesotriton alpestris (Salamandridae, Caudata), inferred from mtDNA sequences. Mol Phylogenet Evol. 2007;45:211–26.

    Article  CAS  PubMed  Google Scholar 

  61. Crottini A, Andreone F, Kosuch J, Borkin LJ, Litvinchuk SN, Eggert C, et al. Fossorial but widespread: the phylogeography of the common spadefoot toad (Pelobates fuscus), and the role of the Po Valley as a major source of genetic variability. Mol Ecol. 2007;16:2734–54.

    Article  PubMed  Google Scholar 

  62. Canestrelli D, Nascetti G. Phylogeography of the pond frog Rana (Pelophylax) lessonae in the Italian peninsula and Sicily: multiple refugia, glacial expansions and nuclear-mitochondrial discordance. J Biogeogr. 2008;35:1923–36.

    Article  Google Scholar 

  63. Salvi D, Harris DJ, Kaliontzopoulou A, Carretero MA, Pinho C. Persistence across Pleistocene ice ages in Mediterranean and extra-Mediterranean refugia: phylogeographic insights from the common wall lizard. BMC Evol Biol. 2013;13:147.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Kindler C, Böhme W, Corti C, Gvoždík V, Jablonski D, Jandzik D, et al. Mitochondrial phylogeography, contact zones and taxonomy of grass snakes (Natrix natrix, N. megalocephala). Zool Scripta. 2013;42:458–72.

    Article  Google Scholar 

  65. Maura M, Salvi D, Bologna MA, Nascetti G, Canestrelli D. Northern richness and cryptic refugia: Phylogeography of the Italian smooth newt Lissotriton vulgaris meridionalis. Biol J Linn Soc. 2014;113:590–603.

    Article  Google Scholar 

  66. Stojak J, Mcdevitt AD, Herman JS, Searle JB, Wójcik JM. Post-glacial colonization of eastern Europe from the Carpathian refugium: evidence from mitochondrial DNA of the common vole Microtus arvalis. Biol J Linn Soc. 2015;115:927–39.

    Article  Google Scholar 

  67. Wielstra B, Babik W, Arntzen JW. The crested newt Triturus cristatus recolonized temperate Eurasia from an extra-Mediterranean glacial refugium. Biol J Linn Soc. 2015;115:574–87.

    Article  Google Scholar 

  68. Reuther AU, Urdea P, Geiger C, Ivy-Ochs S, Niller H-P, Kubik PW, et al. Late Pleistocene glacial chronology of the Pietrele Valley, Retezat Mountains, Southern Carpathians constrained by 10Be exposure ages and pedological investigations. Quatern Int. 2007;164–165:151–69.

    Article  Google Scholar 

  69. Babik W, Branicki W, Crnobrnja-Isailović J, Cogălniceanu D, Sas I, Olgun K, et al. Phylogeography of two European newt species – discordance between mtDNA and morphology. Mol Ecol. 2005;14:2475–91.

    Article  CAS  PubMed  Google Scholar 

  70. Fijarczyk A, Nadachowska K, Hofman S, Litvinchuk SN, Babik W, Stuglik M, et al. Nuclear and mitochondrial phylogeography of the European fire-bellied toads Bombina bombina and Bombina variegata supports their independent histories. Mol Ecol. 2011;20:3381–98.

    Article  PubMed  Google Scholar 

  71. Wójcik JM, Kawałko A, Marková S, Searle JB, Kotlík P. Phylogeographic signatures of northward post-glacial colonization from high-latitude refugia: a case study of bank voles using museum specimens. J Zool. 2010;281:249–62.

    Google Scholar 

  72. Kerey IE, Meric E, Tunoglu C, Kelling G, Brenner RL, Dogan AU. Black Sea–Marmara Sea Quaternary connections: new data from the Bosphorus, Istanbul, Turkey. Palaeogeogr Palaeoclimatol Palaeoecol. 2004;204:277–95.

    Article  Google Scholar 

  73. Wielstra B, Espregueira Themudo G, Güçlü Ö, Olgun K, Poyarkov NA, Arntzen JW. Cryptic crested newt diversity at the Eurasian transition: The mitochondrial DNA phylogeography of Near Eastern Triturus newts. Mol Phylogenet Evol. 2010;56:888–96.

    Article  CAS  PubMed  Google Scholar 

  74. Jesse R, Shubart CD, Klaus S. Identification of cryptic lineage within Potamon fluviatile (Herbst) (Crustacea: Brachyura: Potamidae). Invertebrate Systematics. 2010;24:348–56.

    Article  Google Scholar 

  75. Ferchaud A-L, Ursenbacher S, Cheylan M, Luiselli L, Jelić D, Halpern B, et al. Phylogeography of the Vipera ursinii complex (Viperidae): mitochondrial markers reveal an east–west disjunction in the Palaearctic region. J Biogeogr. 2012;39:1836–47.

    Article  Google Scholar 

  76. Poulakakis N, Kapli P, Lymberakis P, Trichas A, Vardinoyiannis K, Sfenthourakis S, et al. A review of phylogeographic analyses of animal taxa from the Aegean and surrounding regions. J Zool Sys Evol Res. 2015;53:18–32.

    Article  Google Scholar 

  77. Haley T. A metapopulation of the lizard Anguis fragilis (Squamata: Anguidae) on a local scale in Dorset, Great Britain, as indicated by spatial distribution and movement. Phyllomedusa. 2014;13:91–8.

    Article  Google Scholar 

  78. Mittermeier RA, Gil PR, Hoffmann M, Pilgrim J, Brooks T, Mittermeier CG, et al. Hotspots Revisited: Earth’s Biologically Richest and Most Endangered Terrestrial Ecoregions, Conservation International. Washington: The University of Chicago Press; 2005.

    Google Scholar 

  79. Orme C, Davies R, Burgess M, Eigenbrod F, Pickup N, Olson VA, et al. Global hotspots of species richness are not congruent with endemism or threat. Nature. 2005;436:1016–9.

    Article  CAS  PubMed  Google Scholar 

  80. Rull V. Biotic diversification in the Guayana Highlands: a proposal. J Biogeogr. 2005;32:921–7.

    Article  Google Scholar 

  81. Meegaskumbura M, Bossuyt F, Pethiyagoda R, Manamendra-Arachchi K, Bahir M, Milinkovitch MC, et al. Sri Lanka: An Amphibian Hot Spot. Science. 2002;298:379.

    Article  CAS  PubMed  Google Scholar 

  82. Mittermeier RA, Myers N, Thomsen JB, Da Fonseca GA, Olivieri S. Biodiversity hotspots and major tropical wilderness areas: approaches to setting conservation priorities. Conserv Biol. 1998;12:516–20.

    Article  Google Scholar 

  83. Funk WC, Blouin MS, Corn PS, Maxell BA, Pilliod DS, Amish S, et al. Population structure of Columbia spotted frogs (Rana luteiventris) is strongly affected by the landscape. Mol Ecol. 2005;14:483–96.

    Article  CAS  PubMed  Google Scholar 

  84. Giordano AR, Ridenhour BJ, Storfer A. The influence of altitude and topography on genetic structure in the long-toed salamander (Ambystoma macrodactylum). Mol Ecol. 2007;16:1625–37.

    Article  CAS  PubMed  Google Scholar 

  85. Badgley C. Tectonics, topography, and mammalian diversity. Ecography. 2010;33:220–31.

    Google Scholar 

  86. Rodríguez A, Börner M, Pabijan M, Gehara M, Haddad CFB, Vences M. Genetic divergence in tropical anurans: deeper phylogeographic structure in forest specialists and in topographically complex regions. Evol Ecol. 2015;29:765–85.

    Article  Google Scholar 

Download references

Acknowledgments

We would like to thank (alphabetically) to J.W. Arntzen, P. Balej, L. Blažej, B. Bolfíková, L, Choleva, D. Coğalniceanu, M. Donát, C. Dundarová, I. Ghira, S. Giokas, J. Hájek, J. Hlaváč, M. Homolka, T. Husák, K. Janoušek, M. Jedlička, A. Kecskés, D. Koleška, P. Kornilios, Z. Lajbner, H. Laťková, P. Lymberakis, P. Meduna, O. Mettouris, T. Míšek, E. Mizsei, L. Muller, R. Musilová, P. Pavlík, N. Preradović, M. Rindoš, R. Rozínek, R. Šanda, M. Šandera, M. Šandová, H. Šifrová, E. Stloukal, B. Švecová, O. Tzortzakaki, I. Velikov, P. Vlček, J. Vörös, J. Vukić, B. Wielstra, S.R. Zamfirescu, V. Zavadil for their kind donation of tissue samples or help in the field. We thank J. Kreisinger for helpful discussion about the statistics, D. Senko for GIS calculations and two anonymous reviewers for their comments, which improved previous version of the manuscript. The Geographic Resources Analysis Support System (GRASS) and GIS calculations were performed in the Computing Centre of the Slovak Academy of Sciences using the infrastructure acquired within the projects ITMS 26230120002, ITMS 26210120002 and ITMS 6240120014 supported by the Research & Development Operational Program funded by the European Regional Development Fund (ERDF). We also acknowledge the local authorities of nature conservation. This study was supported by the Comenius University grants UK/20/2014, UK/37/2015 and a grant of the Scientific Grant Agency of the Slovak Republic VEGA 1/0073/14. Daniel Jablonski was also supported by the Societas Europaea Herpetologica (Travel grant 2013). Additional financial support was from the Institute of Vertebrate Biology, Czech Academy of Sciences (RVO 68081766), and the Ministry of Culture of the Czech Republic (DKRVO 2016/15, National Museum, 00023272). The work of K. Ljubisavljević and G. Džukić was supported by the Serbian Ministry of Education, Science and Technological Development (grant173043).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Daniel Jablonski.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

DJab, DJan, PM, JM, and VG conceived the ideas, designed the project and interpreted the obtained data; DJab, VG performed statistical analyses; DJab, DJan., and VG. wrote the manuscript; all authors collected field samples and revised the manuscript. All authors read and approved the final manuscript.

Authors’ information

Daniel Jablonski is a PhD candidate at the Comenius University in Bratislava where he studies evolutionary biology and historical biogeography of amphibians and reptiles.

Additional files

Additional file 1: Table S1.

A list of samples, their coordinates, locality numbers in maps (Figs. 2, 3, 4 and 5), and GenBank accession numbers. Sample IDs in bold were already used earlier (Gvoždík et al. [14, 15]). (PDF 121 kb)

Additional file 2: Table S2.

Samples used in regression analyses of nucleotide diversity (π) and terrain ruggedness index (TRI), and their assignment to particular topographic mountain units (mountain ranges). (PDF 27 kb)

Additional file 3: Figure S1.

A map of demarcated topographic units as defined for regression analyses of nucleotide diversity (π) and terrain ruggedness index (TRI). (PDF 131 kb)

Additional file 4: Table S3.

Values of nucleotide diversity (π), area size of particular topographic units (in km2), and terrain ruggedness index (TRI). The third quartile (TRI Q3), and median (TRI Q3 median) and modus (TRI Q3 modus) of data above TRI Q3 were used in regression analyses. (PDF 21 kb)

Additional file 5: Table S4.

Average uncorrected p-distances calculated among the main evolutionary lineages within each of the four Anguis species distributed in the Balkans. The highest values are in bold. (PDF 13 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jablonski, D., Jandzik, D., Mikulíček, P. et al. Contrasting evolutionary histories of the legless lizards slow worms (Anguis) shaped by the topography of the Balkan Peninsula. BMC Evol Biol 16, 99 (2016). https://doi.org/10.1186/s12862-016-0669-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-016-0669-1

Keywords