Skip to main content

Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones



The Iberian Peninsula is recognized as an important refugial area for species survival and diversification during the climatic cycles of the Quaternary. Recent phylogeographic studies have revealed Iberia as a complex of multiple refugia. However, most of these studies have focused either on species with narrow distributions within the region or species groups that, although widely distributed, generally have a genetic structure that relates to pre-Quaternary cladogenetic events. In this study we undertake a detailed phylogeographic analysis of the lizard species, Lacerta lepida, whose distribution encompasses the entire Iberian Peninsula. We attempt to identify refugial areas, recolonization routes, zones of secondary contact and date demographic events within this species.


Results support the existence of 6 evolutionary lineages (phylogroups) with a strong association between genetic variation and geography, suggesting a history of allopatric divergence in different refugia. Diversification within phylogroups is concordant with the onset of the Pleistocene climatic oscillations. The southern regions of several phylogroups show a high incidence of ancestral alleles in contrast with high incidence of recently derived alleles in northern regions. All phylogroups show signs of recent demographic and spatial expansions. We have further identified several zones of secondary contact, with divergent mitochondrial haplotypes occurring in narrow zones of sympatry.


The concordant patterns of spatial and demographic expansions detected within phylogroups, together with the high incidence of ancestral haplotypes in southern regions of several phylogroups, suggests a pattern of contraction of populations into southern refugia during adverse climatic conditions from which subsequent northern expansions occurred. This study supports the emergent pattern of multiple refugia within Iberia but adds to it by identifying a pattern of refugia coincident with the southern distribution limits of individual evolutionary lineages. These areas are important in terms of long-term species persistence and therefore important areas for conservation.


Evidence from phylogeographic studies suggests that southern Europe, including the peninsulas of Iberia, Italy and the Balkans and areas near the Caucasus and the Caspian Sea, have functioned as refugial areas for species survival during periods of adverse climatic conditions [1]. Recent studies also emphasize the important role that these refugial areas had in shaping the evolutionary history of species that have persisted within these regions for several ice ages. As suggested by early studies [24] the topographic complexity and geographic mosaic of habitats in southern refugial peninsulas have favoured the occurrence of multiple disjunct refugia, allowing the persistence of isolated populations within them during glacial periods. Within the Iberian Peninsula, complex species histories have been revealed for a variety of taxa, with some showing remarkable patterns of phylogeographic concordance [see [5] and references therein] involving deep genetic subdivisions, high haplotype richness and distinct hybrid zones. Not only has the Iberian Peninsula sourced the northern redistribution of species after ice ages, but it has also facilitated diversification through patterns of repeated population fragmentation, contraction, expansion and admixture. Phylogeographic analyses for the golden-striped salamander, Chioglossa lusitanica [68] and Schreiber's Lizard, Lacerta schreiberi [9, 10] provide good examples of the complexity that most likely typifies many species within this major peninsular glacial refugium.

However, even though the Iberian Peninsula is the best studied glacial refugium, most phylogeographic studies have focused on species that either have a narrow distribution within the region (e.g Chioglossa lusitanica [6], Lacerta schreiberi [9], Lissotriton boscai [11]) or involve species groups that, although distributed across the entire region, generally have a genetic structure that relates to older cladogenetic events (e.g Podarcis spp. [12, 13], Alytes spp [14], Oryctolagus cunniculus [15, 16]). In order to better understand the complex phylogeographic history of Iberian species, and the way they have responded to Pleistocene climatic oscillations, it is important to study species with distributions encompassing the entire Iberian Peninsula. For this purpose we have used the ocellated lizard, Lacerta lepida (Daudin, 1802), as a model to study the impact of Pleistocene climatic changes in generating and structuring intraspecific genetic diversity on this regional scale. Lacerta lepida is typically Mediterranean, with a distribution encompassing all the Iberian Peninsula, and with phylogeographic structure across the region [17]. Several mitochondrial lineages that appear to have non-overlapping geographic ranges were recently described, suggesting a history of allopatric differentiation in multiple refugia during the Plio-Pleistocene [17]. Here we assess the broader phylogeographic patterns within Lacerta lepida with the specific aims to i) clarify the distribution of mtDNA phylogroups; ii) identify refugial areas within these phylogroups during the glacial periods; iii) date the main demographic and evolutionary events within Lacerta lepida; and finally iv) identify secondary contact zones between the different phylogroups.

It is generally accepted that phylogeographic histories recovered using only mtDNA as a marker are constrained to reveal one genealogy that reflects the maternally inherited natural history of an organism. Relationships among phylogroups inferred through mtDNA might be discordant with inferences made based on nuclear genes (Harrison 1991; Avise 2000) and discordances have been illustrated in several studies [e.g. [1824]]. Particularly within the Iberian Peninsula, several studies have emphasized the importance of using different types of markers to fully recover the complex evolutionary and demographic scenarios that most likely characterize the species that have persisted there across the Quaternary. For example, in Lacerta schreiberi evidence for gene flow and ancestral introgression between apparently allopatric mtDNA lineages was only detected by the use of nuclear markers [10]. Here we use both mtDNA and nDNA derived genealogies, since their contrasting molecular and population properties (principally uniparental versus biparental mode of inheritance and contrasting population sizes) are mutually informative when there have been opportunities for secondary contact, gene flow and hybridization between diverging populations.



Lizards were captured under licence during the years 2005, 2006 and 2007, sampling most of the distribution of Lacerta lepida in Portugal, Spain, and France. Sampling intensity was concentrated in regions known to contain high genetic divergence within western and south-eastern part of Iberia [17]. Lizards were captured using tomahawk traps or by hand, and tissue samples were taken by clipping 1 cm of the tail tip that was subsequently preserved in 100% ethanol. After tissue sampling, animals were immediately released back into the wild in the place of capture. Geographic coordinates of sampling sites were recorded with a GPS.

DNA extraction, amplification and sequencing

Total genomic DNA was extracted from ethanol-preserved muscle tissue using a salt extraction protocol [25, 26]. A fragment of 627 base pairs (bp) of the mitochondrial DNA (mtDNA) cytochrome b (cytb) gene was amplified using the truncated version of primer L14841 [27] (CYTBF, 5'-CCA TCC AAC ATC TCA GCA TGA TGA AA-3') and the modified version of primer MV16 [28] (CYTBR, 5'- AAA TAG GAA GTA TCA CTC TGG TTT-3') to increase specificity for Lacerta lepida, using information from published sequences in GenBank. Polymerase chain reactions (PCRs) were performed in a total volume of 25 μl, containing 0.5 U of Taq polymerase (BIOTAQ™), 4 mM of MgCl2, 0.4 mM of each nucleotide (Bioline), 0.4 μM of each primer, 2 μl of 10 × NH4 reaction buffer (Bioline) and approximately 50 ng of DNA. PCR amplifications were conducted as follows: DNA was initially denaturated at 94°C for 3 min followed by 35 cycles of denaturation at 94°C for 45 s, annealing at 51°C for 45 s and extension at 72°C for 45 s, plus a final extension step at 72°C for 5 min. Negative controls (no DNA) were included for all amplifications. PCR products were visualized on a 2% agarose gel and purified by filtration through QIAquick® columns (Qiagen) following the manufacturer's recommendations. Purified PCR products were sequenced in both directions using the above primers, Taq polymerase BigDye Terminator v3.1™ (Applied Biosystems) and 30-90 ng of PCR product.

Intron 7 of the β-fibrinogen gene (β-Fibint7) has been successfully used as a nuclear marker in several vertebrate phylogeographic and phylogenetic studies [e.g. [13, 2932]]. Specifically it was recently employed for a phylogenetic study of the genus Lacerta [17] where it revealed sufficient variation within Lacerta lepida to be useful for phylogeographical inference. Initial β-fibint7 amplifications were performed using primers FIB-B17U (5'- GGA GAA AAC AGG ACA ATG ACA ATT CAC - 3') and FIB-B17L (5' - TCC CCA GTA GTA TCT GCC ATT AGG GTT - 3') [29] and the conditions described in Paulo et al. [17]. However, due to low amplification and sequencing success a nested PCR approach as suggested by Sequeira et al. [32] was subsequently adopted. A fragment of 788bp was first amplified from genomic DNA using primers FIB-B17U and FIB-B17L (PCRa). The product of this reaction (1 μl) was then used as a template for a subsequent PCR of 691bp (PCRb) using primer BFXF [32] and BFX8 [13]. Both amplifications were performed in a total volume of 25 μl, and included reagents in the same concentrations as those specified for cytb gene fragment. PCR cycle conditions were the same as described for cytb fragment but the primer annealing temperatures were 55°C and 56°C, for PCRa and PCRb respectively. Negative controls (no DNA) were included for all amplifications. Purified PCRb products were then sequenced with primers BFBX and BFX8 using identical sequencing conditions as for the mtDNA cytb sequencing. All PCRs and sequencing reactions were performed in a MJ Research thermocycler (PTC-240 DNA Engine Tetrad 2) and sequences were obtained using an ABI 3700 capillary sequencer.

Sequence analysis

DNA sequences were aligned by eye using BioEdit Sequence Alignment Editor 7.01 [33]. β-fibint7 alleles of heterozygous individuals were inferred using Phase version 2.1 [34, 35]. Ten replicate Phase runs were conducted using 1000 burnin steps and 1000 iterations. The Phase probability threshold was first set to 0.90, but not all genotypes were resolved at this threshold. Garrick et al. [36], suggest that lowering the Phase threshold to 0.60 reduces the number of unresolved haplotypes with little or no increase in the number of false positives. In contrast, omitting unresolved haplotypes may have a substantial impact on downstream analyses, such as the estimation of Tajima's D and Fu's Fs statistics, as unresolved haplotypes usually represent rare alleles [36]. We therefore reduced the Phase threshold to 0.60 to maximise recovery of haplotypes. Several tests implemented in the software RDP3 [Recombination Detection Program, [37]] were used to detect evidence for recombination in β-fibint7: RDP [38], GENECONV [39], Maximum Chi Square [40, 41], Chimaera [40] and Sister Scanning [42]. Due to the small size of the fragment used in the analyses (315 bp), the window size for the recombination detection methods was set to 20bp.

Phylogenetic analysis and haplotype network construction

Phylogenetic relationships among cytb haplotypes were inferred using Bayesian inference as implemented in Mrbayes v3.1.2 [43, 44]. Data was divided into two partitions: (1st+2nd) and (3rd) codon positions. The most appropriate substitution model for each partition was selected using Modeltest v3.7 [45] and the Akaike Information Criterion (AIC) [46]. The model selected for both partitions was the same (GTR + I) but with a different proportion of invariable sites (pinv(1st+2nd) = 0.79 and pinv(3rd) = 0.09). The different partitions were allowed to evolve at different rates, with unlinked topology and unlinked parameters for the nucleotide substitution models. Four Markov chains were run for 10 million generations. To avoid local optima we used two independent runs, and to improve swapping of states between heated and cold chains the heating parameter was decreased to 0.05. Trees were sampled every 1000 generations and the first 100 trees were discarded as burn-in. Posterior probabilities were obtained from the 50% majority-rule consensus tree of the retained trees.

Intraspecific gene genealogies were inferred using the median-joining (MJ) [47] and the statistical parsimony (SP) [48] network construction approaches. The MJ network was computed with the program Network 4.5.0 ( and the SP network was inferred using the program Tcs 1.21 [49]. For the MJ approach the parameter ε was set to 0, preventing less parsimonious pathways from being included in the analysis. The SP network was inferred with a parsimony confidence limit of 95%, allowing the inclusion of less parsimonious alternatives that fall within this confidence limit. Ambiguities within networks were resolved where possible following the criteria of Crandall & Templeton [50]. The phylogenetic and β-Fibint7 network analyses were rooted using gene sequences of the African sister species Lacerta pater (Lataste, 1880) [17, 51] (Genebank accession numbers: AF: 378967 and EU: 365413 for cytb and β-Fibint7 genes respectively).

Estimation of divergence times

Divergence times within and between phylogroups were estimated using Beast version 1.4.8 [52] and the cytb dataset. Beast performs Bayesian statistical inferences of parameters, such as divergence times, by using Markov Chain Monte Carlo (MCMC) as a framework. Input files were generated with Beauti version 1.4.8 [53]. The nucleotide substitution model and its parameter values were selected according to the results of Modeltest version 3.7 [45], with upper and lower bounds around the values defined as 120% and 80% respectively [54]. According to the Bayesian Information Criterion (BIC) and the hierarchical Likelihood Ratio Tests (hLRT's), the model of nucleotide substitution identified as the best fit to the data is the HKY model [55] with a gamma distribution (Γ) for substitution rates across sites (shape parameter, α = 0.2889) and no category of invariable sites. An uncorrelated lognormal relaxed molecular clock was used with mean mutation rate of 0.01 mutations/site/million years and standard deviation of 0.0015, assuming a normal distribution, as prior information. The mutation rate used was based on the provisional calibration of a reptile-specific molecular clock presented in Paulo et al. [9] who used published cytochrome b data sets from two studies of Canary island reptiles (Gallotia spp.) [56, 57] calibrated with the geological age of El Hierro. The standard deviation included in our study encompasses a faster mutation rate of 0.0115 and a slower mutation rate of 0.0085 mutations/site/million years. This takes into account potential underestimation of the mutation rate due to violation of the assumption of immediate island colonization in the clock calibration of Paulo et al. [9], or overestimation due to the longer generation time and larger body size of Lacerta spp. when compared to Gallotia spp. [58]. Two runs were executed for 106 generations, sampling every 500 generations and discarding the first 10% as burn in. Results of the two runs were displayed and combined in Tracer v1.4 [59] to check for stationarity and ensure that ESSs were above 200. For all analyses one sequence of Lacerta pater (Genebank accession number: AF378967) was included as an outgroup.

In a second approach to estimate divergence times within phylogroups the method of Saillard et al. [60] was employed where each extant haplotype descending from the most recent common ancestor (MRCA) represents a time interval between the present and the MRCA. Average distances from the MRCA within each phylogroup are calculated from the number of mutational steps separating each haplotype from the MRCA. The absolute timing of divergence is then calculated by multiplying the observed values by the average mutational changes per lineage per million years (Myr). Three mutation rates were used: 0.01 mutations/site/million years, representing the average mutation rate; a faster mutation rate of 0.0115, and a slower mutation rate of 0.0085.

Neutrality tests and demographic analyses

In order to detect departures from a constant population size under the neutral model, Tajima's D [61], Ramos-Onsins & Rozas' R 2 [62] and Fu's Fs [63] tests were applied to both mtDNA and nDNA datasets. Both R 2 and Fs statistics have been shown to be the best statistical tests to detect population growth (R 2 has been suggested to behave better for small sample sizes whereas Fs is better for bigger ones) [62]. Population expansions have also been shown to leave particular signatures in the distribution of pairwise sequence differences [64, 65]. We capitalized upon this by employing statistics based on the mismatch distribution to test for demographic expansions. The observed distribution of pairwise differences between haplotypes within phylogroups was compared with the expected results under a sudden-demographic and a spatial-demographic expansion model. Statistically significant differences between observed and simulated expected distributions were evaluated with the sum of the square deviations (SSD) and Harpending's raggedness index (hg) [66, 67]. For the β-Fibint7 gene we used two datasets: the first including only haplotypes inferred by Phase with probability higher than 0.90, and the second including all haplotypes inferred with the less stringent Phase probability threshold of 0.60. Tests were performed with Arlequin version 3.11 [68] for Tajima's D, Fu's Fs, SDD and hg, and with DNAsp version 4.50 [69] for R 2 and expected values for the mismatch distribution. The demographic history of each phylogroup was also inferred using a coalescent-based approach. The model used to infer past population dynamics was the Bayesian Skyline Plot (BSP, [70]) as implemented in Beast version 1.4.8 [52]. For each phylogroup we carried out two independent runs of 10 million generations each. Trees and parameters were sampled every 1000 iterations, with a burn in of 10%. Results of each run were visualized using Tracer v1.4 [59] to ensure stationarity and convergence had been reached, and that effective sample sizes (ESS) were higher than 200.

Geographical distribution of alleles and inference of refugial areas

Using predictions from coalescent theory (haplotypes at the tips of a tree are younger than interior haplotypes to which they are connected) ancestral and derived haplotypes within each phylogroup were identified, thus obtaining a temporal framework for inferring haplotype origin within phylogroups. The null hypothesis of random geographic distribution of haplotypes was tested statistically using Geodis version 2.5 [71] with significance testing by permuting the data 106 times. When non-random associations of haplotypes with geography were detected, the geographic distribution of ancestral versus derived haplotypes (interior and tip haplotypes, respectively) was further explored to identify possible refugial areas and the directionality of previously detected demographic and spatial expansions. Coalescent theory predicts refugial areas to exhibit a high frequency of ancestral haplotypes [72], and therefore the average number of mutations of a haplotype from the MRCA is expected to be significantly lower in refugial areas than would be expected by chance. In contrast, recently colonized areas are predicted to show a high incidence of derived haplotypes [72], and therefore the average number of mutations of a haplotype from the MRCA is expected to be significantly higher than would be expected by chance. To test the hypothesis that refugial areas are located in the south of a phylogroup's distribution, with the north having been more recently colonized, we estimated the average distances of alleles from the MRCA within both the north and south of the geographic range of a phylogroup. We then tested whether these values were significantly high or low by randomizing the geographic states of alleles to generate a null distribution of mean haplotype distance from the MRCA within both the south and north of a phylogroup's distribution.


A total of 353 lizards were sampled from 129 different sites across the distribution of Lacerta lepida. Sampling sites and number of samples per site are shown in Figure 1 and Additional file 1, Table S1, respectively. A total of 321 cytb and 104 β-Fibint7 sequences were obtained.

Figure 1
figure 1

Map of the Iberian Peninsula showing the current distribution of Lacerta lepida and sampled localities. a) Distribution of the recognized continental subspecies of Lacerta lepida (L. lepida. iberica, L. lepida lepida and L. lepida nevadensis). b) Sampling localities. Numbers are as in Additional file 1, Table S1. Shaded areas denote altitude gradients, with darker areas representing higher altitudes.

Sequence variability

All cytb sequences represented uninterrupted open reading frames, with no gaps or premature stop codons, suggesting they are functional mitochondrial DNA copies. The chromatograms of 9 sequences were polymorphic at several nucleotide positions. These samples were from individuals collected at sampling sites 59 (1 sample), 79 (2 samples), 80 (1 sample), 84 (1 sample), 86 (2 samples), 88 (1 sample) and 89 (1 sample). As it was not possible to correctly identify the corresponding mitochondrial sequence of each sample these sequences were eliminated from further analysis. One hundred and forty five (145) unique haplotypes were obtained from the 312 sequences analysed. Of a total of 627 sites sequenced, 168 were variable, from which 129 were parsimony informative. Pairwise genetic distances (uncorrected p-values) between haplotypes ranged from 0.16% to 13.2% (Table 1). β-fibint7 sequences were trimmed to 315 bp in order to eliminate gaps within the sequences. From the 315 bp, 15 sites were variable of which 6 were parsimony informative. The haplotypes inferred by Phase for heterozygous individuals were consistent across runs. Twenty unique haplotypes (B1 to B20) were identified among the 208 alleles analysed and recombination was not detected by any of the tests applied. When using a Phase probability threshold of 0.90, 8 genotypes were unresolved (Additional file 1, Table S1), and this included the loss of only a single haplotype (B7) that was additionally recovered with a probability threshold of 0.60.

Table 1 Average pairwise genetic distances between Lacerta lepid a mitochondrial phylogroups.

Phylogenetic and network relationships among haplotypes

The results of the phylogenetic and network analyses both support the existence of two main phylogroups: phylogroup N and phylogroup L (Figure 2). The geographic distribution of phylogroup N is coincident with the Betic Mountains in south-eastern Spain and is associated with the subspecies Lacerta lepida nevadensis (Buscholz, 1963), while phylogroup L occupies the remaining area of the species distribution. Within phylogroup L there are four geographically structured monophyletic groups (L1, L2, L3 and L5), albeit with only moderate Bayesian support for the monophyly of L5 (Figure 2). The genealogical relationships among cytb haplotypes inferred by the two approaches for network construction (MJ and SP) are highly congruent. At the 95% confidence limit TCS produced three unconnected networks (results not shown). When the confidence limit was reduced (to 92% to include group L2 and to 65 mutational steps for group N) the three SP networks were connected in one single network (Figure 3), with relationships identical to those from the MJ analysis, with 12 loops, of which 7 were easily resolved by applying the criteria of Crandall and Templeton [50].

Figure 2
figure 2

Phylogenetic relationships among cytochrome b sequences (627 bp). Fifty percent majority-rule consensus phylogram from the Bayesian inference analysis. Numbers above branches represent posterior probabilities. Haplotype numbers are the same as in Figure 3 and as in Additional file, Table S1. Colours represent the geographic distribution of haplotypes as in Figure 5.

Figure 3
figure 3

Statistical Parsimony network of Lacerta lepida cytochrome b haplotypes from 312 samples. Open circles with no numbers represent unsampled or extinct haplotypes. L1, L2, L3, L4, L5 and N represent different mitochondrial phylogroups. The ancestral haplotype within each phylogroup is marked with an asterisk. Phylogroup N connects to the main network through 65 mutations, represented by an interrupted line.

Within phylogroup L, five geographically distinct groups of haplotypes can be identified (Figure 3 and Figure 4), which include the four mitochondrial phylogroups (L1-L4) identified by Paulo et al. [17] and a new lineage of haplotypes (L5) possessing a geographically distinct distribution from other haplotypes within phylogroup L4. Average genetic distances (uncorrected p distances) between these phylogroups range from 1.1% (between phylogroup L4 and L5) to 3.28% (between phylogroups L2 and L5; and L2 and L1). Phylogroup L1 is distributed mainly across the Central Mountain system between the Douro and Tagus river basins in Spain. Phylogroup L2 is distributed in southern Portugal, occupying the regions of Algarve and Baixo Alentejo. This phylogroup is clustered together in the network with phylogroup L3 (Figure 3) forming a monophyletic group in the phylogeny (Figure 2). L3 is distributed several hundred kilometres (300 to 450 km) to the north of L2 occupying the north-western corner of Iberia, mainly the regions to the north of the Douro River in Portugal and the regions of Asturias and Galicia in Spain. The geographic region between phylogroups L3 and L2 is occupied by two other phylogroups, L4 and L5. Phylogroup L5 is restricted to central Portugal, occupying the region between the Tagus and Douro rivers. L4 has the broadest distribution, occupying the remaining areas of southern Portugal and Spain passing through the Ebro valley to reach the Atlantic and Mediterranean coasts of France.

Figure 4
figure 4

Distribution of Lacerta lepida mitochondrial phylogroups based on 627 bp of the cytochrome b gene. Colours are the same as in Figures 2 and 3. Filled red circles represent populations where divergent haplotypes from two or more phylogroups were detected in sympatry. Numbers correspond to sampling localities as in Figure 1 and Additional file 1, Table S1.

The root of the network is located along the branch that connects the very divergent lineages L and N, allowing for the inference of the most recent common ancestor (MRCA) for each phylogroup (Figure 3). Although this identification is straightforward for clade N (haplotype 126), the networks reveal two probable ancestral haplotypes within clade L (haplotypes 127 and 104), which are connected to haplotype 126 through a loop. SP and MJ networks constructed with 0-fold degenerate sites only, thus reducing homoplasy within the data set (data not shown) result in the collapse of this reticulation, and haplotype 126 (lineage N) connects unambiguously to haplotype 127 (lineage L), supporting 127 as the ancestral haplotype within phylogroup L. Divergent mitochondrial haplotypes from different phylogroups occurring in sympatry were detected in several populations. The admixed populations were: populations 6, 23, 24 and 25 (phylogroups L2 and L4); populations 15, 17, 18, 19, 20, 21 22 and 96 (phylogroups L4 and L5); population 33 (phylogroups L4 and N); population 44 (phylogroups L1 and L4) and population 56, 58, 78 and 82 (phylogroups L1, L3 and L5) (Additional file 1, Table S1).

The relationships among β-fibint7 haplotypes inferred by the two network construction approaches (MJ and SP) were identical, resulting in a single network with 4 unresolved reticulations (Figure 5). Haplotype B15 is inferred to be the ancestral haplotype as it connects unambiguously to the outgroup. Haplotype B15 is restricted to the south-west of the species distribution (sampling sites 11, 56, 72, 94 and 112). Haplotype B1 is the most common haplotype and has the widest distribution within the group, occurring in 83% of samples. This haplotype is connected to several low frequency haplotypes, generating a star-like genealogy, suggestive of a possible past range expansion and for which signatures of expansion were formally detected by mismatch distribution analysis and neutrality tests (Table 2). Within the 50 southernmost samples, 16 alleles are found, representing 80% of the nuclear allele diversity. From those 16 alleles, 8 are restricted to the southern area.

Table 2 Results from mismatch distribution and neutrality tests for cytb mtDNA phylogroups and for the β-fibint 7 nuclear gene.
Figure 5
figure 5

Statistical Parsimony network of Lacerta lepida β-Fibrinogen intron 7 alleles from 104 samples. Open circles with no numbers represent unsampled or extinct alleles, and the filled black circle represents the outgroup (Lacerta pater). Pie chart shading represents the proportion of each allele found within each mitochondrial phylogroup. Colours in pie charts are the same as those used to represent mitochondrial phylogroups in Figure 3. Dashed lines represent ambiguities in the network.

The nuclear data failed to recover the phylogroups detected by the mtDNA dataset. Nevertheless, when the nuclear dataset is grouped according to the mtDNA phylogroups previously identified, structure in the distribution of alleles can be detected. For this analysis each of the mtDNA phylogroups was considered as a geographic region and Geodis was used to test for geographical structure amongst the nuclear genetic variation. Nuclear haplotypes are shared among some of the mtDNA phylogroups (with the exception of phylogroup N): all haplotypes from L3 (B7 and B20) occur either in L1 (B7) or in L5 (B20); all haplotypes from L2 (B4, B6, B13 and B14) occur in L4 (with B13 also occurring in L5); all haplotypes from L5 (B5, B13 and B20) occur either in L4 (B5 and B13), either in L3 (B20) or in L2 (B13) and finally haplotypes from L1 (B7 and B17) occur either in L3 (B7) or in L4 (B17). Additionally all phylogroups share haplotype B1 (the most common and one of the most ancestral haplotypes) and all phylogroups apart from L3 and L5 also share haplotype B15 (the most ancestral haplotype). It is important to note that only geographically close phylogroups share nuclear haplotypes. It is also possible to identify private haplotypes within the geography of several mtDNA phylogroups: B2 was only detected in phylogroup L1; haplotypes B3, B9, B10, B14, B16 and B18 in L4 and haplotypes B8, B11, B12 and B19 in phylogroup N.

Divergence times

Mean ages and 95% highest posterior density (HPD) of mtDNA phylogroups are shown in Table 3. Divergence within the group is estimated to have started approximately 9.4 million years ago (Ma) (5.58-13.66) in the mid-late Miocene, corresponding to the cladogenetic event between phylogroups N and L. Although divergence within phylogroup L is estimated to have started in the Plio-Pleistocene (1.96 Ma; 1.13-2.91) the current diversity within each phylogroup seems to have Pleistocene origins, with diversification times for all phylogroups estimated to be younger than 1.0 Ma. The oldest split within group L relates to the divergence of the monophyletic lineage composed of phylogroups L2 and L3 from the remaining phylogroups.

Table 3 Divergence time estimates in million years (Ma) from the most recent common ancestor (mrca) of all group members from each Lacerta lepida phylogroup estimated using the method of Saillard et al. (2000) using three different mutation rates; and divergence time estimates for the main nodes recovered in the phylogenetic analysis using a mean mutation rate of 2% in Beast.

Neutrality tests and demographic analyses

Significant deviations from neutrality that could reflect past population expansion events were detected for almost all phylogroups with Tajima's D, with the exception of phylogroup L3. When more powerful statistics were applied, a similar pattern of significant deviations from neutrality were observed, with the exception of phylogroups L1 and L3 (Table 2). The distributions of pairwise differences within each phylogroup were also found to be consistent with sudden-expansion and spatial-expansion models (as seen by the SDD and hg p values in Table 2), with signatures of population growth being exhibited by the bell shaped mismatch distributions (Figure 6). Historical demographic reconstructions (BSPs) also show a trend of population growth in all phylogroups (Figure 7).

Figure 6
figure 6

Mismatch distribution of mtDNA haplotypes for each of the 6 Lacerta lepida phylogroups. The expected frequency is based on a population growth-decline model, determined using DnaSP v4.50 [69] and is represented by a continuous line. The observed frequency is represented by a dotted line.

Figure 7
figure 7

Bayesian skyline plots showing the historical demographic trends for each Lacerta lepida mitochondrial phylogroup detected using cytochrome b sequences. Along the y-axis is the expressed population size estimated in units of Neμ (Ne: effective population size, μ: mutation rate per haplotype per generation). The y-axis is in a log-scale. Solid lines represent median estimates and shaded areas represent confidence intervals.

Geographical distribution of alleles and inference of refugial areas

Within phylogroups, statistically significant associations between genetic variation and geographic distribution were detected for L1, L3, L5 and N (Table 2). More detailed geographic sampling within phylogroups L3 and L5 allowed for the testing of hypotheses of southern refugial origin and northern expansion. The average number of mutations of haplotypes from the MRCA in southern regions of phylogroups L3 and L5 was significantly lower than that expected by chance (2.25 and 3.03 for L3 and L5, respectively) (Figure 8). In contrast, northern haplotypes exhibit a significantly higher average number of mutations from the MRCA than expected by chance (3.85 and 4.04 for L3 and L5, respectively) (Figure 8). These results reveal significant associations of ancestral haplotypes with southern ranges and derived haplotypes with northern regions for these two phylogroups.

Figure 8
figure 8

Observed and null model expected distances of haplotypes from the MRCA in southern and northern regions of phylogroups L3 and L5. Null model distributions were generated by randomising the geographic states of alleles (see text for details). Arrows indicate observed distances. For both phylogroups the average distance of southern haplotypes from the MRCA is significantly less than would be expected by chance, and the average distance of northern haplotypes from the MRCA is significantly greater than would be expected by chance. L3 southern populations: 78, 81, 82, 83, 85, 93 and 94; L3 northern populations: 87, 89, 90, 91 and 92; L5 southern: 15, 17, 18, 19, 20, 21, 22, 53, 95, 96, 97, 98, 99 and 100; L5 northern populations: 56, 58, 78, 79, 82, 101, 102, 103, 104, 105, 106 and 108.


Detailed analysis of the distribution of mitochondrial genetic variation within Lacerta lepida across the Iberian Peninsula revealed a complex phylogeographic history for the species. The cytb genealogy clearly defines 6 geographic phylogroups with generally non-overlapping geographic distributions. The strong association of mtDNA genetic variation with geography suggests a history of allopatric divergence in different refugia within the Iberian Peninsula, a pattern that has been described for several taxa within the region [see [5], and references therein]. Although this pattern of differentiation of distinct evolutionary units in allopatry was less evident from the analysis of the nuclear data, the distribution of nuclear haplotypes is not in conflict with the mtDNA phylogroups. Failure of the nuclear gene genealogy to reveal concordant genetic structure with the mitochondrial genealogy can be expected if one takes into account the fact that nuclear genes take on average four times longer to reach monophyly than mitochondrial ones [73]. In fact, most of the intraspecific differentiation within Lacerta lepida is of relatively recent origin, with the majority of phylogroups (L1 to L5) estimated to have diversified within the Pleistocene, increasing the probability of mitochondrial lineages not being reciprocally monophyletic for nuclear markers. Nevertheless, phylogroup N is estimated to have diverged from the remaining phylogroups during the Miocene (5.6-13.7 Ma) representing a much older cladogenetic event within the group. This older split therefore allows for more time for lineage sorting at the nuclear level. This is evident in the composition of nuclear haplotypes of phylogroup N, which are almost all private with the exception of the ancestral haplotypes (which are shared among almost all phylogroups). Thus, although mtDNA lineages have not reached monophyly at the nuclear level, some level of nuclear differentiation between mtDNA lineages is detected by the existence of private alleles, and this is most pronounced for phylogroups that represent older cladogenetic events (L1, L4 and N). Incomplete lineage sorting seems therefore a plausible explanation for the discrepancies in mitochondrial and nuclear gene genealogies. Recent phylogeographic and phylogenetic studies focusing on lizards in the Iberian Peninsula reveal similar discrepant patterns between mtDNA and nuclear genealogies [10, 13]. In the case of Podarcis wall lizards Pinho et al. [13] showed that, despite a complete lack of monophyly at the nuclear level, most of the species show historical reproductive isolation and that the lack of monophyly is mainly a result of shared ancestral polymorphism.

Although it is plausible that incomplete lineage sorting is responsible for the failure of nuclear genealogy to recover the mitochondrial phylogroups in Lacerta lepida, current gene flow, mainly male-biased, at zones of secondary contact can also be invoked to explain the patterns observed. Evidence for this comes from the observation that geographically close phylogroups share more derived haplotypes. For example allele B20 occurs only in the very divergent phylogroups L3 and L5 near the zone of contact, but it was not detected in the ancestral phylogroup L4, suggesting that nuclear gene flow may be occurring between the L3 and L5 mtDNA lineages. The same is true for allele B4 which, although more widespread within phylogroup L4, is also found in individuals of L2 near the zone of contact, but was not detected in phylogroup L3, which is phylogenetically closer to phylogroup L2.

Historical biogeography of Lacerta lepida

Divergence within Lacerta lepida is estimated to have started in the Miocene, approximately 9.4 Mya (5.58-13.66), with divergence of phylogroup N from phylogroup L. Within phylogroup L estimated divergence times are much younger and inferred to have been initiated in the Plio-Pleistocene, approximately 1.96 Mya (1.13-2.91 Mya). Interestingly, although divergences between the mitochondrial phylogroups have a Plio-Pleistocene (within group L) or a late Miocene (for group N) origin, haplotype diversities within phylogroups indicate a strong influence of later Pleistocene events, with diversification within phylogroups starting between 0.92 and 0.45 Ma. The importance of the Pleistocene climatic oscillations in promoting species differentiation in the Iberian Peninsula has been emphasized by previous studies [see [5], for a recent review] and this clearly also seems to be the case for Lacerta lepida.

Phylogroup N

Phylogroup N is distributed across the Betic Mountain range in south-western Spain and its distribution approximately coincides with the described subspecies Lacerta lepida nevadensis. The existence of clear significant morphological differentiation between L. l. nevadensis and the remaining subspecies of Lacerta lepida [7476], together with high levels of allozyme differentiation [76] and mitochondrial genetic differentiation (our study and [17]), suggests that the two lineages (N and L) might in fact reflect two different species. Paulo et al. [17] have inferred that divergence between phylogroups L and N were initiated by overseas dispersal between what was then the Iberian mainland and the Betic Massif that at that time existed as an island between Iberia and North Africa. Under this scenario subsequent contact between the phylogroups would have been initiated after the merging of the Betic Massif with Iberian mainland, due to the closing of the Betic corridor 7.8-7.6 Ma [see [77, 78] for a detailed explanation of the kinematics of the western Mediterranean basin]. The Betic Mountain range is a region of high endemism for plants and animals, and its importance as a refugium for other taxa has been highlighted before [see [5]]. Interestingly, some of these taxa show similar divergence times and distribution as phylogroup N (e.g. Salamandra salamandra longirostris [79] and Alytes dickhilleni [80]).

Phylogroups L2 and L3

The monophyletic group composed of phylogroups L2 and L3 is estimated to have started diverging from the remaining phylogroups in the early Pleistocene, approximately 1.5 Mya (0.82-2.27). Interestingly these two phylogroups have a disjunct distribution, with phylogroup L2 occupying the south of Portugal whereas L3 occupies the north-western parts of the Iberian Peninsula. The intervening region between phylogroups L2 and L3 is occupied by phylogroup L5. A vicariant event during the Middle Pleistocene (0.82-2.27 Mya) triggering divergence between the L2-L3 lineage and the remaining populations of Lacerta lepida seems probable. It is noteworthy that most phylogeographic studies within Iberia reveal similar phylogenetic breaks associated with the same period (e.g. Chioglossa lusitanica [6, 7], Oryctolagus cuniculus [81], Lacerta schreiberi [9, 82]), suggesting a common vicariant history. Such a vicariant event was most likely climatically-mediated as no apparent geographical barrier exists within western Iberian Peninsula.

The western Algarve region in southern Portugal has been indicated as the evolutionary centre for several species and also as a key refugial area [83, 84]. The region harboured relicts of temperate forests during the Last Glacial Maximum [85], probably providing suitable conditions for species survival through glacial periods. The high genetic distance and disjunct distribution found between phylogroups L2 and L3 is most likely the consequence of fragmentation of a once more continuous range in response to a cooling climate. Subsequent climate amelioration during an interglacial period probably resulted in this distribution gap being colonized by L5, expanding from another nearby refugium. A similar distribution pattern of genetic variation is found within Discoglossus galganoi across Portugal, with two closely related phylogroups distributed in the south and north and a less related phylogroup bisecting their distribution [86]. The distribution of Discoglossus galganoi phylogroups and the evolutionary relationship between them is remarkably similar to that found within Lacerta lepida.

Phylogroups L1, L4 and L5

The refugia for Lacerta lepida phylogroups L1, L4 and L5 were probably located in the south-eastern side of the Guadalquivir basin. Support for this comes from the distribution of ancestral haplotypes 127, 140 and 8 within phylogroup L4 that are found only in this region (localities 10, 12, 32, 36, 37, 38, 39, 40, 41). The two very divergent haplotypes within L4 (60 and 61) are also restricted to this region. The widespread distribution of the most frequently sampled haplotype 1 suggests that the spatial and demographic expansion detected within L4 was of a "leading edge" type [87, 88], with few individuals rapidly colonizing adjacent regions, leading to a decrease in genetic diversity on the newly colonized areas. The different phylogroups are most likely the result of three different expansions from the southern refugia dominated by different ancestral haplotypes, followed by further divergence in allopatry.

Refugial areas, range expansions, and secondary contact

The lower mutation rate and slower fixation rate of nuclear genes, when compared to the mitochondrial genome, mean that nuclear genealogies may be more indicative of older demographic events [89]. Therefore, the analysis of the distribution of ancestral nuclear alleles can be useful in the identification of refugial areas. Haplotype B15 is the root of the nuclear gene network, representing the most ancestral allele in the dataset. The distribution of B15 is primarily in the southern region of Iberia, where it occurs with higher frequency than in northern regions (70% of samples with haplotype B15 are from southern latitudes). Furthermore the southern region of Iberia (to the south of river Tagus) shows higher nuclear haplotype richness, with 90% of the identified haplotypes occurring in this region from which 67% do not occur further north. The high nuclear haplotype richness found in the southern regions of Iberia, together with the high frequency of B15 in this region suggest that southern populations are older and the source of subsequent northern expansions.

The geographical and genealogical relationships of mitochondrial haplotypes within phylogroups provide further indications as to the probable refugial areas from where range expansions subsequently occurred. The same pattern of older (interior) haplotypes being typically found within southern sampling sites is also evident for the mitochondrial dataset. For example within phylogroup L2, the ancestral haplotype 96 and the related interior haplotype 93 are most frequent in southern sites 25, 63, 64, 66, 69 and 70, coincident with the Algarve and southern Alentejo region in Portugal. Within phylogroup L3 the ancestral haplotype 51 and the related interior haplotype 46 occur either in the southern region of the phylogroup distribution, to the south of Douro river (sites 82 and 78), or in sites located at the centre of the phylogroup distribution (sites 87, 89 and 92). Within phylogroup L4, the ancestral haplotype 127 and the closely related haplotypes 140 and 8 occur again only in the south-eastern limit of this phylogroup (sites 10, 12, 32, 36, 37, 38, 39, 40 and 41), near the source of the Guadalquivir river in Jaen, Spain. Finally the ancestral haplotype of phylogroup L5 (haplotype 13) and the related interior haplotypes 14 and 17 are most frequent within southern sites of this phylogroup (sites 22, 95, 96 and 99), along the river Tagus basin. These patterns are statistically significant for phylogroups L3 and L5 (Figure 8), which also reveal an excess of derived haplotypes in the north of their geographic ranges. These repeated pattern of ancestral haplotypes occurring more frequently towards the southern range of phylogroups, and in some cases derived haplotypes being more frequent within more northern sites, suggests that southern regions within phylogroup ranges have probably functioned as refugia during past adverse climatic conditions, from which subsequent northern expansions occurred.

The evolutionary history of Lacerta lepida is consistent with a history of population fragmentation and allopatric divergence in southern refugia, followed by recent demographic and spatial range expansions. Demographic expansions are supported by the neutrality and mismatch distribution tests that reveal signatures of demographic and spatial expansions in almost all phylogroups and by the BSP analysis, which show a trend of population expansion in all phylogroups between 0.1 to 0.15 million years ago (Figure 7). Spatial expansions are supported by the latitudinal variation observed for the distributions of ancestral and derived haplotypes. This history of diversification in allopatry followed by range expansions has lead to the establishment of at least four secondary contact zones where very divergent mitochondrial haplotypes belonging to different phylogroups were found in sympatry (Figure 4). Further analyses of these zones may provide insights into the mechanisms involved in speciation and divergence within this lizard species.


Mitochondrial and nuclear gene genealogies in Lacerta lepida provide evidence for a history of isolation and divergence in allopatry resulting in the diversification of six genetically and geographically distinct lineages. Although diversification within the group is largely concordant with the onset of the major glaciations at the beginning of the Pleistocene (approximately 2 Mya), an earlier event associated with the Miocene was also identified. This event (approximately 9 Mya), which marks the divergence of lineages N and L, seems to be associated with geological events related to the evolution of the Mediterranean basin. Signatures of recent demographic and spatial expansions were apparent in all phylogroups, with several phylogroups having established zones of secondary contact. Analyses of the distribution of ancestral and derived alleles within each phylogroup, and inferences related to the biogeography of Lacerta lepida, allowed the identification of probable refugia within the Iberian Peninsula, suggesting southern refugial areas within phylogroups. Our work highlights the importance of these areas for the long-term conservation and management of diversity with Lacerta lepida across its geographic range.


  1. 1.

    Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society B: Biological Sciences. 2004, 359: 183-195. 10.1098/rstb.2003.1388.

    CAS  Article  Google Scholar 

  2. 2.

    Cooper SJB, Hewitt GM: Nuclear DNA sequence divergence between parapatric subspecies of the grasshopper Chorthippus parallelus. Insect Molecular Biology. 1993, 2: 185-194. 10.1111/j.1365-2583.1993.tb00138.x.

    CAS  Article  PubMed  Google Scholar 

  3. 3.

    Hewitt GM: Some genetic consequences of ice ages, and their role, in divergence and speciation. Biological Journal of the Linnean Society. 1996, 58: 247-276.

    Article  Google Scholar 

  4. 4.

    Hewitt GM: After the Ice: Parallelus meets Erythropus in the Pyrenees. Hybrid zones and the evolutionary process. Edited by: Harrison RG. 1993, New York: Oxford University Press

    Google Scholar 

  5. 5.

    Gomez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of Southern European Refugia. Edited by: Weiss S, Ferrand N. 2007, Dordrecht: Springer

    Google Scholar 

  6. 6.

    Alexandrino J, Froufe E, Arntzen JW, Ferrand N: Genetic subdivision, glacial refugia and postglacial recolonization in the golden-striped salamander, Chioglossa lusitanica (Amphibia: Urodela). Molecular Ecology. 2000, 9: 771-781. 10.1046/j.1365-294x.2000.00931.x.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Alexandrino J, Arntzen JW, Ferrand N: Nested clade analysis and the genetic evidence for population expansion in the phylogeography of the golden-striped salamander, Chioglossa lusitanica (Amphibia: Urodela). Heredity. 2002, 88: 66-74. 10.1038/sj.hdy.6800010.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Sequeira F, Alexandrino J, Rocha S, Arntzen JW, Ferrand N: Genetic exchange across a hybrid zone within the Iberian endemic golden-striped salamander, Chioglossa lusitanica. Molecular Ecology. 2005, 14: 245-254.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Paulo OS, Dias C, Bruford MW, Jordan WC, Nichols RA: The persistence of Pliocene populations through the Pleistocene climatic cycles: evidence from the phylogeography of an Iberian lizard. Proceedings of the Royal Society B: Biological Sciences. 2001, 268: 1625-1630. 10.1098/rspb.2001.1706.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Godinho R, Crespo EG, Ferrand N: The limits of mtDNA phylogeography: complex patterns of population history in a highly structured Iberian lizard are only revealed by the use of nuclear markers. Molecular Ecology. 2008, 17: 4670-4683. 10.1111/j.1365-294X.2008.03929.x.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Martínez-Solano I, Teixeira J, Buckley D, Garcia-Paris M: Mitochondrial DNA phylogeography of Lissotriton boscai (Caudata, Salamandridae): evidence for old, multiple refugia in an Iberian endemic. Molecular Ecology. 2006, 15: 3375-3388. 10.1111/j.1365-294X.2006.03013.x.

    Article  PubMed  Google Scholar 

  12. 12.

    Harris DJ, Sá-Sousa P: Molecular phylogenetics of Iberian Wall lizards (Podarcis): Is Podarcis hispanica a species complex?. Molecular Phylogenetics and Evolution. 2002, 23: 75-81. 10.1006/mpev.2001.1079.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Pinho C, Harris DJ, Ferrand N: Non-equilibrium estimates of gene flow inferred from nuclear genealogies suggest that Iberian and North African wall lizards (Podarcis spp.) are an assemblage of incipient species. BMC Evolutionary Biology. 2008, 8: 63-10.1186/1471-2148-8-63.

    Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Martínez-Solano I, Gonçalves HA, Arntzen JW, García-París M: Phylogenetic relationships and biogeography of midwife toads (Discoglossidae: Alytes). Journal of Biogeography. 2004, 31: 603-618. 10.1046/j.1365-2699.2003.01033.x.

    Article  Google Scholar 

  15. 15.

    Branco M, Ferrand N, Monnerot M: Phylogeography of the European rabbit (Oryctolagus cuniculus) in the Iberian Peninsula inferred from RFLP analysis of the cytochrome b gene. Heredity. 2000, 85: 307-317. 10.1046/j.1365-2540.2000.00756.x.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Branco M, Monnerot M, Ferrand N, Templeton AR: Postglacial dispersal of the european rabbit (Oryctolagus cuniculus) on the Iberian Peninsula reconstructed from nested clade and mismatch distribution analyses of mitochondrial DNA variation. Evolution. 2002, 56: 792-803.

    Article  PubMed  Google Scholar 

  17. 17.

    Paulo OS, Pinheiro J, Miraldo A, Bruford MW, Jordan WC, Nichols RA: The role of vicariance vs. dispersal in shaping genetic patterns in ocellated lizard species in the western Mediterranean. Molecular Ecology. 2008, 17: 1535-1551. 10.1111/j.1365-294X.2008.03706.x.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Ujvari B, Dowton M, Madsen T: Population genetic structure, gene flow and sex-biased dispersal in frillneck lizards (Chlamydosaurus kingii). Molecular Ecology. 2008, 17: 3557-3564.

    CAS  PubMed  Google Scholar 

  19. 19.

    Lindell J, Méndez-de la Cruz FR, Murphy RW: Deep biogeographical history and cytonuclear discordance in the black-tailed brush lizard (Urosaurus nigricaudus) of Baja California. Biological Journal of the Linnean Society. 2008, 94: 89-104. 10.1111/j.1095-8312.2008.00976.x.

    Article  Google Scholar 

  20. 20.

    Dowling DK, Friberg U, Lindell J: Evolutionary implications of non-neutral mitochondrial genetic variation. Trends in Ecology & Evolution. 2008, 23: 546-554. 10.1016/j.tree.2008.05.011.

    Article  Google Scholar 

  21. 21.

    Zink RM, Barrowclough GF: Mitochondrial DNA under siege in avian phylogeography. Molecular Ecology. 2008, 17: 2107-2121. 10.1111/j.1365-294X.2008.03737.x.

    CAS  Article  PubMed  Google Scholar 

  22. 22.

    Leaché AD, McGuire JA: Phylogenetic relationships of horned lizards (Phrynosoma) based on nuclear and mitochondrial data: Evidence for a misleading mitochondrial gene tree. Molecular Phylogenetics and Evolution. 2006, 39: 628-644. 10.1016/j.ympev.2005.12.016.

    Article  PubMed  Google Scholar 

  23. 23.

    McGuire JA, Linkem CW, Koo MS, Hutchison DW, Kristopher A, David L, Orange I, Lemos-Espinal J, Riddle BR, Jaeger JR: Mitochondrial introgression and incomplete lineage sorting through space and time: phylogenetics of Crotaphytid lizards. Evolution. 2007, 61: 2879-2897. 10.1111/j.1558-5646.2007.00239.x.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Thorpe RS, Surget-Groba Y, Johansson H: The relative importance of ecology and geographic isolation for speciation in anoles. Philosophical Transactions of the Royal Society B: Biological Sciences. 2008, 363: 3071-3081. 10.1098/rstb.2008.0077.

    Article  Google Scholar 

  25. 25.

    Aljanabi SM, Martinez I: Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Research. 1997, 25: 4692-4693. 10.1093/nar/25.22.4692.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Sunnucks P, Hales DF: Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae). Molecular Biology and Evolution. 1996, 13: 510-524.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Kocher TD, Thomas WK, Meyer A, Edwards SV, Paabo S, Villablanca FX, Wilson AC: Dynamics of mitochondrial-DNA evolution in animals - amplification and sequencing with conserved primers. Proceedings of the National Academy of Sciences of the United States of America. 1989, 86: 6196-6200. 10.1073/pnas.86.16.6196.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Moritz C, Schneider CJ, Wake DB: Evolutionary relationships within the Ensatina-Eschscholtzii complex confirm the ring species interpretation. Systematic Biology. 1992, 41: 273-291.

    Article  Google Scholar 

  29. 29.

    Prychtko TM, Moore WM: The utility of DNA sequences of an intron from the beta-fibrinogen gene in phylogenetic analysis of woodpeckers (Aves: Picidae). Molecular Phylogenetics and Evolution. 1997, 8: 193-204. 10.1006/mpev.1997.0420.

    Article  Google Scholar 

  30. 30.

    Dolman G, Phillips B: Single copy nuclear DNA markers characterized for comparative phylogeography in Australian wet tropics rainforest skinks. Molecular Ecology Notes. 2004, 4: 185-187. 10.1111/j.1471-8286.2004.00609.x.

    CAS  Article  Google Scholar 

  31. 31.

    Godinho R, Mendonca B, Crespo EG, Ferrand N: Genealogy of the nuclear beta-fibrinogen locus in a highly structured lizard species: comparison with mtDNA and evidence for intragenic recombination in the hybrid zone. Heredity. 2006, 96: 454-463. 10.1038/sj.hdy.6800823.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Sequeira F, Ferrand N, Harris DJ: Assessing the phylogenetic signal of the nuclear β-Fibrinogen intron 7 in salamandrids (Amphibia: Salamandridae). Amphibia-Reptilia. 2006, 27: 409-418. 10.1163/156853806778190114.

    Article  Google Scholar 

  33. 33.

    Hall TA: BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows95/98/NT. Nucleic Acids Symposium Series. 1999, 41: 95-98.

    CAS  Google Scholar 

  34. 34.

    Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data. The American Journal of Human Genetics. 2001, 68: 978-989. 10.1086/319501.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    Stephens M, Scheet P: Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. The American Journal of Human Genetics. 2005, 76: 449-462. 10.1086/428594.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Garrick R, Sunnucks P, Dyer R: Nuclear gene phylogeography using PHASE: dealing with unresolved genotypes, lost alleles, and systematic bias in parameter estimation. BMC Evolutionary Biology. 10: 118-

  37. 37.

    Martin DP, Williamson C, Posada D: RDP2: recombination detection and analysis from sequence alignments. Bioinformatics. 2005, 21: 260-262. 10.1093/bioinformatics/bth490.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Martin DP, Rybicki E: RDP: detection of recombination amongst aligned sequences. Bioinformatics. 2000, 16: 562-563. 10.1093/bioinformatics/16.6.562.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Padidam M, Sawyer S, Fauquet CM: Possible emergence of new geminiviruses by frequent recombination. Virology. 1999, 265: 218-225. 10.1006/viro.1999.0056.

    CAS  Article  PubMed  Google Scholar 

  40. 40.

    Posada D, Crandall KA: Evaluation of methods for detecting recombination from DNA sequences: computer simulations. Proc Natl Acad Sci USA. 2001, 98: 13757-13762. 10.1073/pnas.241370698.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Smith JM: Analyzing the mosaic structure of genes. Journal of Molecular Evolution. 1992, 34: 126-129.

    CAS  PubMed  Google Scholar 

  42. 42.

    Gibbs MJ, Armstrong JS, Gibbs AJ: Sister-Scanning: a Monte Carlo procedure for assessing signals in recombinant sequences. Bioinformatics. 2000, 16: 573-582. 10.1093/bioinformatics/16.7.573.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Huelsenbeck JP, Ronquist F: MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics. 2001, 17: 754-755. 10.1093/bioinformatics/17.8.754.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.

    CAS  Article  PubMed  Google Scholar 

  46. 46.

    Akaike H: Information theory and an extension of the maximum likelihood principle. Proceedings of the 2nd International Symposium on Information Theory; Akadémia Kiado, Budapest. Edited by: Petrov BN, Csaki F. 1973, 267-281.

    Google Scholar 

  47. 47.

    Bandelt HJ, Forster P, Rohl A: Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution. 1999, 16: 37-48.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    Templeton AR, Crandall KA, Sing CF: A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992, 132: 619-633.

    CAS  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Molecular Ecology. 2000, 9: 1657-1659. 10.1046/j.1365-294x.2000.01020.x.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Crandall KA, Templeton AR: Empirical tests of some predictions from coalescent theory with applications to intraspecific phylogeny reconstruction. Genetics. 1993, 134: 959-969.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Perera A, Harris DJ: Genetic variability in the Ocellated lizard Timon tangitanus in Morocco. African Zoology. 45: 321-329.

  52. 52.

    Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology. 2007, 7: 214-10.1186/1471-2148-7-214.

    Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Rambaut A, Drummond A: BEAUTi v1.4.2. Bayesian Evolutionary Analysis Utility. 2007

    Google Scholar 

  54. 54.

    Emerson BC: Alarm bells for the molecular clock? No support for Ho et al.'s model of time-dependent molecular rate estimates. Systematic Biology. 2007, 56: 337-345. 10.1080/10635150701258795.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Hasegawa M, Kishino H, Yano TA: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution. 1985, 22: 160-174. 10.1007/BF02101694.

    CAS  Article  PubMed  Google Scholar 

  56. 56.

    Thorpe RS, McGregor DP, Cumming AM, Jordan WC: DNA evolution and colonization sequence of island lizards in relation to geological history: mtDNA, RFLP, Cytochrome B, Cytochrome Oxidase, 12s RRNA sequence, and Nuclear RAPD analysis. Evolution. 1994, 48: 230-240. 10.2307/2410090.

    CAS  Article  Google Scholar 

  57. 57.

    González P, Pinto F, Nogales M, Jiménez-asensio J, Hernández M, Cabrera VM: Phylogenetic relationships of the Canary Islands endemic lizard genus Gallotia (Sauria: Lacertidae), inferred from mitochondrial DNA sequences. Molecular Phylogenetics and Evolution. 1996, 6: 63-71. 10.1006/mpev.1996.0058.

    Article  PubMed  Google Scholar 

  58. 58.

    Martin AP, Palumbi SR: Body size, metabolic rate, generation time, and the molecular clock. Proceedings of the National Academy of Sciences of the United States of America. 1993, 90: 4087-4091. 10.1073/pnas.90.9.4087.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Rambaut A, Drummond A: Tracer v1.4. 2007, []

    Google Scholar 

  60. 60.

    Saillard J, Forster P, Lynnerup N, Bandelt HJ, Nørby S: mtDNA variation among Greenland Eskimos: the edge of the Beringian expansion. The American Journal of Human Genetics. 2000, 67: 718-726. 10.1086/303038.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.

    CAS  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Ramos-Onsins SE, Rozas J: Statistical properties of new neutrality tests against population growth. Molecular Biology and Evolution. 2002, 19: 2092-2100.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147: 915-925.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Slatkin M, Hudson RR: Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991, 129: 555-562.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Molecular Biology and Evolution. 1992, 9: 552-569.

    CAS  PubMed  Google Scholar 

  66. 66.

    Harpending HC, Sherry ST, Rogers AR, Stoneking M: The genetic structure of ancient human populations. Current Anthropology. 1993, 34: 483-496. 10.1086/204195.

    Article  Google Scholar 

  67. 67.

    Harpending HC: Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Human Biology. 1994, 66: 591-600.

    CAS  PubMed  Google Scholar 

  68. 68.

    Excoffier L, Laval G, Schneider S: Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online. 2005, 47-50.

    Google Scholar 

  69. 69.

    Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003, 19: 2496-2497. 10.1093/bioinformatics/btg359.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Drummond AJ, Rambaut A, Shapiro B, Pybus OG: Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution. 2005, 22: 1185-1192. 10.1093/molbev/msi103.

    CAS  Article  PubMed  Google Scholar 

  71. 71.

    Posada D, Crandall KA, Templeton AR: GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Molecular Ecology. 2000, 9: 487-488. 10.1046/j.1365-294x.2000.00887.x.

    CAS  Article  PubMed  Google Scholar 

  72. 72.

    Wakeley J: Coalescent theory: An introduction. 2008, Greenwood Village, Colorado: Roberts & Company Publishers

    Google Scholar 

  73. 73.

    Moore WS: Inferring phylogenies from mtDNA variation: mitochondrial-gene trees versus nuclear-gene trees. Evolution. 1995, 49: 718-726. 10.2307/2410325.

    Article  Google Scholar 

  74. 74.

    Mateo JA, Castroviejo J: Variation morphologique et revision taxonomique de l'espece Lacerta lepida Daudin, 1802 (Sauria, Lacertidae). Bulletin du Museé de Histoire Naturele de Paris. 1990, 12: 691-706.

    Google Scholar 

  75. 75.

    Mateo JA, López-Jurado LF: Variaciones en el color de los lagartos ocelados; aproximacion a la distribuicion de Lacerta lepida nevadensis Buchholz 1963. Revista Espanola de Herpetologia. 1994, 8: 29-35.

    Google Scholar 

  76. 76.

    Mateo JA, López-Jurado LF, Guillaume CP: Variabilité électrophorétique et morphologique des lézards ocellés (Lacertidae): un complexe d'espèces de part et d'autre du détroit de Gibraltar. Comptes Rendus de L'Academie des Sciences Serie iii-Sciences de la Vie-Life Sciences. 1996, 319: 737-746.

    Google Scholar 

  77. 77.

    Rosenbaum G, Lister GS, Duboz C: Relative motions of Africa, Iberia and Europe during Alpine orogeny. Tectonophysics. 2002, 359: 117-129. 10.1016/S0040-1951(02)00442-0.

    Article  Google Scholar 

  78. 78.

    Rosenbaum G, Lister GS, Duboz C: Reconstruction of the tectonic evolution of the western Mediterranean since the Oligocene. Journal of the Virtual Explorer. 2002, 8: 107-130.

    Google Scholar 

  79. 79.

    Garcia-Paris M, Alcobendas M, Alberch P: Influence of the Guadalquivir river basin on mitochondrial DNA evolution of Salamandra salamandra (Caudata: Salamandridae) from southern Spain. Copeia. 1998, 1998: 173-176. 10.2307/1447714.

    Article  Google Scholar 

  80. 80.

    Arntzen JW, Garcia-Paris M: Morphological and allozyme studies of midwife toads (genus Alytes), including the description of two new taxa from Spain. Contributions to zoology. 1995, 65: 5-34.

    Google Scholar 

  81. 81.

    Biju-Duval C, Ennafaa H, Dennebouy N, Monnerot M, Mignotte F, Soriguer R, El Gaied A, El Hili A, Monoulou JC: Mitochondrial DNA evolution in Lagomorphs: origin of systematic heteroplasmy and organization of diversity in European rabbits. Journal of Molecular Evolution. 1991, 33: 92-102. 10.1007/BF02100200.

    CAS  Article  Google Scholar 

  82. 82.

    Paulo OS, Jordan WC, Bruford MW, Nichols RA: Using nested clade analysis to assess the history of colonization and the persistence of populations of an Iberian Lizard. Molecular Ecology. 2002, 11: 809-819. 10.1046/j.1365-294X.2002.01484.x.

    CAS  Article  PubMed  Google Scholar 

  83. 83.

    Mesquita N, Hanfling B, Carvalho GR, Coelho MM: Phylogeography of the cyprinid Squalius aradensis and implications for conservation of the endemic freshwater fauna of southern Portugal. Molecular Ecology. 2005, 14: 1939-1954. 10.1111/j.1365-294X.2005.02569.x.

    CAS  Article  PubMed  Google Scholar 

  84. 84.

    Fritz U, Barata M, Busack SD, Fritzsch G, Castilho R: Impact of mountain chains, sea straits and peripheral populations on genetic and taxonomic structure of a freshwater turtle, Mauremys leprosa (Reptilia, Testudines, Geoemydidae). Zoologica Scripta. 2006, 35: 97-108. 10.1111/j.1463-6409.2005.00218.x.

    Article  Google Scholar 

  85. 85.

    Zagwijn WH: Migration of vegetation during the Quaternary in Europe. Courier Forschungsinstitut Senckenberg. 1992, 153: 9-20.

    Google Scholar 

  86. 86.

    Martínez-Solano I: Phylogeography of Iberian Discoglossus (Anura: Discoglossidae). Journal of Zoological Systematics & Evolutionary Research. 2004, 42: 298-305. 10.1111/j.1439-0469.2004.00257.x.

    Article  Google Scholar 

  87. 87.

    Hewitt GM: The subdivision of species by hybrid zones. Speciation and its consequences. Edited by: Otte D, Endler J. 1989, Sunderland, Massachusetts: Sinauer Associates, 85-110.

    Google Scholar 

  88. 88.

    Hewitt GM: Postglacial distribution and species substructure: lessons from pollen, insects and hybrid zones. Evolutionary patterns and processes. Edited by: Lees DR, Edwards D. 1993, London: Academic Press, 14: 97-123.

    Google Scholar 

  89. 89.

    Avise JC: Molecular Markers, Natural History and Evolution. 2004, Sunderland, Massachusetts: Sinauer Associates, 2

    Google Scholar 

Download references


This work was supported by a grant from FCT, Portugal (POCI/BIA-BDE/59288/2004). AM was supported by a PhD grant from FCT (SFRH/BD/16996/2004). We thank S. Carranza (University of Barcelona) and J. Pleguezuelos (University of Granada) for providing several lizard samples. We thank Juan Pablo de la Vega, Luis Garcia, Gabriel Marin, Eugenia Zarza-Franco, Rui Osório, Pedro Figueiredo, Nuno Valente and Matthew Osborne for helping with fieldwork. We also thank Mário Miraldo and António Firmino for assembling the lizards' traps. All samples were collected with appropriate licenses in each country (Instituto de Conservação da Natureza in Portugal and the Regional Ministry of Environment of the Government of Andalucia, Asturias, Galicia, Extremadura, Castilla y Leon, Murcia, Castilla la Mancha and Madrid in Spain).

Author information



Corresponding author

Correspondence to Andreia Miraldo.

Additional information

Authors' contributions

AM conceived the study, carried out fieldwork, molecular genetic work, data analysis and drafted the manuscript. OSP participated in the study design and helped to draft the manuscript. BCE and GMH participated in the study design, helped to draft the manuscript and coordinated the overall development of the study. All authors read and approved the final manuscript.

Electronic supplementary material


Additional file 1: Table S1. This file includes one table (Table S1) with information about sampling localities, number of samples per locality and haplotypes detected in each locality. (DOC 223 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

Miraldo, A., Hewitt, G.M., Paulo, O.S. et al. Phylogeography and demographic history of Lacerta lepida in the Iberian Peninsula: multiple refugia, range expansions and secondary contact zones. BMC Evol Biol 11, 170 (2011).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Iberian Peninsula
  • Mismatch Distribution
  • Refugial Area
  • Much Recent Common Ancestor
  • Ancestral Haplotype