Skip to main content

Biogeographic and demographic history of the Mediterranean snakes Malpolon monspessulanus and Hemorrhois hippocrepis across the Strait of Gibraltar



The contribution of North Africa to the assembly of biodiversity within the Western Palaearctic is still poorly documented. Since the Miocene, multiple biotic exchanges occurred across the Strait of Gibraltar, underlying the high biogeographic affinity between the western European and African sides of the Mediterranean basin. We investigated the biogeographic and demographic dynamics of two large Mediterranean-adapted snakes across the Strait and assess their relevance to the origin and diversity patterns of current European and North African populations.


We inferred phylogeographic patterns and demographic history of M. monspessulanus and H. hippocrepis, based on range-wide multilocus data, combined with fossil data and species distribution modelling, under present and past bioclimatic envelopes. For both species we identified endemic lineages in the High Atlas Mountains (Morocco) and in eastern Iberia, suggesting their persistence in Europe during the Pleistocene. One lineage is shared between North Africa and southern Iberia and likely spread from the former to the latter during the sea-level low stand of the last glacial stage. During this period M. monspessulanus shows a sudden demographic expansion, associated with increased habitat suitability in North Africa. Lower habitat suitability is predicted for both species during interglacial stages, with suitable areas restricted to coastal and mountain ranges of Iberia and Morocco. Compiled fossil data for M. monspessulanus show a continuous fossil record in Iberia at least since the Pliocene and throughout the Pleistocene.


The previously proposed hypothesis of Pleistocene glacial extinction of both species in Europe is not supported based on genetic data, bioclimatic envelopes models, and the available fossil record. A model of range retraction to mountain refugia during arid periods and of glacial expansion (demographic and spatial) associated to an increase of Mediterranean habitats during glacial epochs emerges as a general pattern for mesic vertebrates in North Africa. Moreover, the phylogeographic pattern of H. hippocrepis conforms to a well-established biogeographic partition between western and eastern Maghreb.

Peer Review reports


The Mediterranean basin is a global biodiversity hotspot and represents the major ecoregion within the Western Palearctic realm [1, 2]. The prominent role of southern Mediterranean peninsulas as Pleistocene glacial refugia and differentiation centres for Western Palaearctic thermophilic species has been corroborated by a vast amount of fossil and phylogeographic data [3, 4]. Far less is known about the occurrence and location of refugia in North Africa and their contribution to the assembly of Mediterranean biotas on both sides of the Gibraltar Strait [4,5,6,7]. Moreover, while for Europe it is well established that glacial periods represented a phase of demographic (and range) contraction for thermophilic species, with expansion phases associated with interglacial periods [3, 4, 8], less is known regarding the demographic and range dynamics of North African species during Pleistocene climatic cycles. These species might have experience favourable environmental conditions during glacial stages, as a consequence of the expansion of Mediterranean habitats [9], but this is still not well documented.

The biogeographic affinities between the European and North African portions of the Mediterranean basin are particularly accentuated across the Strait of Gibraltar, where the distance between the two regions is relatively short (~ 14 km), leading to the inclusion of Iberia and the northern part of Maghreb in the same biogeographic sub-region, the Atlantic-Mediterranean sub-centre [7, 10]. Indeed, a significant portion of the Iberian Peninsula and the Maghreb share similar climates, integrated within the Mediterranean biome, that are suitable for similar taxa in both regions [1, 7]. Since the Miocene, the complex history of intermittent connections at the Strait of Gibraltar has created opportunities for dispersal and vicariance events, that have contributed to the biotic exchange between North Africa and Western Europe [11,12,13,14]. The Messinian Salinity crisis (5.9–5.33 Ma) is often seen as a major period for biotic exchanges between Africa and Europe, due to their wide land connection, although during the Tortonian (~ 10–6 Ma), tectonic uplift and volcanism caused progressive land emergence and reduction in seaways that might have acted as filter bridges for earlier terrestrial fauna exchange [14, 15]. The opening of the Strait at the end of the Messinian Salinity crisis (5.3 million year ago, Mya) [16] established a biogeographical barrier for terrestrial organisms and isolated the Iberian and North African biotas. However, Plio-Pleistocene transmarine colonisations across the Strait have also been reported, often associated with glacial sea- level low stands [11, 17,18,19]. For example, in reptiles many species show shallow genetic differentiation between “European” and “African” populations, suggesting that they have independently crossed the Strait of Gibraltar, in both directions and on different occasions, throughout the Pleistocene [14, 20,21,22]. For many of these species the relative role of Iberia and North Africa as glacial or interglacial refugia is still unclear, as well as their role as sources of colonization toward the other side of the Strait.

Two large Mediterranean colubrid snakes, the Montpellier snake, Malpolon monspessulanus, and the horseshoe whip snake, Hemorrhois hippocrepis, present fossil and mitochondrial genetic evidence of a complex biogeographic history across the Strait of Gibraltar and suggest a North African origin of the current European populations [19]. However, the pattern of shallow divergence found in both species did not allow inference of putative refugia, or to reach a firm conclusion concerning the timing of colonization of Europe. The lack of a dense sampling and multilocus data may have prevented the identification of areas harbouring high genetic diversity and stable climatic suitability through time (putative refugia) as well as the estimate of the timing of the demographic and spatial expansions of both species across the Strait. It remains unclear whether both species became extinct in Europe during Pleistocene glaciations and only recently recolonized Europe from North Africa as suggested by Carranza et al. (2006) [19], or instead they have been present both in Europe and North Africa during the Pliocene and Pleistocene as appears conceivable based on fossil records [23,24,25].

Here we use an integrative approach based on an extended geographical and genetic sampling, using both mitochondrial and nuclear DNA sequence data, climate modelling techniques and fossil data to reconstruct the demographic and evolutionary history of both species and to clarify the role of North Africa and Europe as putative refuge areas and colonization source for the species during the Pleistocene. The main goal is to understand the biogeographic and evolutionary dynamics of these Mediterranean adapted species across the Strait of Gibraltar during Pleistocene climatic oscillations, and their relevance in the assembly of current European and North African biotas.



In total, we analysed 59 specimens of Malpolon monspessulanus and 64 specimens of Hemorrhois hippocrepis, covering exhaustively the distribution range of each species. The dataset comprised newly generated sequences from specimens of the collections of Centro de Investigação em Biodiversidade e Recursos Genéticos (CIBIO/InBIO) and sequences used in [19] and retrieved from GenBank. Detailed information about samples codes, sampling locality and GenBank accession numbers can be found in Fig. 1 and Additional file 1: Table S1.

Fig. 1

Study area and sampling of Hemorrhois hippocrepis (a) and Malpolon monspessulanus (b) (Photos by D. Salvi). Dashed lines represent species native ranges (IUCN 2008—The IUCN Red List of Threatened Species). White squares: samples used in Carranza et al. (2006) [18]; black circles: original samples from the current study

DNA extraction, amplification and sequencing

Total genomic DNA was extracted from muscle tissue, using either a standard saline method [26] or the Qiagen Dneasy® Blood & Tissue extraction kit. We amplified two mitochondrial gene fragments, the cytochrome b (cyt-b; ⁓ 300 base pairs, bp) and the ribosomal 12S rRNA (12S; ⁓ 380 bp), and two nuclear gene fragments, the melanocotyin receptor 1 (MC1R; ⁓ 650 bp) and the brain-derived neurotrophic factor (BDNF; ⁓ 600 bp). These markers were successfully used in previous phylogenetic and phylogeographic studies on the same species and other Palaearctic snakes. For the cyt-b we amplified the same fragment used in [19], in order to combine our data with those from this study. Additionally, for selected samples we amplified a longer cyt-b fragments (“cyt-b long”; ⁓ 700 bp) in order to increase phylogenetic resolution. Sequenced were checked and edited in Geneious R7 [27]. Sequences of protein coding genes (cyt-b, MC1R, and BDNF) were translated into amino acids to check for the presence of stop codons. Details regarding primers and PCR protocols are provided in Table 1.

Table 1 Details on gene amplifications: primers’ name, sequence and references, and PCR conditions used

Phylogeographic analysis

Multiple sequence alignments were performed using MAFFT 7 [33] with default settings, except for the 12S alignment, for which we used the Q-INS-i-strategy that takes in to account the secondary structure of the RNA. Haplotype reconstruction for nuclear gene fragments was performed in PHASE 2.1 [34] as implemented in DNAsp 5 [35]. Haplotype (h) and nucleotide diversity (π) were also calculated in DNAsp for the concatenate mtDNA alignment (cyt-b and 12S) and each nuclear marker (phased sequences). Phylogenetic relationships between haplotypes were inferred using the parsimony network approach implemented in TCS 1.21 with a 95% connection limit [36]. Network approaches are the most appropriate for intraspecific gene evolution, particularly when few characters for phylogenetic analysis are available due to shallow levels of divergence. Haplotype networks were inferred for each gene alignment and for the concatenated mtDNA alignment.

The demographic histories of both snakes were inferred using the Extended Bayesian Skyline Plots (EBSP) implemented in BEAST 1.8.0 [37]. EBSP is a coalescent based method that estimates current and past effective population sizes (Ne) from genetic data from multiple independent loci, which significantly improves the power of detecting past demographic fluctuations and reduces estimation errors [37, 38]. EBSP were performed with the concatenated mtDNA dataset, plus the two nuclear gene datasets, implementing the best-fit models of nucleotide substitution selected by jModelTest 2.1.10 [39] (Additional file 1: Table S2). In order to estimate the time of demographic events within Malpolon monspessulanus and Hemorrhois hippocrepis we used the rate of evolution of 1.3% substitution/million years (s/my) estimated for mitochondrial DNA in colubrids by Daza, Smith, Páez, & Parkinson (2009) [40]. We defined a lognormal distribution on the ucld.mean parameter with a mean of 0.013 substitutions/my and a 95% confidence interval of [0.009–0.0179], thereby covering substitution rates estimated in previous studies on colubrids [41]. BEAST analyses were run twice for 150 million generations sampling every 5000 steps, (burn-in = 10%); run convergence was assessed using Tracer 1.5 [42] and EBSP were drawn using the package ggplot [43].

Predictive modelling

We evaluated the suitability of Iberian and North African environments for Malpolon monspessulanus and Hemorrhois hippocrepis under present, Holocene (6 thousand years ago, kya), Last Glacial Maximum (LGM; 23–18 kya) and Last Interglacial (LIG; 140–120 kya) bioclimatic envelopes, using the maximum-entropy algorithm implemented in MaxEnt 3.3.3e [44]. This presence-only modelling technique fits with the elusive nature of snakes and has a good performance with both small and large sample sizes [45]. To prepare the environmental layers we used bioclimatic data from WorldClim version 1.4 [46]. We downloaded 19 bioclimatic variables available for current and past conditions at a resolution of 5 arc-min ~ 10 km. We cut each bioclimatic raster to a study area that accommodates both species native distribution and putative range shifts across Pleistocene: Xmin = − 17.65, Ymin = 14.69, Xmax = 27.60, and Ymax = 49.1. Among all bioclimatic variables, we selected six variables that have high biological significance for snakes [47, 48] and show a Spearman's correlation lower than 0.77: bio6, Minimum Temperature of Warmest Month, bio7, Temperature Annual Range; bio 8, Mean Temperature of Wettest Quarter; bio9, Mean Temperature of Driest Quarter; bio 15, Precipitation Seasonality; bio16, Precipitation of Wettest Quarter. Using different sets of not-autocorrelated variables gave consistent results, thus suggesting that the models are robust to variable selection. In order to project suitability of the species in the past, we downloaded the same bioclimatic variables for LIG, LGM and Holocene periods from the CCSM4 (CC) Global Climate Model as well as from the MIROC-ESM (MIROC) and MPI-ESM-P (MPI) models for LGM and Holocene.

We collected 10,792 presence data points for the two species (M. monspessulanus: 7489; H. hippocrepis: 3303) from field observations [49], the GBIF [50], atlas of moroccan herpetofauna [49, 51] and from the collections of the Centre d’Ecologie Fonctionelle et Evolutive (CEFE, Montpellier), the CIBIO/InBIO, the Muséum National d’Histoire Naturelle (MNHN, Paris) and the Natural History Museum (NHM, London). To account for sampling bias, first we filtered occurrence data in order to have only one record in each 10 km square (corresponding to the resolution of bioclimate variables). Second, we filtered the occurrence dataset taking in to account the bioclimatic data along the study area. We performed a principal component analysis (PCA) using the values of the six bioclimatic variables at each locality with a presence record. Based on the PCA result, we realized a K-means cluster analysis to group each of these localities according to four different bioclimatic clusters. Having assign a bioclimatic cluster, we filter the presence data in order to maximize bioclimatic variance. All this analysis was performed in R with packages factoextra and ggplot2 [52,53,54]. The final dataset of occurrence records included a total of 325 records of H. hippocrepis and 367 records of M. monspessulanus (Additional file 1: Tables S3 and S4 and Fig. S1 for geographic details). In order to evaluate the degree of spatial autocorrelation we calculated the Nearest neighbour index (NNI) in QGIS 2.8.3 [76] for H. hippocrepis (NNI = 0.59; Z-score = − 14.25) and M. monspessulanus (NNI = 0.65; Z-score = − 12.62), that revealed a low degree of clustering in the distribution of records.

MaxEnt analyses were run with random selection of 70% presence records; the remaining 30% were used as test data on 25 replicates models. The importance of each variable for model building was assessed with jackknife analyses. Curve plots with a relationship between species presence and bioclimatic variables values were obtained. Models were also tested with received-operated characters (ROC) for the agreement between observed species presence and projected distribution. The area under the curve (AUC) of the ROC analysis provide a measure of the model performance, ranges from 0.5 (randomness) to 1 (perfect discrimination). A good model accuracy is considered with an AUC ≥ 0.7. All projections produced by MaxEnt were posteriorly edited in QGIS.

Fossil record

In order to obtain insights on the occurrence of both species in Iberia and North Africa during the last million years, we gathered fossil data for both species from the Fossil Fishes, Amphibians, Reptiles, Birds database, fosFARbase [55]. For both regions, we also retrieved a few records identified as Malpolon sp. or Coluber sp. (hippocrepis was originally included within Coluber Linnaeus, 1758) for which their attribution to either of the two study species cannot be ruled out (Hugues-Alexandre Blain, personal communication). These latter records were treated separately from those identified at the specie level.


Phylogeographic analysis

For most of the collected samples (Fig. 1) we were able to obtain sequences of the short cyt-b fragment and of nuclear markers, whereas amplification rate was lower for the long cytb-fragment and for the 12S (Additional file 1: Table S1), likely because of the poor preservation of specimens, many of which were road-killed animals.

Both species showed high mtDNA haplotype diversity and low nucleotide diversity (Mm, h = 0.913 ± 0.036, π = 0.0054 ± 0.00091; Hh, h = 0.924 ± 0.038, π = 0.006 ± 0.00086). Diversity estimates were significantly lower at nuclear genes (MC1R: Mm, h = 0.131 ± 0.082, π = 0.00022 ± 0.00014; Hh, h = 0.053 ± 0.04, π = 0.00009 ± 0.00009; BDNF: Mm, h = 0.434 ± 0.0095, π = 0.00156 ± 0.00035; Hh, h = 0.484 ± 0.065, π = 0.00092 ± 0.00015). Haplotype diversity at mitochondrial and nuclear genes is higher within North Africa in both species (Figs. 2, 3; Additional file 1: Fig. S2, 3).

Fig. 2

Statistical parsimony networks recovered in Hemorrhois hippocrepis: a cyt-b haplotype network; b mtDNA (cyt-b + 12S) haplotype network; c BDNF haplotype network; d MC1R haplotype network

Fig. 3

Statistical parsimony networks recovered in Malpolon monspessulanus: a cyt-b haplotype network; b mtDNA (cyt-b + 12S) haplotype network; c BDNF haplotype network; d MC1R haplotype network

For both H. hippocrepis and M. monspessulanus, the cyt-b and the concatenated mtDNA alignments were resolved into a single statistical parsimony network (Figs. 2a, b, 3a, b; Additional file 1: Fig. S2, 3). In the cyt-b network of both species we found one haplotype with high frequency, representing 42% and 73% of the samples for H. hippocrepis and M. monspessulanus, respectively. Singletons account for 6 out of 13 haplotypes for H. hippocrepis and 8 out of 10 haplotypes for M. monspessulanus (Figs. 2a, 3a). For the concatenated mtDNA alignments the most common haplotype found in H. hippocrepis represents 25% of the samples, whereas two equally common haplotype represent each 21% of M. monspessulanus samples. Singletons account for 9 out of 14 haplotypes for H. hippocrepis and 9 out of 13 haplotypes for M. monspessulanus (Figs. 2b, 3b). Also, BDNF and MC1R networks show a highly frequent haplotype, representing 67% of the samples for BDNF and 97% for MC1R in H. hippocrepis, and 73% of the samples in BDNF and 93% for MC1R in M. monspessulanus (Figs. 2c, d, 3c, d).

The phylogeographic pattern of both species shows exclusive mtDNA haplogroups in Iberia (Iberian haplogroup, blue colours; Figs. 2a, b, 3a, b; Additional file 1: Fig. S2, 3) and the High Atlas Mountains of Morocco (red–orange colours), whereas one haplogroup was shared across the Strait of Gibraltar (Ibero-Maghrebian haplogroup, green colours), with the most frequent haplotype being distributed in northern Morocco and southern Iberia. In H. hippocrepis one additional divergent haplogroup is found in eastern Algeria and Tunisia (yellow–brown colours). The number of exclusive mitochondrial haplotypes is slightly higher in North Africa than in the Iberian Peninsula for both H. hippocrepis and M. monspessulanus. Nuclear haplotype networks show a few haplotypes with a shallow phylogeographic structure. Exclusive nuclear haplotype, in both markers, are only present in North Africa and concentrated in mountain regions (Figs. 2, 3).

The historical demographic reconstructions based on mitochondrial and nuclear loci show a pattern of demographic expansion for both species during the last glacial period, that is particularly marked in M. monspessulanus with a sudden increase at ⁓ 40 kya (Fig. 4). For both species, a single demographic change obtained the highest posterior probability density.

Fig. 4

Extended Bayesian Skyline Plots of the historical demographic reconstruction: a Hemorrhois hippocrepis and b Malpolon monspessulanus

Predictive modelling

Species distribution models for H. hippocrepis and M. monspessulanus recovered AUC values of 0.96 ± 0.02 and 0.94 ± 0.02, respectively, indicating good performance of the models (Additional file 1: Table S5). For both snakes the Precipitation of Wettest Quarter (bio16) was the most informative variable, followed by the Mean Temperature of Wettest Quarter (bio8) (Additional file 1: Table S5). The present‐day habitat suitability reconstructions showed continuous high suitability across regions with Mediterranean climates, which fits well with the known species distributions (Fig. 5d, h). A lower fit is apparent for M. monspessulanus in the northernmost region of Iberia, characterised by oceanic climate, where the species is currently absent (see the species range in Fig. 1b), but the model shows some suitability (Fig. 5h). In both snakes, projections during LIG showed a lower habitat suitability in North Africa, with a fragmentation of suitable areas that were limited to coastal and mountain areas of Morocco, and with two isolated areas in western and eastern Maghreb (Fig. 5a, e). On the other hand, models predict an increase of suitable areas during the LGM, especially in the Maghreb (Fig. 5b, f), which was particularly evident for M. monspessulanus. Habitat suitability models at the LGM show similar patterns among the three projections in both species (MIROC, MPI and CCSM; Fig. 5b, f and Additional file 1: Fig. S4). Projections for Holocene recover areas of suitability for both snakes similar to present conditions in all models (Fig. 5c, g and Additional file 1: Fig. S4).

Fig. 5

Species distribution models (SDMs) of (ad) Hemorrhois hippocrepis and (eh) Malpolon monspessulanus for (a, e) Last Interglacial (LIG), (b, f) Last Glacial Maximum (LGM), (c, g) Holocene, and (d, h) Present conditions. For LGM and Holocene the SDM is based on the MIROC circulation model (see Additional file 1: Fig. S4 for SDMs based on the CC and MPI circulation models)

Fossil record

Compiled fossil data for M. monspessulanus show a continuous fossil record in Iberia since the Pliocene and throughout the Pleistocene (Fig. 6; Additional file 1: Tables S6, 7). The oldest possible fossil record for this species may be either from the Early Pliocene (4.9–5 Mya) of Puerto de la Cadena (Spain) (but this record is based on a few vertebrae) or the Late Pliocene (3.5–3.6 Mya), and 14 records from Iberia are referred to the Early Pleistocene (0.8–2.6 Mya). The fossil record of H. hippocrepis is poor with the oldest remains assigned to this species found in the Early Pleistocene (2.5 Mya) of Ahl al Oughlam (Morocco) and seven records from the Late Pleistocene from Iberia (0.127–0.01 Mya). All other records from Pliocene to Middle-Pleistocene (0.13–3.6 Mya) come from Iberia and are identified at the genus level.

Fig. 6

Fossil data available for Hemorrhois hippocrepis (grey circles) and Malpolon monspessulanus (dark circles) in the Pliocene, Early Pleistocene, Middle Pleistocene, and Late Pleistocene


The integrative approach used in this study, combining multi-locus genetic data, fossil records and climate suitability models, provided new insights into the evolutionary and demographic history of the snakes Hemorrhois hippocrepis and Malpolon monspessulanus, with implications on the biotic exchanges across the Strait of Gibraltar during Quaternary glaciations. Results support a pattern of demographic and range contraction in North Africa during interglacial epochs followed by glacial expansion, which is unusual across Mediterranean thermophilic organisms. Contrary to the previous hypothesis of glacial extinction for both species in Iberia [19], evidence supports a long history of the two species in the Iberian Peninsula with persistence during Pleistocene climatic oscillations. Finally, denser sampling in North Africa and longer gene fragments allowed identification of areas that during the Pleistocene served either as refugia or as source for the dispersal towards the southern Iberian Peninsula.

Glacial expansion for the two Mediterranean snakes during the Pleistocene

Historical demographic reconstructions indicate that both species underwent a demographic expansion during the Late Pleistocene (Fig. 4). According to our molecular calibration, the demographic expansion started around 40 kya in M. monspessulanus and sometime before in H. hippocrepis, a timeframe that coincides with the last glacial period. This finding is not consistent with demographic and phylogeographic patterns observed in Europe. In most temperate and Mediterranean species studied to date, range-wide patterns of genetic variation have been explained by phases of range and demographic contraction associated to unfavourable glacial climates, while expansion phases coincided with more suitable interglacial conditions [3, 56,57,58,59]. In the last decade, exceptions to this general ‘Expansion–Contraction’ model have been described, especially for coastal and insular species where demographic and range expansions have been associated to the glacial increase of climatically suitable coastal lowlands [60,61,62,63]. While far less case studies are available for North Africa, for mesic vertebrates such as H. hippocrepis and M. monspessulanus, that is, species adapted to arid conditions but still requiring some moisture, a model of expansion and contraction was proposed by Brito et al. (2014) [9], following the cycles of wet and arid phases in the area. According to Brito et al. (2014) [9], mesic species would expand their ranges under wet conditions associated with the Sahara Desert contraction, whereas, during unfavourable arid periods, mountain areas would act as refuge for these species by preserving mesic habitats. These predictions are in line with the results of phylogeographic, demographic and modelling analyses on H. hippocrepis and M. monspessulanus, suggesting a close association between the demographic and spatial expansion of these snakes with the expansion of Mediterranean habitats in North Africa during the last glacial period. Species distribution models indicate that in this region the amount of climatically suitable areas available to H. hippocrepis and M. monspessulanus was substantially higher during the last glacial period compared to the interglacials (Fig. 5 and Additional file 1: Fig. S4). The historical demographic reconstructions show that in both species a sudden demographic expansion coincided with favourable conditions in the Maghreb during the last glacial period, slightly earlier in H. hippocrepis than in M. monspessulanus. These snakes also show the highest haplotype diversity, including divergent and private haplotypes, in proximity of the High Atlas and Algerian mountains (Figs. 2, 3 and Additional file 1: Fig. S2–3) suggesting a key role of mountain areas of North Africa as long term refugia during Pleistocene climate oscillations. A remarkably similar pattern was found in a recent study on the North African snake Daboia mauritanica [47], that inferred range contractions and allopatric diversification during warmer interglacial periods and range expansion during colder glacial periods. Therefore, mounting evidence supports a demographic model for Mediterranean species in North Africa with allopatric isolation following retraction in mountain refugia during arid periods, and (demographic and spatial) expansion associated to an increase of Mediterranean habitats during the last glacial epoch. In this regard, it must be noted that the history of spatial and temporal climatic variability of North Africa is complex and still debated, and that a strict association of glacial stages with wet conditions, and of interglacial stages with arid conditions, is too simplistic [64,65,66]. Palaeoclimatic records show several wet periods in North Africa during the last glacial period, supporting a scenario of expansion of Mediterranean habitats, but with at least one period of aridification at ⁓ 60 kya [64, 65, 67]. On the other hand, the last interglacial might have been a short humid period [66], with arid conditions being prevalent just before it, from ⁓ 130 to ⁓ 170 kya [64]. Further phylogeographical and palaeoclimatic data are necessary to understand, on a finer time scale, the history of the North African biota during Pleistocene climatic oscillations.

Recent colonization or long-term persistence in the Iberian Peninsula?

A previous study based on mtDNA data recover a phylogeographic pattern of shallow divergence across the range of H. hippocrepis and M. monspessulanus, especially in the Iberian Peninsula [19]. This is unusual among Mediterranean taxa that usually show deep and complex phylogeographic structure in this region resulting from multiple isolated refugia during Pleistocene Ice Ages [68]. Such an unusual pattern of low genetic diversity in Iberia and limited differentiation between Maghreb and Iberian populations was explained as a very recent arrival (Late Pleistocene) in the Iberian Peninsula of these two snakes. According to this scenario, fossils of H. hippocrepis and M. monspessulanus from the Pliocene of Iberia derived from ancient populations of these species that underwent extinction during earlier glaciation events [19].

Based on the results of this study we can assess two additional hypotheses. First, the low mitochondrial diversity observed in the previous study in these two species could be explained by a selective process acting on the mitochondrial genome (i.e., ‘mitochondrial sweep’) [69]. A ‘mitochondrial sweep’ was suggested for two geckos to explain the shallow mitochondrial diversity observed in Europe as compared to higher nuclear polymorphisms observed across their range [70]. However, the nuclear data obtained in this study do not support this hypothesis for H. hippocrepis and M. monspessulanus, as we found very limited genetic diversity also in two nuclear genes. In this regard, while a low rate of evolution may explain the low diversity observed at the BDNF locus, the MC1R locus usually show high polymorphism and phylogenetic structure in many North African and European reptiles [70,71,72,73,74,75].

A second hypothesis is that one or more lineages endemic to Iberia were overlooked in the previous study, and these would provide evidence for the long-term persistence during Plio- Pleistocene glaciations of both species. Our results strongly support this hypothesis of a long history of the two species in Iberian Peninsula. For both species we unveiled two main haplogroups in Iberia, one shared with Northwest Africa (Ibero-Maghrebian haplogroup), and another one endemic to eastern and northern Iberia (Iberian haplogroup). The latter testify that at least part of the current species’ diversity in Europe evolved in situ and persisted during Pleistocene climatic oscillations. This is particularly evident in analyses based on longer gene fragments (Additional file 1: Fig. S2-3), suggesting that these endemic Iberian lineages may have been undetected in the previous study [19] because of a lack of phylogeographic resolution of the short gene fragment used. The long-term persistence in Europe of the two snakes is also in line with fossil data. Malpolon monspessulanus show an extensive fossil record (in many cases represented by several skeletal elements) in Iberia since the Pliocene (4.9–5 Mya) and the Early Pleistocene (1–2.6 Mya), and throughout Middle and Late Pleistocene [25, 76,77,78] (Fig. 6; Additional file 1: Table S7). For H. hippocrepis the fossil record is less rich and unfortunately older records from Iberia are only identified at the genus level. Finally, species distribution models indicate suitable areas in Iberia during last glacial period for both species, particularly in southeast Iberia, which are coincident with the occurrence of Iberian endemic haplogroups and with the majority of fossil remains (Figs. 2, 3, 5, 6), thus providing solid arguments against the hypothesis of a glacial extinction of both species in Europe [19].

While the endemic Iberian haplogroups indicate the persistence of these two species in Europe during Pleistocene glaciations, the Ibero-Maghrebian haplogroups found in southern Iberian and Northwest African populations of both species suggest recent gene flow across the Strait of Gibraltar. In line with the finding of Carranza et al. (2006) [19], the higher diversity of these Ibero-Maghrebian haplogroups is found in Morocco, supporting a recent colonization of the Iberian Peninsula from Northwest Africa. Trans-marine dispersal of these two species was likely facilitated by the reduced extent of the sea channel (despite a possibly faster marine flow) at the Strait of Gibraltar during sea-level low stands associated with glacial stages [11, 79]. Indeed, overseas dispersion across the Strait of Gibraltar have been inferred for several terrestrial organisms throughout the Pleistocene (see [14] for a review). According to our historic demographic reconstruction, during the last glacial stage a demographic expansion took place in both species that could have been associated with their colonization of southern Iberia across the Strait. Interestingly, in both species the observed phylogeographic pattern suggest that these Maghrebian lineages, after their arrival in Iberia, established a wide area of admixture with endemic Iberian lineages in the southern and south-eastern region.

To answer the main question on the biogeographic origin of the European populations of these two snakes, this study shows that the phylogeographic pattern of both species is more complex than previously thought and suggests that both long-term persistence in the east and north Iberia and recent colonization of south Iberia from North Africa have contributed to the current diversity of these two snakes within the Iberian Peninsula.

Insights into North African phylogeography

The inventory of phylogeographic patterns in North Africa is far from complete [7], however as more and more studies accumulate, some general biogeographic patterns start to emerge, suggesting that the genetic diversity of many Maghrebian taxa have been shaped by the same biogeographic events [80,81,82,83].

The major biogeographic break between western and eastern Maghreb, documented in different organisms (see [83] for a recent review), is found also in H. hippocrepis corresponding to the Kabylian region (central-eastern Algeria). The climatic suitability model for this species during the LIG shows a clear fragmentation between the western and eastern ranges, with low suitability across central Algeria (Fig. 5). During this period suitable areas in proximity of the Tell and Aurès mountains could have acted as an isolated eastern refuge for the species (Fig. 5).

Both H. hippocrepis and M. monspessulanus show two main lineages in the western Maghreb, one in northern regions from Morocco to central Algeria (Figs. 2, 3; haplogroups with green tones), and another one restricted to the High Atlas (Figs. 2, 3; haplogroups with orange and red tones). This latter lineage, endemic to southern mountains, was not detected by Carranza et al. (2006) [19] because of a lack of sampling in this region. This phylogeographic pattern suggests the occurrence of two isolated refugia in the western Maghreb during unfavourable period of the Pleistocene: one in northern Morocco, and another one in the High Atlas Mountains. The climatic suitability model for both species during the LIG support the occurrence of both refugia (Fig. 5).


In this study the use of denser sampling, longer gene fragments, multilocus analyses and species distribution modelling allowed the reconstruction of a more complex scenario for the evolutionary and demographic history of H. hippocrepis and M. monspessulanus across the Strait of Gibraltar compared to previous studies. The hypothesis of Pleistocene glacial extinction of both species in Europe is not supported based on increased genetic data and a detailed analysis of the available fossil record. Instead, results of phylogeographic, demographic and modelling analyses suggest that, while during the last glacial stage one lineage of both species re-colonised the southern Iberian Peninsula from North Africa, multiple populations of these two snakes survived in northern and eastern refugia in the Iberian Peninsula throughout Pleistocene Ice Ages. A model of glacial expansion (demographic and spatial) associated to an increase of Mediterranean habitats during the last glacial epoch, and of retraction in mountain refugia during arid phases, emerges as a general pattern for mesic vertebrates in North Africa [9, 47]. Finally, the phylogeographic pattern of H. hippocrepis conforms with a well-established biogeographic partition between western and eastern Maghreb [83].

Availability of data and materials

All genetic data in this study, either original or from previous studies, is available on GenBank (, with respective accession numbers listed in Additional file 1: Table S1. Presence data for Hemorrhois hippocrepis and Malpolon monspessulanus used as input in Maxent is available in Additional file 1: Table S3–4. Bioclimatic data used in Maxent was retrieved from Fossil Data for both species was retrieved from


  1. 1.

    Blondel J, Aronson J, Bodiou J-Y, Boeuf G. The Mediterranean region—biological diversity in space and time. Vasa. 2010.

  2. 2.

    Zachos FE, Habel JC. Biodiversity Hotspots. Springer. 2011. 123–208.

  3. 3.

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

    CAS  PubMed  Google Scholar 

  4. 4.

    Hewitt GM. Quaternary phylogeography: the roots of hybrid zones. Genetica. 2011;139(5):617–38.

    PubMed  Google Scholar 

  5. 5.

    Dobson M, Wright A. Faunal relationships and zoogeographical affinities of mammals in north-west Africa. J Biogeogr. 2000;27:417–24.

    Google Scholar 

  6. 6.

    Cosson JF, Hutterer R, Libois R, Sarà M, Taberlet P, Vogel P. Phylogeographical footprints of the Strait of Gibraltar and Quaternary climatic fluctuations in the western Mediterranean: a case study with the greater white-toothed shrew, Crocidura russula (Mammalia: Soricidae). Mol Ecol. 2005;14(4):1151–62.

    CAS  PubMed  Google Scholar 

  7. 7.

    Husemann M, Schmitt T, Zachos FE, Ulrich W, Habel JC. Palaearctic biogeography revisited: evidence for the existence of a North African refugium for Western Palaearctic biota. J Biogeogr. 2014;41:81–94.

    Google Scholar 

  8. 8.

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

    Google Scholar 

  9. 9.

    Brito JC, Godinho R, Martínez-Freiría F, Pleguezuelos JM, Rebelo H, Santos X, et al. Unravelling biodiversity, evolution and threats to conservation in the Sahara-Sahel. Biol Rev. 2014;89(1):215–31.

    PubMed  Google Scholar 

  10. 10.

    De Lattin G. No Grundriß der Zoogeographie. Jena: Verlag Gustav Fischer; 1967.

    Google Scholar 

  11. 11.

    Pleguezuelos JM, Fahd S, Carranza S. El papel del Estrecho de Gibraltar en la conformación de la actual fauna de anfibios y reptiles en el Mediterráneo Occidental. Boletín la Asoc Herpetológica Española. 2008;19:2–17.

    Google Scholar 

  12. 12.

    Jaramillo-Correa JP, Grivet D, Terrab A, Kurt Y, De-Lucas AI, Wahid N, et al. The Strait of Gibraltar as a major biogeographic barrier in Mediterranean conifers: a comparative phylogeographic survey. Mol Ecol. 2010;19(24):5452–68.

    CAS  PubMed  Google Scholar 

  13. 13.

    García-Vázquez D, Bilton DT, Alonso R, Benetti CJ, Garrido J, Valladares LF, et al. Reconstructing ancient Mediterranean crossroads in Deronectes diving beetles. J Biogeogr. 2016;43(8):1533–45.

    Google Scholar 

  14. 14.

    Mendes J, Harris DJ, Carranza S, Salvi D. Biogeographical crossroad across the Pillars of Hercules: evolutionary history of Psammodromus lizards in space and time. J Biogeogr. 2017;44(12):2877–90.

    Google Scholar 

  15. 15.

    Booth-Rea G, Ranero CR, Grevemeyer I. The Alboran volcanic-arc modulated the Messinian faunal exchange and salinity crisis. Sci Rep. 2018;8:13015.

    PubMed  PubMed Central  Google Scholar 

  16. 16.

    Duggen S, Hoernie K, Van den Bogaard P, Rüpke L, Morgan JP. Deep roots of the Messinian salinity crisis. Nature. 2003;422(6932):602–6.

    CAS  PubMed  Google Scholar 

  17. 17.

    Carranza S, Arnold EN, Wade E, Fahd S. Phylogeography of the false smooth snakes, Macroprotodon (Serpentes, Colubridae): mitochondrial DNA sequences show European populations arrived recently from Northwest Africa. Mol Phylogenet Evol. 2004;33(3):523–32.

    CAS  PubMed  Google Scholar 

  18. 18.

    Harris DJ, Batista V, Carretero MA, Ferrand N. Genetic variation in Tarentola mauritanica (Reptilia: Gekkonidae) across the Strait of Gibraltar derived from mitochondrial and nuclear DNA sequences. Amphibia-Reptilia. 1992;2004(25):451–9.

    Google Scholar 

  19. 19.

    Carranza S, Arnold EN, Pleguezuelos J. Phylogeny, biogeography, and evolution of two Mediterranean snakes, Malpolon monspessulanus and Hemorrhois hippocrepis (Squamata, Colubridae), using mtDNA sequences. Mol Phylogenet Evol. 2006;40(2):532–46.

    CAS  PubMed  Google Scholar 

  20. 20.

    Carranza S, Arnold EN. History of West Mediterranean newts, Pleurodeles (Amphibia: Salamandridae), inferred from old and recent DNA sequences. Syst Biodivers. 2004;1(3):327–37.

    Google Scholar 

  21. 21.

    Kaliontzopoulou A, Pinho C, Harris DJ, Carretero MA. When cryptic diversity blurs the picture: a cautionary tale from Iberian and North African Podarcis wall lizards. Biol J Linn Soc. 2011;103(4):779–800.

    Google Scholar 

  22. 22.

    Veríssimo J, Znari M, Stuckas H, Fritz U, Pereira P, Teixeira J, et al. Pleistocene diversification in Morocco and recent demographic expansion in the Mediterranean pond turtle Mauremys leprosa. Biol J Linn Soc. 2016;119(4):943–59.

    Google Scholar 

  23. 23.

    Agustí J, Santos-Cubedo A, Furió M, De Marfá R, Blain HA, Oms O, et al. The late Neogene-early Quaternary small vertebrate succession from the Almenara-Casablanca karst complex (Castellón, Eastern Spain): chronologic and paleoclimatic context. Quat Int. 2011;243(1):183–91.

    Google Scholar 

  24. 24.

    Stoetzel E, Denys C, Bailon S, El Hajraoui MA, Nespoulet R. Taphonomic analysis of amphibian and squamate remains from El Harhoura 2 (Rabat-Témara, Morocco): contributions to Palaeoecological and archaeological interpretations. Int J Osteoarchaeol. 2012;22(5):616–35.

    Google Scholar 

  25. 25.

    Bisbal-Chinesta JF, Blain HA. Long-term changes in composition and distribution patterns in the Iberian herpetofaunal communities since the latest Pleistocene. Quat Sci Rev. 2018;184:143–66.

    Google Scholar 

  26. 26.

    Sambrock J, Fritschi E, Maniatis T. Molecular cloning: a laboratory manual. New York: Cold Spring Harbor LaboratoryPress; 1989.

    Google Scholar 

  27. 27.

    Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.

    PubMed  PubMed Central  Google Scholar 

  28. 28.

    Palumbi SR. Nucleic acids II: the polymerase chain reaction. In: Hillis DM, Moritz C, Mable BK, editors. Molecular systematics. Massachusetts: Sinauer & Associates Inc; 1996. p. 205–47.

    Google Scholar 

  29. 29.

    Moritz C, Schneider CJ, Wake DB. Evolutionary Relationships within the Ensatina eschscholtzii complex confirm the ring species interpretation. Syst Biol. 1992;41(3):273–91.

    Google Scholar 

  30. 30.

    Kocher TD, Thomas WK, Meyer A, Edwards SV, Pääbo S, Villablanca FX, et al. Dynamics of mitochondrial DNA evolution in animals: amplification and sequencing with conserved primers. Proc Natl Acad Sci. 1989;86:6196–200.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Pinho C, Rocha S, Carvalho BM, Lopes S, Mourao S, Vallinoto M, et al. New primers for the amplification and sequencing of nuclear loci in a taxonomically wide set of reptiles and amphibians. Conserv Genet Resour. 2009;2(1):181–5.

    Google Scholar 

  32. 32.

    Vieites DR, Min M-S, Wake DB. Rapid diversification and dispersal during periods of global warming by plethodontid salamanders. Proc Natl Acad Sci U S A. 2007;104(50):19903–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Katoh K, Toh H. Recent developments in the MAFFT multiple sequence alignment program. Brief Bioinform. 2008;9(4):286–98.

    CAS  PubMed  Google Scholar 

  34. 34.

    Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68(4):978–89.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

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

    CAS  Google Scholar 

  36. 36.

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

    CAS  PubMed  Google Scholar 

  37. 37.

    Heled J, Drummond AJ. Bayesian inference of population size history from multiple loci. BMC Evol Biol. 2008;8:289.

    PubMed  PubMed Central  Google Scholar 

  38. 38.

    Ho SYW, Shapiro B. Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol Ecol Resour. 2011;11(3):423–34.

    PubMed  Google Scholar 

  39. 39.

    Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772–772.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Daza JM, Smith EN, Páez V, Parkinson CL. Complex evolution in the Neotropics: the origin and diversification of the widespread genus Leptodeira (Serpentes: Colubridae). Mol Phylogenet Evol. 2009;53(3):653–67.

    PubMed  Google Scholar 

  41. 41.

    Salvi D, Mendes J, Carranza S, Harris DJ. Evolution, biogeography and systematics of the western Palaearctic Zamenis ratsnakes. Zool Scr. 2018;47(4):441–61.

    Google Scholar 

  42. 42.

    Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7(1):214.

    PubMed  PubMed Central  Google Scholar 

  43. 43.

    Wickham H. Ggplot2: elegant graphics for data analysis. US: Springer; 2009.

    Google Scholar 

  44. 44.

    Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190:231–59.

    Google Scholar 

  45. 45.

    Hernandez PA, Graham CH, Master LL, Albert DL. The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography (Cop). 2006;29:773–85.

    Google Scholar 

  46. 46.

    Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.

    Google Scholar 

  47. 47.

    Martínez-Freiría F, Crochet PA, Fahd S, Geniez P, Brito JC, Velo-Antón G. Integrative phylogeographical and ecological analysis reveals multiple Pleistocene refugia for Mediterranean Daboia vipers in north-west Africa. Biol J Linn Soc. 2017;122(2):366–84.

    Google Scholar 

  48. 48.

    Silva-Rocha I, Salvi D, Sillero N, Mateo JÁ, Carretero MA. Snakes on the Balearic Islands: an invasion tale with implications for native biodiversity conservation. PLoS ONE. 2015;10(4):e0221026.

    Google Scholar 

  49. 49.

    Martínez del Mármol G, Harris DJ, Geniez P, de Pous P, Salvi D. Amphibians and Reptiles of Morocco. Chimaira. Frankfurt am Main; 2019.

  50. 50.

    GBIF. Accessed in 2019.

  51. 51.

    Bons J, Geniez P. Amphibians and Reptiles of Morocco (Including Western Sahara) Biogeographical Atlas. Asociacion Herptologica Espanola; 1996.

  52. 52.

    Kassambara A, Mundt F. Factoextra: Extract and Visualize the Results of Multivariate Data Analyses. 2017.

  53. 53.

    R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. Viena: R Foundation for Statistical Computing; 2019.

  54. 54.

    QGIS Association. QGIS Geographic Information System. 2021.

  55. 55.

    Böhme M, Ilg A. fosFARbase. 2003. Accessed March 2021.

  56. 56.

    Hewitt GM. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996;58(3):247–76.

    Google Scholar 

  57. 57.

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

    CAS  PubMed  Google Scholar 

  58. 58.

    Hewitt GM. Biodiversity: a climate for colonization. Heredity. 2004;92(1):1–2.

    PubMed  Google Scholar 

  59. 59.

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

    PubMed  Google Scholar 

  60. 60.

    Canestrelli D, Cimmaruta R, Nascetti G. Phylogeography and historical demography of the Italian treefrog, Hyla intermedia, reveals multiple refugia, population expansions and secondary contacts within peninsular Italy. Mol Ecol. 2007;16(22):4808–21.

    CAS  PubMed  Google Scholar 

  61. 61.

    Bisconti R, Canestrelli D, Colangelo P, Nascetti G. Multiple lines of evidence for demographic and range expansion of a temperate species (Hyla sarda) during the last glaciation. Mol Ecol. 2011;20:5313–27.

    PubMed  Google Scholar 

  62. 62.

    Salvi D, Schembri P, Sciberras A, Harris DJ. Evolutionary history of the Maltese wall lizard Podarcis filfolensis: insights on the ‘Expansion–Contraction’ model of Pleistocene biogeography. Mol Ecol. 2014;23:1167–87.

    PubMed  Google Scholar 

  63. 63.

    Senczuk G, Colangelo P, Harris DJ, Castiglia R, Litsi V, Canestrelli D, et al. Evolutionary and demographic correlates of Pleistocene coastline changes in the Sicilian wall lizard Podarcis wagleriana. J Biogeogr. 2019;46:224–37.

    Google Scholar 

  64. 64.

    Castañeda IS, Mulitza S, Schefuß E, Dos Santos RAL, Damsté JSS, Schouten S. Wet phases in the Sahara/Sahel region and human migration patterns in North Africa. Proc Natl Acad Sci U S A. 2009;106(48):20159–63.

    PubMed  PubMed Central  Google Scholar 

  65. 65.

    Drake NA, Breeze P, Parker A. Palaeoclimate in the Saharan and Arabian Deserts during the Middle Palaeolithic and the potential for hominin dispersals. Quat Int. 2018;2013(300):48–61.

    Google Scholar 

  66. 66.

    Braun K, Nehme C, Pickering R, Rogerson M, Scroxton N. A window into Africa’s past hydroclimates: the sisal_v1 database contribution. Quaternary. 2019;2(1):1–36.

    Google Scholar 

  67. 67.

    Hoffmann DL, Rogerson M, Spötl C, Luetscher M, Vance D, Osborne AH, et al. Timing and causes of North African wet phases during the last glacial period and implications for modern human migration. Sci Rep. 2016; 6.

  68. 68.

    Gomez A, Lunt DH. Refugia within refugia: patterns of phylogeographic concordance in Iberian Peninsula. In: Phylogeography of southern European refugia. 2007; 155–188.

  69. 69.

    Ballard JWO, Rand DM. The population biology of mitochondrial DNA and its phylogenetic implications. Annu Rev Ecol Evol Syst. 2005;36:621–42.

    Google Scholar 

  70. 70.

    Rato C, Carranza S, Perera A, Carretero MA, Harris DJ. Conflicting patterns of nucleotide diversity between mtDNA and nDNA in the Moorish gecko Tarentola mauritanica. Mol Phylogenet Evol. 2010;56(3):962–71.

    CAS  PubMed  Google Scholar 

  71. 71.

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

  72. 72.

    Harris DJ, Rosado D, Xavier R, Salvi D. New genetic lineages within Moroccan day geckos Quedenfeldtia (Sphaerodactylidae) revealed by mitochondrial and nuclear DNA sequence data. Amphibia-Reptilia. 2017;38(1):97–101.

    Google Scholar 

  73. 73.

    Rosado D, Rato C, Salvi D, Harris DJ. Evolutionary history of the Morocco lizard-fingered geckos of the Saurodactylus brosseti complex. Evol Biol. 2017;44(3):386–400.

    Google Scholar 

  74. 74.

    Salvi D, Pinho C, Harris DJ. Digging up the roots of an insular hotspot of genetic diversity: decoupled mito-nuclear histories in the evolution of the Corsican-Sardinian endemic lizard Podarcis tiliguerta. BMC Evol Biol. 2017;17(63):1–22.

    Google Scholar 

  75. 75.

    Sampaio FL, Harris DJ, Perera A, Salvi D. Phylogenetic and diversity patterns of Blanus worm lizards (Squamata: Amphisbaenia): Insights from mitochondrial and nuclear gene genealogies and species tree. J Zool Syst Evol Res. 2015;53(1):45–54.

    Google Scholar 

  76. 76.

    Bailon S. Amphibiens et reptiles du Pliocène et du Quaternaire de France et d’Espagne: mise en place et évolution des faunes. Université de Paris VII; 1991.

  77. 77.

    Blain H-A. Contribution de la paléoherpétofaune (Amphibia & Squamata) à la connaissance de l’évolution du climat et du paysage du Pliocène supérieur au Pléistocène moyen d’Espagne. Muséum National d’Histoire Naturelle de Paris; 2005.

  78. 78.

    Blain H-A. Contribution de la paléoherpétofaune (Amphibia & Squamata) à la connaissance de l’ evolution du climat et du paysage du Pliocène supérieur au Pléistocène moyen d’ Espagne. Treballs del Mus Geol Barcelona. 2009;16:39–170.

    Google Scholar 

  79. 79.

    Gibert J, Gibert L, Iglesias A. The Gibraltar strait: a Pleistocene door of Europe? Hum Evol. 2003;18(3–4):147–60.

    Google Scholar 

  80. 80.

    Barata M, Harris DJ, Castilho R. Comparative phylogeography of northwest African Natrix maura (Serpentes: Colubridae) inferred from mtDNA sequences. African Zool. 2008;43(1):1–7.

    Google Scholar 

  81. 81.

    Guiller A, Madec L. Historical biogeography of the land snail Cornu aspersum: a new scenario inferred from haplotype distribution in the Western Mediterranean basin. BMC Evol Biol. 2010;10(18):2–20.

    Google Scholar 

  82. 82.

    Nicolas V, Mataame A, Crochet P-A, Geniez P, Ohler A. Phylogeographic patterns in North African water frog Pelophylax saharicus (Anura : Ranidae ). J Zool Syst Evol Res. 2015; 53.

  83. 83.

    Salvi D, Perera A, Sampaio FL, Carranza S, Harris DJ. Underground cryptic speciation within the Maghreb: multilocus phylogeography sheds light on the diversification of the checkerboard worm lizard Trogonophis wiegmanni. Mol Phylogenet Evol. 2018;120:118–28.

    PubMed  Google Scholar 

Download references


We thank CIBIO colleagues for their help during sample collection. We thank Pedro Tarroso and Fernando Martinez-Freiria for their valuable suggestions on distribution modelling analyses, Fernando Seixas and Ana Couto for their help with R, Hugues-Alexandre Blain for advice on fossil data interpretation, and Salvador Carranza and Juan Pleguezuelos for fruitful discussion.


Fundação para a Ciência e a Tecnologia (FCT, Portugal): Doctoral Grant SFRH/BD/89820/2012 (L.M.) and IF contract 01627/2014 (D.J.H.); Program “Rita Levi Montalcini” (MIUR, Ministero dell'Istruzione dell'Università e della Ricerca) for the recruitment of young researchers at the University of L’Aquila (D.S).

Author information




DJH and DS conceived the ideas; LM, DJH, and DS collected the specimens; LM generated the data; LM and DS analysed the data; LM and DS led the writing which was revised by all authors. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Daniele Salvi.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The author declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Table S1: Details on samples’ location, gene fragments sequenced and GenBank accession numbers of sequences, for H. hippocrepis and M. monspessulanus. Table S2: Summary statistics and substitution models used for the Extended Bayesian Skyline Plot analyses. Table S3: Occurrence data for H. hippocrepis used as input for Maxent modelling analysis. Table S4: Occurrence data for M. monspessulanus used as input for Maxent modelling analysis. Table S5: Summary statistics of Maxent performance, each bioclimatic variable contribution (%) and output (logistic threshold for the “10th percentile”, “Equal sensitivity and specificity” and “Maximum sensitivity plus specificity”) for models of H. hippocrepis and M. monspessulanus. Table S6: Fossil data of Hemorrhois hippocrepis (source: fosFARbase, Bohme & Ilg, 2003). Mya: million years ago; epoch boundaries follow Gradstein et al. (2005). Table S7: Fossil data of Malpolon monspessulanus (source: fosFARbase, Bohme & Ilg, 2003). Mya: million years ago; epoch boundaries follow Gradstein et al. (2005). Fig. S1: Occurrence points used for Maxent climate modeling: (a) H. hippocrepis; (b) M. monspessulanus. Dashed lines represent species native distribution. Fig. S2: Haplotype network recovered in H. hippocrepis using the 700 bp cyt-b fragment. Fig. S3: Haplotype network recovered in M. monspessulanus using the 700 bp cyt-b fragment. Fig. S4: Species distribution models (SDMs) of H. hippocrepis and M. monspessulanus for Holocene and Last Glacial Maximum (LGM) conditions based on the MPI and CC circulation models.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Machado, L., Harris, D.J. & Salvi, D. Biogeographic and demographic history of the Mediterranean snakes Malpolon monspessulanus and Hemorrhois hippocrepis across the Strait of Gibraltar. BMC Ecol Evo 21, 210 (2021).

Download citation


  • Evolutionary history
  • Iberian Peninsula
  • Maghreb
  • Mediterranean Basin
  • North Africa
  • Pleistocene climatic oscillations