Multilocus phylogeography of the common lizard Zootoca vivipara at the Ibero-Pyrenean suture zone reveals lowland barriers and high-elevation introgression

Background The geographic distribution of evolutionary lineages and the patterns of gene flow upon secondary contact provide insight into the process of divergence and speciation. We explore the evolutionary history of the common lizard Zootoca vivipara (= Lacerta vivipara) in the Iberian Peninsula and test the role of the Pyrenees and the Cantabrian Mountains in restricting gene flow and driving lineage isolation and divergence. We also assess patterns of introgression among lineages upon secondary contact, and test for the role of high-elevation trans-mountain colonisations in explaining spatial patterns of genetic diversity. We use mtDNA sequence data and genome-wide AFLP loci to reconstruct phylogenetic relationships among lineages, and measure genetic structure. Results The main genetic split in mtDNA corresponds generally to the French and Spanish sides of the Pyrenees as previously reported, in contrast to genome-wide AFLP data, which show a major division between NW Spain and the rest. Both types of markers support the existence of four distinct and geographically congruent genetic groups, which are consistent with major topographic barriers. Both datasets reveal the presence of three independent contact zones between lineages in the Pyrenean region, one in the Basque lowlands, one in the low-elevation mountains of the western Pyrenees, and one in the French side of the central Pyrenees. The latter shows genetic evidence of a recent, high-altitude trans-Pyrenean incursion from Spain into France. Conclusions The distribution and age of major lineages is consistent with a Pleistocene origin and a role for both the Pyrenees and the Cantabrian Mountains in driving isolation and differentiation of Z. vivipara lineages at large geographic scales. However, mountain ranges are not always effective barriers to dispersal, and have not prevented a recent high-elevation trans-Pyrenean incursion that has led to asymmetrical introgression among divergent lineages. Cytonuclear discordance in patterns of genetic structure and introgression at contact zones suggests selection may be involved at various scales. Suture zones are important areas for the study of lineage formation and speciation, and our results show that biogeographic barriers can yield markedly different phylogeographic patterns in different vertebrate and invertebrate taxa.


Background
The geographical distribution of intraspecific lineages and the pattern of gene flow among them provides valuable insight into the process of lineage divergence and speciation [1]. Of particular interest are phylogeographic studies at large spatial scales that encompass major geographic barriers to gene flow as well as areas where divergent lineages come into secondary contact. Patterns of introgression at these contact zones can shed light on the degree of reproductive isolation among lineages and the historical factors driving divergence [2]. In Europe, glacial dynamics over the last several million years had a major impact on the phylogeography of animal and plant lineages, which underwent repeated cycles of rapid expansions from southern refugia as species recolonized northern latitudes following glacial maxima [3][4][5]. These range expansions gave rise to a series of suture zones across the continent as lineages that diverged during glacial periods came into secondary contact during interglacials, providing us with unique natural experimental areas in which to study the process of lineage divergence and the evolution of reproductive isolation [5].
One of these regions of secondary contact is the Ibero-Pyrenean suture zone, an important biogeographic region where both inter and intra-specific lineages are known to come into contact in a number of taxa [3][4][5][6]. In this area of southwestern Europe, the Pyrenees stand as a major barrier separating central Europe from the Iberian peninsula, one of the major glacial refugia in southern Europe, and also a geologically and ecologically complex area that has suffered itself severe climatic changes, giving rise to a complex phylogeographic pattern that consists of "refugia within refugia" [7]. Here we examine the phylogeography and patterns of gene flow among lineages of the common lizard Zootoca vivipara (= Lacerta vivipara) in the Ibero-Pyrenean region. The species has a broad Eurasian distribution composed largely of viviparous lineages, yet individuals in this region belong to an oviparous lineage isolated from the nearest viviparous populations in the French Massif Central by a gap of unsuitable habitat [8]. Recent studies of Z. vivipara in this region have documented a sharp contact zone between two divergent maternal lineages corresponding generally to the French and Spanish sides of the Pyrenees, and showing consistent results with mtDNA 16S and cytochrome b gene sequences [8,9] and a maternally inherited marker on the female W chromosome [10]. However, because of limited regional sampling in previous phylogeographic studies, and the lack of biparentally inherited markers examined in this system to date, our understanding of region-wide genetic structure, patterns of introgression among lineages and the relative importance of the Pyrenees and other mountain ranges as barriers to gene flow remain general and tentative.
Inferring the demographic and evolutionary history of lineages at various spatial and temporal scales often requires the use of molecular markers from genomic regions with different modes of inheritance and rates of evolution [11]. Mitochondrial DNA sequence markers have proven useful for studying intraspecific genetic variation because of their relatively high mutation rate, small effective population size, short coalescence times, and matrilineal mode of inheritance [1,12,13]. However, mt-DNA genomes are susceptible to the effect of historical factors such as population expansions [3,14], introgressive hybridization [15,16], incomplete lineage sorting [17] and even selection [18,19], all of which can effectively mask real patterns of gene flow among populations and lineages, potentially confounding inferences of evolutionary and demographic events [20]. Complementing mtDNA datasets with biparentally inherited nuclear DNA (nDNA) marker data is often essential to obtain a more complete picture of lineage history [11,21].
Here we undertake a multilocus phylogeographic survey across the region using both mtDNA sequences and genome-wide AFLP loci and an improved geographic sampling in order to document the presence and distribution of major genetic lineages and test the role of the Pyrenees and the Cantabrian Mountains in driving lineage divergence in Z. vivipara. We also use this largescale phylogeographic context to assess patterns of introgression and test hypotheses related to the history of contact zones.

Mitochondrial DNA structure and diversity
Sequencing of 200 individuals of Z. vivipara from 48 localities (Table 1) revealed 14 cyt-b and 29 ND2 haplotypes, for a total of 33 haplotypes for the concatenated dataset (1,370 bp).
Phylogenetic analysis revealed two major clades corresponding generally to southern France and northern Spain, respectively (Figures 1 & 2A). The Spanish clade is further divided into three major subclades, including a NW Spain clade (Western Cantabrian Mountains), a north-central Spain clade (Eastern Cantabrian Mountains and Basque Country) and a NE Spain clade (Pyrenees) (Figures 1 & 2A). Lizards carrying haplotypes from the latter clade were also found on the French slope of the Pyrenees, where they were found sympatrically with the divergent French clade haplotypes (Figures 2A and  3A). Percent divergence in corrected distances between the two major clades was 2%, and average divergence among northern Spain clades was 0.6% ( Table 2). Dating of the main clades using Bayesian inference in BEAST, indicate that the main France-Spain split (node 1 on  Figure 1A) took place about 0.935 My ago (95% HPD: 0.396-1.591), whereas the divergence among clades in northern Spain ranges from 0,408 My (95% HPD: 0.166-0.713) for node 2, to 0.169 My (95% HPD: 0.056-0.307) for node 3 ( Figure 1A). Patterns of genetic diversity revealed that the NE Spain clade (blue in Figures 1 & 2) showed lower haplotypic and nucleotide diversity values than the rest ( Table 3).
The prominent "starlike" phylogenetic pattern in this clade, with a single high-frequency haplotype ("AA") and various lower-frequency, closely-related haplotypes, is suggestive of a recent population expansion [22]. This was corroborated by results from Fu's test of population expansion, which revealed significant values of F s , consistent with a recent burst of population growth ( Table 3). The pattern of haplotype frequency and distribution of the  NE Spain (blue) clade reflects higher haplotype diversity on the Spanish side of the Pyrenees than on the French side, where the relative frequency of haplotype "AA" is more pronounced. This pattern, together with the evidence for a recent population expansion in this lineage, suggests that the presence of Spanish haplotypes on the French side is the result of a trans-Pyrenean colonization of the French side. We estimated time since the population expansion from the distribution of pairwise differences among blue-clade haplotypes, which yielded a value of τ = 3.00 (95% CI: 0.00-3.83). Applying a mutation rate of 0.01 s/s/Myr per lineage, this value corresponds to a time since the expansion of 54,744 years, with confidence intervals between the present and 69,890 years ago.

AFLP structure and diversity
A total of 34 genome-wide amplified fragments were unambiguously scored for 223 individual lizards using 3 pairs of selective primers. Out of the total 34 loci, 53% to 82% were polymorphic depending on the population ( Table 4). Analysis of AFLP variation revealed marked geographic structure that shows partial congruence but  Table 1). Colours correspond to the major clades from the phylogenetic tree on Figure 1. In the French Pyrenees population six haplotypes with a frequency of one are unlabelled: S1L, YL and ZN (green tones) and AD, AE and DA (blue tones). (B) Genetic variation in genome-wide AFLP markers. Each pie represents the posterior probabilities of assignment to four main genetic clusters (K = 4) obtained from STRUCTURE 3.1, averaged across individuals at each locality. The size of each pie graph is proportional to sample size. Black dots represent sampling localities, and white frames represent areas including several localities, shown in detail on Figure 3. The dashed contour line represents the approximate distribution range of the species in this region.   Pyrenees, and K = 5 revealed small additional groups to those in K = 4 ( Figure 5). The Evanno method for obtaining the "optimal" K value gave the highest probability to K = 2, with higher values being less likely (ΔK values: 222.82 for K = 2, 72.14 for K = 3, 50.72 for K = 4, and 22.72 for K = 5, Additional file 1: Figure S1A). In contrast, a plot of mean values of the estimated Ln probability of K suggests that K = 4 could correspond to the most realistic level of structure (Additional file 1: Figure S1B), as it marks the start of the flattening of the curve. Determining the optimal K is not always straightforward and consideration of both probability scores and biological information is recommended when assessing the "true" number of populations [23]. Interestingly, K = 4 showed a remarkable geographic correspondence with the four main clades in the mtDNA phylogenetic tree (Figures 2 and 5), strongly suggesting that this K value for AFLP data is biologically meaningful. In any case, the geographic congruence between K = 4 AFLP clusters and the mtDNA phylogenetic tree provides a unique opportunity to compare patterns of divergence between the two types of markers and examine the behaviour of each one at Z. vivipara contact zones in the region. With respect to genetic differentiation between AFLP genetic clusters (using K = 4), Φ PT between NW Spain (yellow lineage in Figure 2) and all remaining populations was 0.30 (P < 0.001), with Φ PT values between the yellow cluster and other eastern clusters ranging from 0.32 to 0.41 ( Table 5).
None of the 34 loci showed significant departures from neutrality in the BAYSCAN 2.0 analysis across all comparisons, although one of them (locus "15c194") was identified as an outlier in the K = 2 comparison, and another one (locus "19c84") in the K = 3 and K = 4 comparisons (Additional file 1: Figure S2). Excluding these two loci from the dataset yielded results for the Structure analysis and the PCoA analysis in GENALEX that were indistinguishable from those using all loci. This is due to the fact that most markers in these comparisons had

Introgression at contact zones
Using the geographic distribution of the main four mt DNA clades and the comparable AFLP-defined genetic units corresponding to K = 4, we detected localities with mixed haplotypes (for mtDNA) or assignment probability  values (for ALFPs) in two main areas within the region, one corresponding to the contact zone between blue and green lineages across the Pyrenees, and one between blue and red lineages in Navarre, just east of the Basque Country ( Figure 2). The contact zone in southern France between divergent blue and green mtDNA lineages, apparent in Figures 2A and 3A and previously described in detail by Heulin et al. [24] and references therein, is shifted to the south according to genome-wide AFLP data. Populations containing genotypes assigned to blue and green genetic clusters are restricted to a limited area within a few kilometres of the Spanish border (upper Ossau River), whereas the genotypes of individuals in the remaining area to the north group with the green southern France cluster ( Figure 3B). At the AFLP contact zone, the degree of nuclear introgression varies per population ( Figure 6A), with most individuals in some localities showing high to complete assignment probability to either the blue or green clusters (AYG and BRO, respectively), whereas in others (SOU, GAB and SOP) individuals show a range of assignment probabilities suggesting a hybrid origin. The ES4 population and other sites located further away from the border show complete assignment probability to the southern France (green) cluster.
At a second mtDNA contact zone detected between the red and blue clades at the Ibañeta site in the Spanish side of the western Pyrenees (Figure 2A), the AFLP data revealed that all individuals had a higher probability of assignment to the blue cluster ( Figure 2B). A more detailed examination of the 10 individuals genotyped from this site shows that six of them had >90% probability of assignment to the blue cluster and less than 10% to the red cluster, and the remaining four individuals had 86-48% probability of assignment to the blue cluster and 10-33% to the red cluster ( Figure 6B Figure 2B. sampled carry southern France mtDNA haplotypes (green clade), yet their AFLP profiles cluster with the NE Spain group (blue cluster) (Figure 2), a situation similar to that of many southern France individuals.

Phylogeography and mtDNA lineage history
Our mtDNA dataset revealed four major clades and a remarkable degree of geographic structuring among Zootoca vivipara lineages. The high variation found in the ND2 gene relative to cyt-b, together with a more thorough geographic sampling allowed us to identify several distinct and geographically segregated clades in northern Spain, which help explain their historical interactions with the divergent clade found in southern France and clarify the phylogeographic history of the region. In line with previous phylogeographic studies on the species, our study shows that populations of Z. vivipara at the Ibero-Pyrenean zone do not show an actual "suture zone" between Iberian and northern European lineages [8,24], in sharp contrast to patterns found in other species, like the grasshopper Chorthippus parallelus [25] or the warbler Phylloscopus collybita [26]. We also find no evidence of Iberian clades expanding northward to colonize northern European latitudes, as documented for toads [27], shrews [28], hedgehogs [29], or bears [30]. Instead, the oviparous clades of Z. vivipara on opposite sides of the Pyrenees form a monophyletic clade with respect to all other European lineages, indicating a relatively long time in isolation in this region. However, our mtDNA data do support a major role for the Pyrenees in isolating populations, and the two main divergent lineages (depicted by node "1" in Figure 1A) appear to have diverged on either side of this major barrier within the last one million years, since around the mid Pleistocene. The Spanish clade further differentiated within the last 500 K years into at least three subclades corresponding to separate mountainous regions, namely the western and eastern sections of the Cantabrian Mountains, and the south-facing slope of the Pyrenees. Further structure is also indicated by the small clade in the western-most end of the range, formed by two haplotypes (XJ and X1J) sampled in the Serra do Xistral, Galicia, although further sampling will be necessary to confirm this pattern.
Overall, the general pattern of allopatric differentiation in isolated highlands is consistent with the species ecology,  Figure 2A, and for the AFLP data each vertical bar corresponds to the percent assignment probability to genetic clusters on Figure 2B.
as in the southern part of its distribution the common lizard is found mainly in humid bogs and associated habitats [31,32], which in Spain are often found in montane areas.
Other lizards associated with humid areas, such as members of the genus Iberolacerta, show a deep phylogenetic break between the NW Spain form (I. monticola) and the Pyrenean taxa (I. bonnali, I. aranica and I. aurelioi), suggesting similar patterns of isolation as in Z. vivipara albeit at longer time scales [33,34]. In contrast, species that do not specialize on humid habitats have shown a contrasting lack of structure across similar regions. For instance, populations of Lacerta schreiberi along northern Spain are genetically uniform and expanded from a refuge in central Portugal since the last glacial maximum [35], a pattern also shown by Podarcis bocagei [36] and the salamander Chioglossa lusitanica [37]. The pattern of distinct Z. vivipara lineages associated with different mountain areas along the northern Iberian peninsula is consistent with a pattern of "refugia-withinrefugia" documented for a number of other organisms [7,38,39], and is consistent with a climatically dynamic history which combined with a complex topography has given rise to allopatrically differentiated populations and lineages. In this context, Z. vivipara lineages are younger than those reported in other taxa, such as Psammod romus hispanicus [39], Alytes obstetricans [40], or Lisso triton boscai [41], and more similar in age to those in Alytes cisternasii [42].
AFLP structure: contrasting patterns with mtDNA The pattern of spatial variation in genome-wide AFLP markers revealed both important similarities and differences compared to that of mtDNA sequence. Both datasets revealed the presence of four genetic units that are geographically concordant, yet relative divergence among groups differs in both datasets. AFLP variation forms two main genetic groups that separate NW Spain from the rest, while the two main mtDNA clades separate Spain from France. This difference suggests that at a large geographic scale different neutral and/or selective factors have shaped diversity and structure in both types of markers, yet identifying the specific cause is beyond the scope of our dataset. Potential causes range from neutral factors such as random sorting of loci or chromosomal rearrangements, to the role of natural selection. Even though the AFLP loci we used appear to be largely neutral, these markers are likely to be spread randomly across the nuclear genome, and are therefore more likely to be associated with regions under selection than mtDNA markers. If this were the case, we would expect to find a stronger association between traits under selection and AFLP structure. Available phenotypic data on Z. vivipara in the region provides some support for this hypothesis, and analysis of morphological traits has revealed marked phenotypic differentiation of Cantabrian populations with respect to those in the Basque Country and Pyrenees [43]. This pattern is driven largely by the shape of the head-shield scales, a trait with unknown fitness value, although other traits related to coloration known to be important in related lacertids [44,45] are yet to be analysed in Z. vivipara. A closer association between traits under selection and nuclear instead of mitochondrial neutral makers has been documented in other studies. For example, in a contact zone between central European newts, Babik et al. [46] document a tighter correlation between nuclear DNA and morphological hybrid indices than between mtDNA and morphology. Also, a recent study in the Ensatina salamander complex of North America shows that divergence in neutral nuclear markers is better correlated with reproductive isolation than mtDNA [47]. Additional ecological and phenotypic data will be necessary to properly test the hypothesis that neutral nuclear markers show a better association to adaptive variation than mitochondrial markers.

Contact zone dynamics and patterns of introgression
Our phylogeographic results reveal the presence of two independent contact zones between lineages: one in southern France across the high Pyrenees (between the blue and green clades in Figure 2) and one in the foothills of the southwestern Pyrenees (between the red and blue clades). An additional contact zone in the Basque lowlands between the red and green clades was not detected here due to limited sampling in that area, but a population with haplotypes from both clades has been reported by Heulin et al. (2011) at Moura de Montrol sites. Patterns of introgression in this area have been studied in detail among lineages of Corthippus grasshoppers [48] and Phylloscopus birds [26,49], both forming suture zones among old lineages. Additional spatial sampling will be necessary to properly compare these suture zones to the secondary contact dynamics among the relatively young lineages of Z. vivipara. Another potential contact zone not detected in the present study corresponds to that found in north-central Spain, between the yellow and red lineages, and further sampling would be necessary to confirm its existence.
Relatively good sampling in the central Pyrenees allows us to infer the demographic and introgression history of lineages there. Different lines of evidence (phylogeny, haplotype frequencies and expansion test) support the hypothesis that the presence of NE Spanish (blue) haplotypes in southern France is due to a recent population expansion from the Spanish side of the Pyrenees, which led to the colonization of southern France across the high Pyrenees. The marked predominance of haplotype "AA" there suggests that this incursion took place recently, probably during or since the last glacial maximum, as indicated by time estimates from the mismatch distribution. The relative geographic congruence among mt DNA and AFLP genetic groups allows us to compare patterns of introgression of maternal and bi-parentally inherited loci. Genome-wide AFLP data revealed limited introgression of nDNA into France, with individuals of mixed origin being restricted to the upper Ossau river drainage and the Spain-France border, suggesting that the contact zone among lineages is located within 10 or 15 km of the border. In contrast, Spanish mtDNA haplotypes are found at least tens of kilometres further than nDNA. Interestingly, it appears that Spanish haplotypes expanded into France (mostly towards the W and NW), they were then trapped at river valleys towards the east (i.e. lower Aspe river and upper Ossau river valleys), giving rise to a sharp contact zone there. This apparent barrier to gene flow in mtDNA but not nDNA at river valleys could be due indirectly to Haldane's rule [50], which predicts a lower fitness for hybrids of the heterogametic sex, females in the case of Z. vivipara. A more direct role of selection through differential survivorship of alternative haplogroups has been previously suggested by Heulin et al. (2011), although recent evidence from the same study area in Gabas-piet (Heulin, unpublished data) indicates that further replication across years will be necessary to obtain conclusive evidence. Several studies in diverse organisms have shown that mtDNA can introgress further than nDNA [20,46,[51][52][53][54], and spread relatively quickly. An alternative explanation for the shifted positions of the mtDNA and nDNA clines is that following the initial expansion, the contact zone has moved back south, leaving behind Spanish mtDNA haplotypes in its wake [55,56]. Distinguishing among these alternative scenarios will require further geographic and molecular sampling, including co-dominant nuclear markers.
Additional fine-scale sampling at contact zones among lineages will also be needed to understand introgression dynamics at other areas, such as the Beret site, where all individuals carry French mtDNA and Spanish nDNA, or the contact zone in the western Pyrenees between the blue and red clades detected at the Ibañeta site, where individuals show blue central-Pyrenean nDNA, yet half of them carry red mtDNA haplotypes. Patterns at these individual sites are similar to those found at some sites in southern France, and a more thorough geographic context provided by finer-scale sampling will be needed to help determine whether they belong to active contact zones with ongoing introgression or instead represent populations left in the genetic wake of a contact zone that shifted away.

Conclusions
Our results show that lineages of Z. vivipara in northern Iberia and the Ibero-Pyrenean suture zone show marked levels of differentiation that is congruent with major topographical features in the region. Our data reveal that mountain ranges do not necessarily represent impassable barriers for this species, as demonstrated by a recent high-elevation colonization event across the Pyrenees that has led to secondary contact and partial introgression among divergent lineages. In contrast, we find indirect evidence of contact zones in lowland areas, in the absence of obvious barriers to dispersal. Our results underscore the importance of multilocus datasets in reconstructing the evolutionary history of independent intraspecific lineages [57], and understanding the dynamics of introgression upon secondary contact. The degree of mixing among Z. vivipara lineages in secondary contact is similar to that found among different lizard species [58], suggesting that Z. vivipara lineages in the region may in fact represent separate species. Further research using both molecular and phenotypic data will help determine the extent to which reproductive isolation is underway despite the presence of gene flow [59]. Suture zones are important areas for the study of lineage formation and speciation, as they provide a geographic context for across-taxa comparisons. However, our results on Z. vivipara show that even major biogeographic barriers can yield markedly different phylogeographic patterns in different vertebrate and invertebrate taxa depending on their demographic histories and ecological traits. Combining intraspecific studies of genetic variation at the lineage level with ecological data will be essential to better understand phylogeographic differences among taxa and advance our understanding of diversification mechanisms.

Field sampling
Specimens of Z. vivipara were captured at 49 localities across the species range in Northern Spain and Southern France ( Table 1). The tail tip of each individual captured was clipped and stored in 95% ethanol at −20°C. Individuals were weighed, measured, photographed, and then released at the site of capture. The capture and handling of lizards complied with national and international ethical guidelines and was conducted under the scientific collecting permits issued by Xunta de Galicia, Junta de Castilla y León, Gobierno del Principado de Asturias, Gobierno de Navarra, Gobierno de Aragón and Generalitat de Catalunya in Spain, and Parc National des Pyrénées in France.

MtDNA sequencing and analysis
Genomic DNA was extracted from tail clippings using a Qiagen™ DNeasy Kit and following the manufacturer's protocol. The following regions of the mitochondrial DNA were amplified: a 427 bp fragment of the cytochrome b gene (including 21 bp from the adjacent Glu-tRNA) using primers MVZ04 and MVZ05 [60]; and 964 bp of the NADH dehydrogenase subunit 2 gene (ND2) using primers MetF6 and AsnR2 [61]. Detailed PCR amplification conditions and sequencing details are available in the Additional file 1.
Sequences were automatically aligned using SEQUEN-CHER 4.1. (GeneCodes, Ann Arbor, Michigan) and variable sites were checked visually for accuracy. Coding gene sequences were unambiguously translated into their amino acid sequence and no double peaks were observed in the chromatographs, suggesting sequences were of mitochondrial origin and not nuclear copies. For most analysis, the two gene fragments (cyt-b and ND2) were concatenated into a single sequence of 1370 bp. Haplotype and nucleotide diversity indices were calculated in DNASP V5 [62].
We constructed a haplotype network using a maximum parsimony algorithm as implemented in the program TCS 1.21 [63], setting minimum branch probability at 95%. We used JMODELTEST [64] to determine the model of sequence evolution for each marker (HKY + G for ND2 and HKY + I for Cyt-b) and obtained an optimal partitioning scheme for the dataset (by gene and by codon position) using the sotware package PartitionFinder [65]. We then constructed a phylogeny of mtDNA haplotypes using Bayesian analysis with MRBAYES 3.2.2 (http:// mrbayes.csit.fsu.edu/) [66], and ran two simultaneous parallel runs of four Markov chains each (3 heated and one cold) for 4 million generations and sampled every 100 generations. The first 25% of the trees were excluded as burn-in, and a consensus topology was obtained from the remaining 60,002 samples (30,001 per run). We confirmed MCMC chain convergence using TRACER v1.5 [67], and output parameters such as an average PSRF value of 1.00, and an average value of 0.0035 for the standard deviation of the split frequencies between simultaneous runs.
To estimate divergence times among lineages we used a coalescence approach that uses Bayesian inference and MCMC simulations to generate posterior probability values for divergence times as implemented in the program BEAST [68]. We partitioned our dataset by gene and codon position, and ran the analysis with unlinked substitution and clock models for each partition, yet generated a single tree from both. We conducted preliminary runs using an uncorrelated relaxed lognormal clock to account for rate heterogeneity among lineages, but obtained a value of zero for the ucld.stdev parameter, indicating that our data are clock-like. We therefore used the strict clock in the final analysis, as well as a coalescent model of diversification, and a UPGMA starting tree. After optimizing the priors with preliminary runs, we conducted two final runs of 20 million generations, sampled every 2000 steps. Chain convergence and burn-in were determined with the program TRACER v1.5 [67], and all ESS values in the final runs were above 600. As a prior for the mutation rate we used a value of 0.0228 subst./site/Myr (0.0114 s/s/Myr per lineage) for Cyt b, a rate estimated from a comprehensive phylogenetic study of several squamate groups [69] that used as time calibration points the known age of the Canary Islands and the opening of the Strait of Gibraltar 5.3 Mya, well known geological events that have been associated with speciation events in several reptile and amphibian species [70,71]. We set a lognormal distribution for the prior with a mean of 0.0114 and a log standard deviation of 0.40, so that values of 0.005 and 0.020 s/s/Myr per lineage corresponded to the 5% and 95% quantiles, respectively, as these values encompass the range of rates used for estimating divergence times in other studies on lizards [33,35,38,72].
We tested for past sudden changes in effective population size using Fu's test of neutrality [73], which detects departures from neutrality in scenarios characterized by an excess of rare alleles and young mutations in sequences not subjected to recombination. Significant large negative values of F s (generated with ARLEQUIN 3.5) indicate an excess of recent mutations and reject population stasis (Fu, 1997). To estimate the time elapsed since such a sudden population expansion, we used the parameter τ obtained from the distribution of pairwise differences among Cyt b haplotypes, as τ = 2ut, where t is the time elapsed between initial and current population sizes and u = 2μk, where μ is the mutation rate (0.01 substitutions/ site/My per lineage) and k is the length of the sequence [74]. We estimated the parameter τ using a generalized non-linear least-square approach [75] and computed confidence intervals for an alpha level of 0.05 by a parametric bootstrap method based on 100 replicates as implemented in ARLEQUIN 3.5 [76].

AFLP profiling and analysis
Amplified fragment length polymorphism (AFLP) profiles were generated using the laboratory protocol described in Milá et al. [77], which was modified slightly from Vos et al. [78]. In brief, whole genomic DNA was digested with restriction enzymes EcoRI and MseI (Tru9) and fragments were ligated to oligonucleotide adapters with T4 DNA ligase. A random sub-sample of fragments was obtained through a pre-selective amplification reaction using primers E T and M C , followed by three selective amplifications using primer pairs E TAG /M CGA , E TAG /M CGT , and E TAG / M CTA , with the E primer fluorescently labelled. Ten pairs of selective amplification primers were tested, but only those producing sufficient loci and repeatable and unambiguously scorable profiles were selected for the analysis. Selectively amplified fragments were run in an ABI 3730XL genetic analyser with a LIZ500 size standard. Peaks were visualized using GENEMAPPER 3.7 and scored manually by a single observer (BM). Only unambiguously scorable loci and individuals were included in the analysis and peaks found in less than 2% of individuals were excluded. Methodological error rate was assessed by running a subset of 10 individuals twice from the pre-selective amplification step. The average per-locus genotyping error rate for the AFLP loci selected, measured as recommended by Bonin et al. [79], was low at 1.2%.
To determine whether AFLP markers were selectively neutral, we conducted an outlier analysis using the program BAYESCAN 2.0, which uses differences in allele frequencies between populations and a reversiblejump MCMC algorithm to estimate the locus-specific posterior probability for a model including selection vs. a neutral model [80]. We used 10,000 iterations, a thinning interval of 100, 40 pilot runs, an additional burnin of 5,000, and a value of 10 for the prior odds for the neutral model [80].
We estimated allelic frequencies using Zhivotovsky's [81] Bayesian method with uniform prior distributions and assuming Hardy-Weinberg genotypic proportions. Genetic diversity (H e ) was based on allele frequencies calculated using the method by Lynch and Milligan [82] as implemented in the program AFLP-SURV v. 1.0 [83]. A matrix of pairwise population F ST values using the F ST analogue Φ PT was calculated with GENALEX 6.0 [84]. Φ PT is calculated as V AP /(V AP + V WP ) where V AP is the variance among populations and V WP is the variance within populations. Probability values of pairwise Φ PT were based on 10,000 permutations.
To assess genetic structure among samples we conducted a principal coordinate analysis (PCoA; Orloci [85]) on a genetic distance matrix generated from the binary presence-absence matrix as implemented in GENALEX 6.0. [84].We further examined patterns of population structure using the Bayesian assignment probability test in the program STRUCTURE 3.1 [86]. This program uses a Bayesian approach to generate posterior probabilities of assignment of individuals to each of a given number of groups (K). As recommended for dominant markers, we applied a model of no admixture with correlated allele frequencies [23], and used the Locprior setting, which incorporates sampling locality as a prior without affecting the value of the optimal K [87]. We run five different chains of 1,000,000 steps for each value of K after a burn-in of 100,000 steps. The optimal value of K was estimated by examining the mean Ln likelihood values for different K values from 1 to 10, and by following the method by Evanno et al. [88] as implemented in STRUCTURE HAR-VESTER [89].