Extreme genetic diversity in the lizard Atlantolacerta andreanskyi(Werner, 1929): A montane cryptic species complex

Atlantolacerta andreanskyi is an enigmatic lacertid lizard that, according to the most recent molecular analyses, belongs to the tribe Eremiadini, family Lacertidae. It is a mountain specialist, restricted to areas above 2400 m of the High Atlas Mountains of Morocco with apparently no connection between the different populations. In order to investigate its phylogeography, 92 specimens of A. andreanskyi were analyzed from eight different populations across the distribution range of the species for up to 1108 base pairs of mitochondrial DNA (12S, ND4 and flanking tRNA-His) and 2585 base pairs of nuclear DNA including five loci (PDC, ACM4, C-MOS, RAG1, MC1R). The results obtained with both concatenated and coalescent approaches and clustering methods, clearly show that all the populations analyzed present a very high level of genetic differentiation for the mitochondrial markers used and are also generally differentiated at the nuclear level. These results indicate that A. andreanskyi is an additional example of a montane species complex.


Background
An emerging pattern among European biotas is that the accentuated environmental instability that occurred during the Pleistocene did not lead to increased speciation rates, with many species and populations originating during the Miocene and proceeding through the Quaternary [1,2]. In many species, population fragmentation was triggered by the beginning of the Messinian Salinity Crisis, a short (600 000 years) but crucial period that occurred between 5.9 and 5.3 Mya during which the Mediterranean Sea desiccated almost completely, producing a general and drastic increase in aridity around the Mediterranean Basin [3,4]. As a result of this increased aridity, forests continued to be replaced by more open and arid landscapes forcing the mesic species to retreat to the moister Atlantic-influenced areas and to the mountainous regions, leading to high speciation in some groups [5,6].
Various studies have attempted to unravel the different roles that the global aridification at the end of the Miocene and the Pleistocene glacial cycles have played in the diversity and distribution of European faunas [7]. However, little is known about the effects that these climatic changes had on species living further South, in the African continent. Recent assessments of central African chameleons have uncovered evidence of long-isolated evolutionary histories, with the survival of palaeoendemics leading to considerable diversity [8]. In general, reptiles are excellent model organisms to assess the relative role that the Pre-Quaternary and Quaternary major climatic events have played in the origin, evolution and distribution of species [9]. Available data from some herpetofauna indicate that a similar pattern to the neighboring Iberian Peninsula exists in North Africa, with deep lineages originating at the end of the Miocene (Chalcides [10], Acanthodactylus [11][12][13], Podarcis [2,14,15], Saurodactylus [16], Ptyodactylus [17], Salamandra [18], Pleurodeles [19]). However, the lack of informative nuclear markers in most of these studies may prevent the recovery of the true evolutionary history of the group [eg. 20,21], and makes it difficult to ascertain if these lineages correspond to species complexes or not. Since there is a strong likelihood of discordance between gene trees and species trees [22][23][24], information from different genetic markers (mitochondrial and nuclear) is thus necessary for delimiting evolutionary lineages, as well as for establishing phylogenetic relationships.
Despite being key concepts in the fields of systematic and evolutionary biology, recognizing and delimiting species are highly controversial issues ([e.g. 25,26]). Recognizing species is not only a taxonomic challenge, but is also essential for other biological disciplines such as biogeography, ecology and evolutionary biology [27], and has serious consequences for conservation biology and the design of effective conservation plans [28,29]. Delimiting species is also the first step towards discussing broader questions on evolution, biogeography, ecology or conservation. Recently, thanks to intellectual progress made in the field with the aim of identifying a common element among all the different species concepts, a single, more general, concept of species known as General Lineage Species Concept has been suggested [30]. This unified species concept emphasizes the common element found in many species concepts, which is that species are separately evolving lineages. Therefore, properties like reciprocal monophyly at one or multiple loci, phenotypic diagnosability, ecological distinctiveness, etc. are not part of the species concept but are used to assess the separation of lineages and to species delimitation [31]. This separation between species conceptualization and species delimitation and the proposal of a unified species concept has concentrated efforts in the development of new approaches for species delimitation, as for example with "integrative taxonomy" [32,33, among others]. Under this new approach, species delineation is regarded as an objective scientific process that results in a taxonomic hypothesis. Therefore, the level of confidence in the taxonomic hypothesis supported by several independent character sets is much higher than for species supported by only one character [34]. Such an integrative view is especially useful in the case of taxonomic groups that are morphologically conservative, where cryptic species have probably been overlooked [17,35,36].
Normally, high altitude species carry signatures of the expansion and contraction cycles occurred during glacial and interglacial periods [37][38][39]. Because of this, they are of particular interest to study historical responses to climate change, since they are adapted to a small window of environmental changes, and usually present low tolerance to high temperatures [40]. In Europe, high altitude species often seem to have persisted through glacial periods by short movements to lower altitudes rather than to the classic "southern refugia" of lowland species. In this way current ranges may primarily reflect postglacial expansions [41]. However, it is not clear if the same phenomenon occurs in African montane taxa.
Atlantolacerta andreanskyi (Werner, 1929) is a lacertid lizard endemic to the western and central parts of the High Atlas Mountains of Morocco. It is restricted to areas above 2400 m [42,43], where it is frequently found in the vicinity of small watercourses or plateaus in the top of the mountains that retain some water from rain or snowmelt. Habitat is normally screes and areas with boulders, meadows and, in particular, the base of cushion-like thorny plants in these places [42; personal observation]. Although A. andreanskyi had initially been placed in several different genera within the subtribe Lacertini [44][45][46][47][48], recent phylogenetic analyses based on mitochondrial DNA and a combination of mitochondrial and nuclear markers [49,50] suggest that A. andreanskyi is a member of the subtribe Eremiadini, and apparently sister to the remaining Eremiadini. This position would conform to this species lacking the synapomorphies that characterize most other Eremiadini, namely a derived condition of the ulnar nerve and the presence of a fully developed armature in the hemipenis, which has folded lobes when retracted. It is also distinctive within the Eremiadini regarding the presence of enlarged masseteric scale [49]. Because of its phylogenetic position, without close relationship to any other genus of Eremiadini and its distinctive morphology it was recently placed in a new monotypic genus, Atlantolacerta [49]. Atlantolacerta andreanskyi is distributed across 440 Km (straight line) of mountainous terrain, with the different populations presenting an apparently disjunct distribution ([42,43; see Figure 1]). As with many montane species, the situation observed in A. andreanskyi is similar to an archipelago, with the different "islands" being represented by mountaintops disconnected due to areas of unsuitable habitat below 2400 m. As a result of this scenario, minimal gene flow is currently expected between the different populations; however, it is not known how the different climatic events occurred during the Miocene and Pleistocene have affected this species. Even though some aspects of the biology of A. andreanskyi are already well known [e.g. 51,52], the genetic structure of the different populations, as well as the relationships between the different populations have never been assessed before.
Therefore, in order to shed some light on the previous questions and attempt to assess the evolutionary history of the species and identify the number of lineages, we sampled the distribution area of the species and performed several combined phylogenetic reconstructions and clustering analyses, using both mtDNA and nuclear markers.
Analyses of the concatenated mtDNA data were mostly congruent (Figure 2A). Seven well-supported lineages were recovered from these analyses (pp > 0. 95 and BP > 70%), corresponding to the populations from J. Awlime, J. Sirwa, Tizin Tichka, J. Azourki, Outabati, J. Ayache, and Oukaimeden and nearby Toubkal that were grouped together. Regarding the relationships among these clades, we could distinguish three main groups, Oukaimeden and Toubkal with J. Sirwa from the southern end of the distribution range; J. Ayache with Outabati from the northern distribution, and Tizin Tichka with J. Azourki from the central distribution range. The population from J. Awlime, from the extreme South of the range, is a genetically distinct lineage related to the northern group, although, both ML and BI analysis weakly support this topology (see Figure 2A).
All the populations present a low level of diversity in the mitochondrial DNA (uncorrected genetic distances 0-0.5% for the ND4 + tRNA-His and 0 -0.2% for the 12S; see Table 1), and a very high level of genetic divergence between populations (5.5 -16.5% in the ND4 + tRNA-His and 2.5 -6.6% in the 12S).

Nuclear genealogies
A total of 77 specimens of A. andreanskyi were sequenced for five nuclear genes. The ACM4 was 447 bp long, presenting 47 haplotypes and 34 polymorphic sites, 33 of them parsimony informative; C-MOS was 534 bp long, with 32 haplotypes and 21 polymorphic sites, all of them parsimony informative; MC1R was 635 bp long, with 57 haplotypes and 36 variable sites, 35 of them parsimony informative; PDC was 441 bp long, with 60 haplotypes and 29 variable sites, 26 of them parsimony informative; RAG1 was 528 bp long, with 38 haplotypes and 19 variable sites, 18 of them parsimony informative. The differences in the genetic distances between the lineages are congruent with the geographic distance between them, supporting the grouping of the lineages in three main groups as seen in the analysis of mitochondrial sequences.
The concatenated analyses of the 5 unphased nuclear markers are congruent with the results obtained in the mitochondrial DNA tree, although with some differences ( Figure 2B). Despite recovering the three main groups observed in the mtDNA analysis, according to the  nuclear markers the J. Awlime population is not sister to the northernmost populations but branches off inside a polytomy with the westernmost lineages at the base of the tree. It is possible to distinguish some of the lineages, although in some cases they are not monophyletic. The J. Ayache population is monophyletic but makes Outabati paraphyletic. The same happens with Tizin Tichka, which makes the population from J. Azourki paraphyletic. The population from Oukaimeden is polyphyletic.

Concatenated analysis (mtDNA and nDNA)
The results of the ML and BI analyses of the mtDNA and nDNA ( Figure 2C) support the same seven lineages as recovered in the mitochondrial analysis, although in this case J. Awlime is sister to the central and northern lineages (Tizin Tichka, J. Azourki, Outabati, and J. Ayache) instead of being sister to only the northernmost lineages ( Figure 2A). As in the mtDNA analysis (Figure 2A), the relationship of J. Awlime with the central and northern lineages is very poorly supported. This result was expected, given the higher resolving power of the mtDNA that contributed with 241 variable sites versus the 150 from the nDNA.

Nuclear networks
As show in Figure 3 and Table 2, there is a moderate degree of haplotype sharing between populations, with most of them lacking private alleles for the nuclear genes analyzed.

Clustering analysis and individual assignment
In our study, the obtained K differs with the combination between the ancestry model and the allele frequency model.

Species tree and divergence time estimates
The results of the clustering analysis with K = 6 were used to define the species for the species tree analysis in STARBEAST. The tree inferred with information from mitochondrial and nuclear markers (phased) (figure 2D) recovered the same topology as in Figure 2C, with all the relationships between the lineages supported by previous analyses. The divergence time estimates were calculated for the six populations (Table 1). High effective sample sizes were observed for all parameters in all BEAST analysis (posterior ESS values > 1000 for all four analyses) and assessment of convergence statistics in Tracer indicated that all analyses had converged. Maximum clade credibility tree for ND4 + tRNA-His was identical in topology to those produced by Bayesian and ML analyses. According to the inferred dates resulted from BEAST ( Figure 2D

Discussion
Extreme mtDNA diversity in A. andreanskyi Several recently published analyses of North African herpetofauna have revealed high levels of endemism and cryptic species [12,14,15,17]. In this analysis, the surprising result was the extreme diversity of mitochondrial DNA found between almost all the populations analyzed The genetic differentiation observed between populations (2.8% -6.6% in 12S and 5.5% -16.5% in ND4 + tRNA-His) is similar and, in some cases, higher than the divergence found between Iberolacerta species (between 7.4% and 8.2% in the cytochrome b gene, [53]), a lacertid genus with most of its species occurring in the mountains of the Iberian Peninsula [41,54]. Initially considered one species, there are now seven recognized species of Iberolacerta in the Iberian Peninsula. Genetic differentiation between these species is lower than between the different populations of A. andreanskyi.
Although the mitochondrial phylogeny supports the existence of seven distinct groups, the clustering analysis only supports the existence of six lineages (J. Sirwa, Tizin Tichka, J. Azourki, Outabati, J. Ayache and a lineage formed by Oukaimeden, Toubkal and J. Awlime). Toubkal samples were always part of the same lineage as Oukaimeden, although, they show some divergence at least at the mitochondrial DNA level (1.7% in ND4 + tRNA-His and 0.3% in 12S). This is not unexpected, as these populations are geographically very close and are part of the High Atlas Mountains, where interconnectivity between populations could occur. The mitochondrial phylogenetic analyses supported the existence of a seventh isolated lineage, J. Awlime, however clustering analysis and the nuclear phylogeny did not support the distinctiveness of this population, possibly because of the small sampling size.
Unfortunately, despite multiple attempts to sample in this remote region, only three individuals were captured. The analyses also could not recover the genetic relationship between J. Awlime and the other populations, because its position in the trees fluctuated between the two main groups (North and South), without support in any of the trees.

Non-reciprocal monophyly in nuclear markers and species delimitation
In the phylogenetic analyses of the concatenated nuclear loci, some of the lineages supported by mtDNA data were not monophyletic. This was observed only between the geographically closest lineages, as in the case of Oukaimeden and J. Sirwa; Tizin Tichka and J. Azourki; and Outabati and J. Ayache, that presumably were in contact more recently than the others. This may be due to the larger effective population size of the nuclear DNA compared to the mitochondrial DNA and the consequent stronger effect of the incomplete lineage sorting at each single nuclear loci [55]. Additionally, the slow evolutionary rate of some of these markers may be a factor. The conjugation of these two effects probably explains the absence of concordance in the single nuclear gene trees (not show), although the same general topology was recovered in the concatenated nuclear phylogeny. Reciprocal monophyly is one of the primary criteria to delimit species [31,56]. Although it is possible to delimit species without observing monophyly in gene trees, since a considerable amount of time must pass after the beginning of divergence of species until they show reciprocal monophyly at a sample of multiple loci [57,58]. Pinho et al. [59] have shown that Podarcis from the Iberian Peninsula and North Africa have a similar pattern (between mtDNA and nuclear) but in a smaller time window and using faster evolving nuclear loci and, in contrast to our case, some populations are in contact. Although we are aware that the determination of K, in STRUCTURE, is only an ad hoc guide to describe consistence between models and the data [60], the program has been commonly used for this end [61]. Several methods based on Bayesian clustering have been developed [62][63][64], however, STRUCTURE is the most widely used, and various studies show its efficiency in assigning individuals to their population of origin [65][66][67][68] and its ability to construct an appropriate clustering hypothesis [61]. However, in the present example the analysis was limited because it was based only in haplotype information. The obtained K differ with the combination model used, but in most of the combinations the analysis supports a K = 6 corresponding to the geographical populations and to the results recovered by the other analyses. This analysis also placed the samples from the J. Awlime population together with the Oukaimeden lineage, possibly due to the limited haplotype sampling. Similarly, the concatenated phylogenetic tree, based on all the genes, supports the existence of 7 lineages giving once more a low support to the relationship between J. Awlime and the other lineages.
The networks of the individual nuclear loci show high percentage of private alleles in some of the lineages, which fluctuate depending on the gene.

Dating the trees
All the lineages are grouped in two main clusters, the northern group composed by J. Ayache, Outabati, J. Azourki and Tizin Tichka; and the southern group that includes Oukaimeden and J. Sirwa. The divergence obtained for these two lineages was around 7.6 Mya, (4.3-11.9), which coincides approximately with the time of the final closing of the Rifian Strait (7.2 Mya; [3]). During the Miocene, tectonic activity in the region was intense and included the uplift of the Atlas Mountains that occurred around 9.0 Mya [69,70]. It was more or less at the same time that Podarcis invaded North Africa (7.5 ± 1.2 Mya, [2]) and the Iberian clade of Iberolacerta started to fragment (6.1 Mya, [1]). The split of the six lineages must have occurred later, probably during the Quaternary Glaciations (4.3 ± 3; 2.4 ± 2; 2.9 ± 2 Mya). However, the confidence intervals obtained were very large, increasing the time window for the events and the associated error. Determination of the time of the speciation events is important to understand the evolutionary biogeography of species [71]. However, it is difficult to estimate ages in phylogenies without several sources of error. Clearly the lineages of A. andreanskyi are pre-Pleistocenic and, as found in Central African chameleons [8] can be considered paleoendemics. However, without better calibration points it is difficult to date the split of the lineages more precisely than this.

Conclusions
Phylogeographic assessments of several taxa in northwest Africa have indicated the presence of cryptic diversity in organisms ranging from scorpions [72] to mammals [73], and reptiles are not an exception [e.g. 11,17,74]. What is exceptional in the case of A. andreanskyi are the high levels of mitochondrial divergence between almost every sampled populations, ranging from 5.5 up to 16.5% (ND4 + tRNA-His) between populations separated by low geographic distances (for example just 60 Km between Oukaimeden and J. Sirwa and 45 Km between Oukaimeden and Tizin Tichka). Six of the eight analyzed populations are highly distinct based on both mtDNA and multiple nuclear markers. This raises the issue not of whether A. andreanskyi is a species complex, but just how many species may occur within the group. Presumably, far more than the six possible species identified in this study, since, probably, many populations remain unsampled. However, preliminary morphological analyses suggest that all the different populations included in the present study are very homogeneous (unpublished data). This may imply the presence of cryptic diversity, but definitive conclusions should wait until a complete morphological study is carried out (work in progress).
Current models of reptiles species accessed for the region indicate low levels of diversity across much of the High Atlas Mountains [75]. Indeed only a few species are recorded at altitudes above 2000 m; typically A. andreanskyi, Quedenfeldtia species (Q. trachyblepharus and Q. moerens), Chalcides montanus and Vipera monticola [e.g. 42,76]. However, the finding of high genetic diversity in A. andreanskyi indicates that unidentified lineages occur, and that the other high mountain species should also be assessed as possible cryptic species candidates. Our results are also essential from a conservation point of view, as many forms may actually have smaller ranges than currently thought, and small isolated populations on high mountains have been identified as those of high concern under typical global warming scenarios [77]. Given these results it is necessary to increase the sampling in order to understand the relationship of J. Awlime with the other populations and try to find new populations. Furthermore it is very important to conduct a through morphological study to determine if there is phenotypic variation, and then to revise the taxonomy of the genus Atlantolacerta.

Species concept and integrative approach
Although the present study does not include a taxonomic revision of the genus Atlantolacerta, like many other works in which some of the authors of the present manuscript have participated [35,78,79], we advocate for the use of the General Lineage Species Concept proposed by de Queiroz [30]. Two lines of evidence have been defined on the basis of alleged independence of their respective datasets: mitochondrial DNA and nuclear DNA. In the present study, we have decided to retain as "putative species" only these lineages that were recovered as monophyletic in the phylogenetic analysis of the mtDNA data and that were supported by the analysis of the nDNA using STRUCTURE v.2.3.2 [60]. Within the framework of an integrative approach, and pending the inclusion of morphological data, this would correspond to Integration by total congruence (ITC). However, it is important to take into account that in the absence of a thorough morphological analysis we do not consider the molecular data presented here enough to revise the taxonomy of the genus Atlantolacerta.

DNA extraction, amplification, and sequencing
A total of 92 individuals from eight different populations distributed across the entire range of Atlantolacerta andreanskyi were sampled for this study: 14 from Oukaimeden, 15 from Tizin Tichka, 14 from Jebel Ayache, 15 from Jebel Azourki, 14 from Outabati, 15 from Jebel Sirwa and 2 from Toubkal and 3 from J. Awlime ( Figure 1 and Table 3). Specimens were caught by hand, identified on the basis of external features, measured and photographed for later morphological studies. Tail tips where collected and stored in 96% ethanol, after which individuals were released in the same place where they were caught.
Genomic DNA was extracted from ethanol-preserved tissue samples using standard high-salt protocols [80]. A total of 89 specimens of Atlantolacerta andreanskyi plus three outgroups (Podarcis hispanica, Podarcis carbonelli and Podarcis bocagei) were sequenced for two mitochondrial regions: partial 12S rRNA (12S) and partial NADH dehydrogenase 4 (ND4) and flanking tRNA (tRNA-His) and 77 specimens for five nuclear gene fragments, recombination-activating gene 1 (RAG1), acetylcholinergic receptor M4 (ACM4), melanocortin receptor 1 (MC1R), oocyte maturation factor Mos (C-MOS) and phosducin (PDC). Primers used for both amplification and sequencing were: 12Sa and 12Sb [81] for the 12S following the PCR conditions described in Harris and Arnold [82], ND4 and Leu for ND4 + tRNA-His, PCR conditions described in Arévalo et al. [83]; L2408 and H2920 for RAG1 following the PCR conditions from Vidal and Hedges [84]; tg-F and tg-R [85] for ACM4 with PCR conditions following Gamble et al. [86]; MC1RF and MC1RR for MC1R following PCR conditions described in Pinho et al. [87]; Lsc1 and Lsc2 for C-MOS following the PCR conditions from Godinho et al. [88]; and PHOF2 and PHOF1 for PDC, following PCR conditions described in Bauer et al. [89]. PCRs were carried out in 25 μl volumes, containing 5.0 μl of 10 reaction Buffer, 2.0 mM of MgCl2, 0.5 mM each dNTP, 0.2 μM each primer, 1 U of Taq DNA polymerase (Invitrogen), and approximately 100 ng of template DNA. Finally, PCR products were purified using exosap IT and the resulting amplified fragments were sequenced on an Applied Biosystem DNA Sequencing Apparatus. Chromatographs were checked manually, assembled and edited using Bioedit 7.0.1 [90]. Sequences were aligned for each gene independently using the online version of MAFFT v.6 [91] with default parameters (gap opening penalty = 1.53, gap extension = 0.0) and FFT-NS-1 algorithm. Coding gene fragments (ND4, C-MOS, ACM4, RAG1, PDC and MC1R) were translated into amino acids and no stop codons were observed, suggesting that the sequences were all functional. Heterozygous individuals were identified based on the presence of two peaks of approximately equal height at a single nucleotide site. SEQPHASE [92] was used to convert the input files, and the software PHASE v2.1.1 to resolve phased haplotypes [93]. Default settings of PHASE were used except for phase probabilities that were set as ≥ 0.7 [94]. All polymorphic sites with a probability of < 0.7 were coded in both alleles with the appropriate IUPAC ambiguity code. Phased nuclear sequences were used for the structure analysis; networks and species tree analysis, and the unphased sequences for the phylogenetic analyses (see below). DnaSP [95] was used to calculate the number of haplotypes (h) and mutations (η). Mega v.3.0 [96] was used to estimate uncorrected p-distances and to obtain the number of variable and parsimony informative sites.

Phylogenetic analyses
Phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian (BI) methods. JModelTest [97] was used to select the most appropriate model of sequence evolution under the Akaike Information Criterion [98]. ML analyses were performed with RAxML v.7.0.4 [99] with 100 random addition replicates. A GTR + I + G model was used and parameters were estimated independently for each partition (by gene). Reliability of the ML tree was assessed by bootstrap analysis [100] including 1000 replications. Bayesian analyses were performed with MrBayes v.3.1.2 [101] with best fitting models applied to each partition by gene and all parameters unlinked across partitions. The models selected for the different partitions were: 12S, GTR + I + G; ND4, GTR + G; tRNA-His, GTR + I + G; ACM4, HKY + I; C-MOS, GTR + I + G; MC1R, HKY + I + G; PDC, GTR + I + G; and RAG1, GTR + I. Two independent runs of 5x10 6 generations were carried out, sampling at intervals of 1000 generations producing 5000 trees. Convergence and appropriate sampling were confirmed examining the    standard deviation of the split frequencies between the two simultaneous runs and the Potential Scale Reduction Factor (PSRF) diagnostic. Burn-in was performed discarding the first 1250 trees of each run (25%) and a majority-rule consensus tree was generated from the remaining trees. In both ML and BI alignment gaps were treated as missing data and the nuclear gene sequences were not phased.

Nuclear Networks
The genealogical relationships between the populations were assessed with haplotype networks for all the individual nuclear genes, constructed using statistical parsimony [102] implemented in the program TCS v 1.21 [103] with a connection limit of 95%. This analysis was made with the phased sequences. Haplotypes were colored taking into account the population of origin.

Population structure -Clustering analyses
A model-based Bayesian clustering method was applied to all haplotypes using STRUCTURE v.2.3.2 [60,104,105].
In this analysis, individuals are probabilistically assigned to either a single cluster (the population of origin), or more than one cluster (if there is admixture). STRUC-TURE was run with haplotype information from the nuclear fragments independently. We ran our data with the all parameters combinations between the Ancestry Model and the Allele Frequency Model to compare the results. The genetic structure was forced to vary from K = 2 to K = 10 clusters, the latter corresponding to the number of geographic populations sampled plus two. STRUCTURE ran for 550 000 steps, of which the first 50 000 were discarded as burn-in. For each value of K ten independent replicates of the Markov Chain Monte Carlo (MCMC) were conducted. To detect the true number of clusters (K) we followed the graphical methods and algorithms outlined in Evanno et al. [61], with the comparison of the average posterior probability values for K (log likelihood; ln L) using the online version, STRUCTURE HARVESTER v0.6.5 (available at: http:// taylor0.biology.ucla.edu/struct_ harvest/, April 2011).

Species tree, and divergence time estimates
Here we applied the coalescent-based species-tree approach implemented in STARBEAST [106] an extension of BEAST v1.6.1 [107] to test the origin and diversification patterns in Atlantolacerta, and to compare these results to those obtained from the ML and BI analyses of the concatenated dataset. This analysis needs a priori information regarding the species/populations delimitation and the species/populations assignation of the individuals in order to reconstruct the topology of the species tree. For this approach, we used the results obtained from previous clustering analyses to define the groups of individuals to be used as "species" (populations) in STARBEAST [106]. The clustering analysis supported the existence of six lineages, as Oukaimeden, Toubkal and J. Awlime were included in the same lineage. All five nuclear gene fragments, 12S and the fragment consistent of the ND4 and flanking tRNA-His were included in the analyses as 7 independent partitions. The phased dataset was used for the nuclear loci.
The input file was formatted with the BEAUti utility included in the software package. We performed two independent runs of 1.5 x 10 8 generations, sampling every 15 000 generations, from which 10% were discarded as burn-in. Models and prior specifications applied were as follows (otherwise by default): 12S -GTR + G; ND4 and tRNA-His -HKY + G; MC1R -HKY + I; ACM4 -HKY + I; C-MOS -GTR + I + G; RAG1 -HKY + I; PDC -GTR + I; Relaxed Uncorrelated Lognormal Clock (estimate); Yule process of speciation; random starting tree; alpha Uniform (0, 10).
For all analyses implemented in BEAST, convergence for all model parameters was assessed by examining trace plots and histograms in Tracer v1.5 [108] after obtaining an effective sample size (ESS) > 200. The initial 10% of samples were discarded as burn-in. Runs were combined using LogCombiner, and maximum credibility trees with divergence time means and 95% highest probability densities (HPDs) were produced using Tree Annotator (both part of the BEAST package). Trees were visualized using the software FigTree v1.3.1 [109].
Several studies have already calculated divergence rates for reptiles, and particularly for lacertids [2,15,49]. Pinho et al. [15] used well-known and dated independent geological events in the Aegean [110] to estimate a maximum and minimum mutation rate for the ND4 mitochondrial fragment (and flanking tRNA-His) for the lacertid lizards of the genus Podarcis (0.0278 and 0.0174 mutation/site/million years, respectively). However, this was the only information available for our data, since we did not have any fossils or calibrations for nuclear markers. It is important to bear in mind that, in the absence of accurate calibration points in the phylogeny from external and independent data (fossil records, known biogeographic events, or paleoclimatic reconstructions) or as a result of the heterogeneity in the evolutionary rate between the calibrated and uncalibrated taxa, temporal estimates by means of molecular data could be a potential source of inference error, and, therefore, they should be treated with caution [111]. Despite the limitations of molecular clocks [111,112], divergence time estimates can still provide a proxy for the temporal window of evolutionary diversification in species groups of interest. Therefore and taking into account our data limitations and availability, we used BEAST v.1.6.1 [107] to estimate dates of the cladogenetic events using only ND4 and flanking tRNA-His. We used a phylogeny pruned arbitrarily to include one representative from each of the major lineages uncovered with the concatenated analysis (6 specimens in total, we excluded J. Awlime population, because of the lack of support of the branch in previous analyses). This method excludes closely related terminal taxa because the Yule tree prior (see below) does not include a model of coalescence, which can complicate rate estimation for closely related sequences [113]. Analyses were run four times for 5x10 7 generations with a sampling frequency of 10 000. Models and prior specifications applied were as follows (otherwise by default): GTR + G for 12S; HKY + G for ND4 and tRNA-His; HKY + I for MC1R; HKY + I for ACM4; GTR + G + I for C-MOS; HKY + I for RAG1; GTR + I for PDC; Relaxed Uncorrelated Lognormal Clock (estimate); Yule process of speciation; random starting tree; alpha Uniform (0, 10); ucld.mean of ND4 Normal (initial value: 0.0226, mean: 0.0226, Stdev: 0.0031).
Authors' contributions MB carried out the molecular laboratory work, analyzed the data and drafted a preliminary version of the manuscript. All authors participated in the conception and design of the study, collection of samples, writing and approval of the final manuscript.