- Research article
- Open Access
Integrative analyses of speciation and divergence in Psammodromus hispanicus(Squamata: Lacertidae)
© Fitze et al; licensee BioMed Central Ltd. 2011
- Received: 9 May 2011
- Accepted: 30 November 2011
- Published: 30 November 2011
Genetic, phenotypic and ecological divergence within a lineage is the result of past and ongoing evolutionary processes, which lead ultimately to diversification and speciation. Integrative analyses allow linking diversification to geological, climatic, and ecological events, and thus disentangling the relative importance of different evolutionary drivers in generating and maintaining current species richness.
Here, we use phylogenetic, phenotypic, geographic, and environmental data to investigate diversification in the Spanish sand racer (Psammodromus hispanicus). Phylogenetic, molecular clock dating, and phenotypic analyses show that P. hispanicus consists of three lineages. One lineage from Western Spain diverged 8.3 (2.9-14.7) Mya from the ancestor of Psammodromus hispanicus edwardsianus and P. hispanicus hispanicus Central lineage. The latter diverged 4.8 (1.5-8.7) Mya. Molecular clock dating, together with population genetic analyses, indicate that the three lineages experienced northward range expansions from southern Iberian refugia during Pleistocene glacial periods. Ecological niche modelling shows that suitable habitat of the Western lineage and P. h. edwardsianus overlap over vast areas, but that a barrier may hinder dispersal and genetic mixing of populations of both lineages. P. h. hispanicus Central lineage inhabits an ecological niche that overlaps marginally with the other two lineages.
Our results provide evidence for divergence in allopatry and niche conservatism between the Western lineage and the ancestor of P. h. edwardsianus and P. h. hispanicus Central lineage, whereas they suggest that niche divergence is involved in the origin of the latter two lineages. Both processes were temporally separated and may be responsible for the here documented genetic and phenotypic diversity of P. hispanicus. The temporal pattern is in line with those proposed for other animal lineages. It suggests that geographic isolation and vicariance played an important role in the early diversification of the group, and that lineage diversification was further amplified through ecological divergence.
- Enhance Vegetation Index
- Niche Divergence
- Niche Conservatism
- Messinian Salinity Crisis
- Nuptial Coloration
Species diversity emerges from the combination of both past and ongoing evolutionary and ecological processes driving speciation [1–3]. However, it is challenging to determine the relative contributions of historical and ecological factors in causing genetic differentiation . The traditional classification of modes of speciation (allopatric, peripatric, parapatric, and sympatric) within a spatial context [5, 6] is currently revisited in the light of recent studies that integrate phylogenetic, ecological, and geographical data [3, 7, 8]. In the last decade, evolutionary biologists have focused on discerning the mechanisms leading to reproductive isolation, and the field has witnessed major advances in determining the relative contribution of historical geographic barriers to diversification thanks to the possibility of linking geological and phylogenetic data . In contrast, the elucidation of the contribution to diversification of ecologically-based divergent selection due to environmental differences has been hindered until recently by the difficulty of simultaneously gathering genetic, phenotypic, and environmental data within the framework of a single study . However, the recent accumulation of environmental data and the development of ecological niche modelling allow overcoming these limitations and provide the basis for an integrative approach that combines phylogenetic and biogeographic data in order to explain the origin and large-scale distribution patterns of biodiversity . Integrative analyses provide new insights on the factors driving diversification and speciation, and allow disentangling the effects of environment from those of historical barriers [8–12]. In particular, it is possible to test explicitly whether diversity can evolve in allopatry and under similar ecological conditions (i.e. niche conservatism; ) or whether different ecological environments (i.e. niche divergence) that promote divergent natural selection are at the root of diversification [13, 14].
Here, we integrate phylogenetic, phenotypic, geographic, and environmental data to investigate the contributions of historical geographic barriers and environmental differences to speciation and divergence in the Spanish sand racer (Psammodromus hispanicus). The current distribution of P. hispanicus includes the Iberian Peninsula, and the French Mediterranean coast with an upper altitudinal limit at 1700 m a.s.l. . The broad distribution of this group, which inhabits regions with very distinct habitats as well as areas with complex geological histories, makes it a suitable model to investigate how vicariant events due to geographical barriers and niche divergence due to selection have influenced diversification. P. hispanicus Fitzinger, 1826 consists of two subspecies, namely P. hispanicus hispanicus Fitzinger, 1826 and P. hispanicus edwardsianus (Dugès, 1829). The Iberian Peninsula hosts a second species of the same genus, P. algirus (Linnaeus, 1758), consisting of two divergent Eastern (E) and Western lineages, the latter including African, as well as Northwestern (NW) and Southwestern (SW) Iberian clades . Other species of the genus Psammodromus are P. blanci (Lataste, 1880) from Algeria, Morocco, and Tunisia, and P. microdactylus (Boettger, 1881) endemic to Morocco.
Based on a representative sampling of P. hispanicus in Spain, we first reconstructed phylogenetic relationships among sampled populations using both mitochondrial (mt) and nuclear markers, and dated major cladogenetic events, comparing patterns observed in P. hispanicus with those of P. algirus. Second, we investigated differences between molecular lineages in phenotypic traits using multivariate analyses. Third, we performed ecological niche modelling and applied different procedures to assess niche divergence and the spatial structure of shared environmental conditions among lineages, in order to investigated the evolutionary and ecological processes that promoted genetic and phenotypic differentiation in the Spanish sand racer. More specifically, we tested whether niche divergence and/or allopatric speciation may explain the observed diversity. Under the niche divergence hypothesis we predicted that two closely related taxonomic groups would live in habitats characterized by different environmental conditions. Under the allopatric speciation hypothesis we predicted genetic, but not necessarily ecological divergence, and thus that environmental niches of sister species should be more similar than under ecological speciation. Finally, we speculate which geological events may have led to the observed diversity.
Phylogenetic Relationships within Psammodromus
Phylogenetic Relationships and Phylogeography of P. hispanicus
Molecular differentiations in mtDNA between northern and southern groups of a) P. hispanicus edwardsianus and b) P. hispanicus hispanicus Western lineage.
a) P. hispanicus edwardsianus
among populations within groups
b) P. hispanicus hispanicus Western lineage
among populations within groups
Testing northward range expansion
Phenotypic Differences within P. hispanicus
Phenotypic differences between the three P. hispanicus lineages
b. Means ± SE per lineage
c. Univariate ANOVAs
Femoral pores (#)
12.1 ± 0.1
9.9 ± 0.1
11.2 ± 0.2
F 2,208 = 102.03
Throat scales (#)
20.4 ± 0.2
17.9 ± 0.2
18.9 ± 0.3
F 2,208 = 40.55
1.9 ± 0.1
0.7 ± 0.1
2.0 ± 0.4
F 2,208 = 31.71
2.492 ± 0.038
2.199 ± 0.039
2.337 ± 0.053
F 2,208 = 12.46
1.065 ± 0.005
1.043 ± 0.006
1.098 ± 0.009
F 2,208 = 12.36
Anal scale width (mm)
0.063 ± 0.001
0.061 ± 0.001
0.066 ± 0.001
F 2,208 = 3.77
Body mass (g)
1.693 ± 0.03
1.877 ± 0.04
1.866 ± 0.07
F 2,208 = 7.12
Ventral scales (#)
24.6 ± 0.2
25.8 ± 0.3
25.1 ± 0.4
F 2,208 = 6.55
0.482 ± 0.003
0.500 ± 0.002
0.498 ± 0.004
F 2,208 = 7.49
46.39 ± 0.33
49.81 ± 0.43
48.31 ± 0.72
F 2,208 = 18.04
Collar scales (#)
0.3 ± 0.1
1.3 ± 0.2
0.7 ± 0.2
Ha2,208 = 29.25
0.5 ± 0.1
4.3 ± 0.3
5.2 ± 0.5
Ha2,208 = 106.59
d . Post - hoc analyses
edwardsianus vs Central
edwardsianus vs Western
Central vs Western
Femoral pores (#)
t170 = 9.13, P < 0.0001
t145 = 5.09, P < 0.0001
t101 = 6.82, P < 0.0001
Throat scales (#)
t170 = 9.03, P < 0.0001
t145 = 4.38, P = 0.0002
t101 = 2.92, P = 0.038
t170 = 7.60, P < 0.0001
t145 = 0.22, P = 0.830
t101 = 6.29, P < 0.0001
t170 = 5.20, P < 0.0001
t145 = 2.33, P = 0.105
t101 = 1.89, P = 0.366
t170 = 2.53, P = 0.024
t145 = 3.28, P = 0.012
t101 = 5.13, P < 0.0001
Anal scale width (mm)
t170 = 1.60, P = 0.111
t145 = 1.62, P = 0.429
t101 = 2.76, P = 0.055
Body mass (g)
t170 = 3.51, P = 0.002
t145 = 2.46, P = 0.105
t101 = 0.47, P = 1.000
Ventral scales (#)
t170 = 3.63, P = 0.002
t145 = 1.24, P = 0.537
t101 = 1.68, P = 0.483
t170 = 3.29, P = 0.004
t145 = 3.01, P = 0.025
t101 = 0.21, P = 1.000
t170 = 5.76, P < 0.0001
t145 = 2.43, P = 0.105
t101 = 1.53, P = 0.520
Collar scales (#)
TcBF = 5.50, P < 0.0002
TcBF = 1.79, P = 0.537
TbBF = -2.44, P = 0.322
TcBF = 16.3, P < 0.0003
TcBF = 9.18, P < 0.0001
TbBF = 1.82, P = 0.537
Results from univariate ANOVAs are shown in Table 3b and 3c. In brief, there were statistically significant differences between all three major lineages in the number of femoral pores, in the number of throat scales and in the snout shape. The number of ocelli differed between the Central lineage and the other two lineages, but no differences were present between edwardsianus and the Western lineage. The snout-to-vent length (SVL), SVL ratio, body mass, the number of ventral scales, and the number of collar scales differed between edwardsianus and the Central lineage, and there were no differences between the Western lineage and the two other lineages. There were statistically significant differences between the edwardsianus lineage and the other two lineages in head ratio and nuptial coloration, and no differences between Central and Western lineage.
Ecological Niche Modelling within P. hispanicus
Five out of eight sampled populations (62.5%) that were previously unknown to the authors were located in 10 × 10 km squares where P. hispanicus has not been recorded previously . This indicates that the distribution of P. hispanicus is underestimated, and thus biogeographic modelling may be importantly hindered when using presence/absence data. In this regard, the Spanish Atlas of Amphibians and Reptiles  does not discriminate between presences belonging to the different lineages. Consequently, we used modelling techniques that do not require absence data to link our presence records with environmental predictors and run Ecological Niche Factor Analysis (ENFA) based on a systematic population sampling (1 population per intersection of a ± 250 km grid) that included 22 populations (Figure 4, Additional File 1: Table S1).
We evaluated the model error of the niche models using an independent data-set obtained from the Official Spanish Atlas of Amphibians and Reptiles . The model omission error was 13.7% with respect to the currently known presences of P. hispanicus, which is reasonably small. We further determined the importance of the sampling points for the lineage predictions using jackknife methodology, in order to investigate whether sampling points of a given lineage may be more important than in other lineages and thus whether differing sample sizes among lineages may have biased our predictions and conclusions. The average proportional habitat suitability (HS) difference was similar among lineages (one-way ANOVA: F2,19 = 0.91, P = 0.42), indicating similar population representativeness among lineages and no bias due to unequal sample sizes. These are important conditions for comparing niche divergences/similarities and overlaps among pair combinations of lineages.
Assessment of predictor relevance
Differences between the three lineages in environmental population parameters
Estimates (mean ± SE)
Annual mean temperature
5810.1 ± 297.3
662.6 ± 394.9
-928.9 ± 458.2
266.4 ± 405.4
Max temperature of warmest month
3256.5 ± 78.5
9.956 ± 104.246
21.378 ± 120.969
-31.3 ± 107.0
Mean temperature of driest quarter
4945.4 ± 327.7
-26.469 ± 438.302
-213.645 ± 505.132
240.113 ± 446.93
Mean temperature of wettest quarter
168.728 ± 7.016
22.606 ± 9.319
-4.128 ± 10.814
-18.478 ± 9.568
Min temperature of coldest month
136.538 ± 5.775
12.344 ± 7.671
-20.716 ± 8.902
8.372 ± 7.876
Annual temperature range
1156.6 ± 52.3
-88.0 ± 69.5
174.1 ± 80.7
-86.1 ± 71.4
258.175 ± 4.548
-5.177 ± 6.041
-0.716 ± 7.010
5.892 ± 6.203
5.901 ± 0.054
-0.112 ± 0.071
-0.038 ± 0.083
0.150 ± 0.073
Precipitation of coldest quarter
7.253 ± 0.147
-0.403 ± 0.195
-0.302 ± 0.226
0.705 ± 0.200
Precipitation of warmest quarter
7.155 ± 0.420
0.014 ± 0.558
0.499 ± 0.648
-0.513 ± 0.573
10.283 ± 0.343
0.223 ± 0.456
-1.329 ± 0.529
1.106 ± 0.468
14.077 ± 0.062
-0.171 ± 0.082
-0.082 ± 0.096
0.252 ± 0.085
70.177 ± 7.974
-15.156 ± 10.592
24.713 ± 12.291
-9.557 ± 10.875
Predictive maps of habitat suitability and suitability overlaps
Ecological Niche Divergence within P. hispanicus
Estimation of niche divergence using different comparisons of lumped and split models
false positive rates
edwardsianus - Central
Western - Central
edwardsianus - Western
Finally, our results of ensemble predictions showed that ENFA results of predictive maps, niche divergences and geographic overlaps are robust to inter-model variability arising from different algorithms used for model building (Additional File 2).
Here, we address how geology, climate, and ecology shaped current diversity in Psammodromus hispanicus using a multidisciplinary approach including phylogenetic, phenotypic, phylogeographic, and ecological niche analyses.
Species Status and Phylogenetic Relationships within Iberian Psammodromus
We reconstructed congruent mt- and nuclear-based phylogenies (and the corresponding networks) that confidently recovered three major lineages in P. hispanicus, corresponding to P. hispanicus edwardsianus, P. hispanicus hispanicus Central lineage, and P. hispanicus hispanicus Western lineage. The molecular clock indicates that the age of divergence of the Central and the edwardsianus lineages was about 4.8 (1.5-8.7) Mya, which, together with phylogeographic and phenotypic evidence, strongly suggest that these two lineages reflect independent evolutionary units. The divergence between these two lineages and the Western lineage was estimated to have occurred about 8.3 (2.9-14.7) Mya. Altogether, the rather old age of divergence, the lack of haplotype and geographic overlap, as well as the existence of phenotypic differentiation, and ecological niche divergence (in the Central lineage), allow postulating that both the Central and the Western lineage may be valid species.
There is significant phenotypic differentiation among the three lineages. The Western lineage was phenotypically intermediate between the edwardsianus and Central lineages, which showed higher phenotypic differentiation. The latter two lineages differed in 11 of the 12 studied traits (Table 3: femoral pores, number of ocelli, SVL, SVL ratio, snout shape, body mass, number of ventral scales, head ratio, number of collar and throat scales, and nuptial coloration) whereas the Western lineage differed from the edwardsianus lineage in five traits (femoral pores, number of throat scales, snout shape, head ratio, and nuptial coloration) and from the Central lineage in four traits (femoral pores, number of ocelli and throat scales, and snout shape). The Western lineage showed trait values that were intermediate between those measured in the Central and the edwardsianus lineages in 8 of the 12 measured traits (Table 3a). Moreover, all specimens of the edwardsianus lineage could be distinguished from the other two lineages by the presence of a supralabial scale below the subocular scale. This finding is in line with previous taxonomy, where its presence has been used to distinguish between the two subspecies, P. hispanicus hispanicus and P. hispanicus edwardsianus .
For P. algirus, phylogenetic analyses based on molecular data provide evidence for the existence of four statistically supported lineages, two corresponding to the SW and NW clades of P. algirus , another corresponding to the nominate lineage of P. algirus (from Africa), and the oldest one corresponding to the Eastern lineage of P. algirus . The molecular clock showed that the NW and SW clades split only around 1.00 (0.3 - 1.9) Mya, and that the Eastern lineage split around 3.01 (1.0 - 5.5) Mya. The former divergence estimate is similar to the oldest divergence times inferred for the clades within the edwardsianus and Central lineage (0.86 (0.3 -1.6) and 0.76 (0.2 - 1.5) Mya, respectively; Figure 3), whereas the latter estimate coincided with the split of the clades within the Western lineage (2.80 (0.9 - 5.1) Mya; Figure 3). The younger datings correspond to the Pleistocene, suggesting a role of glaciations at the origin of the different lineages. Interestingly, there is geographic overlap between northern and southern populations of the Central and the edwardsianus lineage, as well as between the SW and NW clades of P. algirus, suggesting that in both cases reproductive isolation among lineages/clades may exist. Within lineages, the recovered trees and network, as well as population genetic analyses suggest northward expansion of P. hispanicus from southern refugia [21, 22], and incipient and ongoing genetic isolation, while the pattern is less clear within P. algirus.
According to our results, the P. algirus group likely had an Iberian origin, and first split into eastern and western lineages (Figure 1). Within the western lineage, a second split separated the Iberian ancestor of the Northwestern + Southwestern clades from Moroccan P. algirus, suggesting that P. algirus colonized Africa from the Iberian Peninsula. These inferred biogeographic patterns confirm previously reported results based on partial 12S rRNA, 16S rRNA and cytochrome b sequences . However, our molecular dating suggests that the Eastern lineage split from the other lineages slightly later than previously estimated (3.6 ± 0.05 mya ). Unfortunately, we cannot determine whether the origin of the genus Psammodromus was Iberian or North African since our phylogeny did not include two key African species, P. blanci (likely the sister group of P. hispanicus; ) and P. microdactylus.
Differentiation of the Western lineage from the ancestor of the other two P. hispanicus lineages occurred in the Miocene when a progressive uplift started to close the East of the Betic Straits, and formed the Guadalquivir basin (Early Messinian; 7.2-5.5 Mya) . During the same geological period, there were hypothesized splits in several other Iberian reptile and amphibian genera such as Lissotriton , Alytes, , and Blanus , producing in some of them [25 and 27] a similar east-west differentiation pattern. The split between the Central and the edwardsianus lineages dates back to the Miocene/Pliocene boundary and thus close to the Messinian salinity crisis and the opening of the Gibraltar Strait. During this period an uplift of the Spanish Central System occurred that led to the current configuration of the Iberian Peninsula's main river drainages , indicating that major geologic and climatic changes occurred. Accordingly, Pliocene diversification has been reported for many Iberian groups including freshwater fishes [29–31] and amphibians [32–34].
The spatial distribution of mt cytb, nad4, nuclear suppressor of SWI4 1, and nuclear clone 17 diversity showed current allopatry for all three P. hispanicus lineages suggesting a vicariant event at their origin (see above). A decrease in cytb diversity with increasing latitude was observed, which likely indicates northward range expansion of all three lineages. The large, negative, and significant test statistics of the neutrality tests in the edwardsianus and Central lineage further supported this result. However, statistical support for a northward range expansion in the Western lineage was low, most likely due to the small sample sizes obtained from the southern Peninsula (N = 7 individuals belonging to three different haplotypes) and the rather homogenous haplotype distribution on the northern Peninsula (mainly haplotype CE7). Similar patterns were observed in the suppressor of SWI4 1 diversity. The observed range expansions may be the result of post-glacial range expansions  from glacial refugia located south of the Guadiana and the Júcar River. According to the molecular clock, range expansions may have occurred around 0.8 (0.2 - 1.5) Mya, which coincides with the Pleistocene glacial and interglacial periods . Similar Pleistocene patterns of range expansion during interglacial periods have been reported for reptile species [37, 38] and also for insects, molluscs, amphibians, mammals, and plants [39, 40].
Niche Modelling within P. hispanicus
We performed ecological niche modelling for each of the three identified lineages to obtain predictive maps of habitat suitability, and assess predictor relevance and niche divergence. The habitat suitability maps fitted the realized niche reasonably well, and the model omission error was only 13.7% with respect to the P. hispanicus presences cited in the official Spanish Atlas of Amphibians and Reptiles  (analysis based on 10 × 10 km resolution) and no sampling bias could be detected.
For all comparisons among the three P. hispanicus lineages, habitat suitability scores in the sampled populations were significantly different among lineages and highest in the lineage populations for which the distribution was modelled. These results indicate that lineages tend to differ in the optimal portions of their niches. According to the spatial prediction models, the mean temperature of driest quarter was one of the most important predictors of the observed distributions of the Central and the Western lineage. In contrast, mean temperature of wettest quarter, minimum temperature, and precipitation of coldest quarter were the most important predictors of the edwardsianus lineage's distribution. This pattern was in line with the finding that the edwardsianus lineage inhabited areas with lower vegetation cover on the generally warmer and drier eastern coast of the Iberian Peninsula. The Central lineage inhabited central peninsular habitats characterized by intermediate vegetation cover, precipitation and temperature in wettest quarter, and with the lowest minimum temperatures of coldest quarter and precipitation seasonality. The Western lineage lived in habitats with the highest vegetation cover, winter precipitation, precipitation seasonality, and winter minimum temperatures, and with the lowest temperatures during the wettest quarter, corresponding to the more humid and climatically more stable Western parts of the Iberian Peninsula.
We used different approaches to assess niche divergence between the three lineages, and analyses revealed the same overall pattern. Interpredictivity, differential model overprediction among hierarchical taxonomic groups, and the extent of geographic overlap of the model predictions, showed that the ecological niche of the Central lineage was most divergent, whereas ecological niches of the Western and edwardsianus lineages were more similar. Based on the reconstructed tree topologies and following the principle of parsimony, we can infer that the ancestor of the Central and edwardsianus lineages likely occupied a niche similar to that of its sister group, i.e. the Western lineage. Hence, niche divergence occurred during the evolution of the Central lineage.
There is suitable climatically suitable habitat for the edwardsianus lineage in the west of the Guadalquivir River, and suitable habitat was predicted for the Western lineage on the Eastern Iberian Peninsula (Figures 8a and 8c). The spatial predictions showed important overlap between the edwardsianus and the Western lineage on the southern and southwestern Iberian Peninsula (Figure 8e), where both lineages share potentially suitable ecological conditions on both sides of the Guadalquivir River. This suggests that a barrier between Málaga and the Guadalquivir River may prevent population mixing and led to vicariant diversification by impeding dispersal, from the betic uplift in the late Miocene/early Pliocene until the present. Earlier findings in amphibians (e.g. Discoglossus galganoi and D. jeanneae ) are in agreement with this hypothesis.
When comparing the Central lineage with the edwardsianus lineage, suitable habitat for both lineages was located in the centre of the Iberian Peninsula (Figure 8d), and there was no connection through HS values in the estimated contact zone (Figure 8a and 8b). In contrast, the Central lineage showed almost no habitat overlap with the Western lineage (and thus, the Central-edwardsianus ancestor) (Figure 8f), which indicates that niche divergence may be an important force preventing the mixing of these lineages.
In summary, these results show that the Western lineage and the ancestor of the Central and edwardsianus lineages may have been geographically isolated due to a barrier  that still may prevent mixing of the Western and edwardsianus lineages, and that niche divergence may have played a limited role in the separation of these two lineages. In contrast, our analyses provide evidence that niche divergence was more prominent in the Central lineage, potentially preventing gene flow with its sister lineage.
Integrating Phylogenetic, Phenotypic, Geological, and Environmental Data
Determining the relative role of historical and ecological factors as evolutionary drivers of diversification is a central question in evolutionary biology. In this work, we performed a multidisciplinary approach to delimit current diversity of P. hispanicus, and to understand its origin and maintenance. Phylogenetic (mt, nuclear, and combined) and phenotypic data allowed us to differentiate three lineages, which showed important differences in phenotypic traits. The early splitting of the Western lineage may coincide with a vicariant event, which is the initiation of the betic uplift and the formation of the Guadalquivir basin at the end of the Tortonian about 7 Mya, but the large confidence intervals hinder more precise dating. Interestingly, ecological niche modelling shows large overlap in suitable habitat between the Western and the edwardsianus lineage, which contrasts with the phylogeographic cytb, nad4, suppressor of SWI4 1, and clone 17 haplotype distribution that show no spatial overlap of the two lineages. The molecular clock together with genetic, geological, and niche-modelling data suggest that an event related to the betic uplift at the Miocene-Pliocene boundary hindered until present gene flow of Psammodromus. The important niche overlap together with the limited gene flow supports a diversification model of niche conservatism in allopatry for the early divergence of the Western lineage and the common ancestor of the Central and the edwardsianus lineage.
The split between the Central and the edwardsianus lineage in the Early Pliocene coincides with the uplift of the Spanish Central System, which resulted in a change of the drainage patterns from internal to external, and the forming of the present river systems in the Iberian Peninsula . These geological changes produced major climate changes ranging from dry climate during the Messinian salinity crisis to more humid habitats in the early Pliocene. These climate changes may be responsible for diversification within P. hispanicus. The finding that the ecological niche of the Central lineage was most divergent with respect to those of the edwardsianus and Western lineages (despite the older age of the Western lineage split), implies that niche divergence and ecologically-based divergent selection were involved in the diversification process. However, since the Central lineage diverged from the edwardsianus lineage when major climatic changes happened, the possibility that a climatic barrier and initial niche conservatism could have been responsible for the initial splitting cannot be fully discarded .
The large overlap of suitable habitat between the edwardsianus and the Western lineage compared to the small overlap with the Central lineage should be reflected in those phenotypic traits that are adaptive. We explored two phenotypic traits that may be under natural selection in lizards (coloration and number of throat scales). Vegetation cover may determine which colours are cryptic and which ones are conspicuous, and thus, background matching to avoid predation may be the cause for the evolution of colour differences [44, 45]. The Western and the Central lineages do not differ in their nuptial coloration, which is greener than that of the edwardsianus lineage. A second example is the number of throat scales. Since smaller scales and more numerous scales reduce skin water exchange , differences in the number of throat scales may have evolved due to precipitation differences. Here we found that the edwardsianus lineage shows an increased number of throat scales with respect to the Western lineage, whereas the Central lineage shows the lowest number of throat scales. When comparing the ecological scenarios derived from phenotypic traits (Table 3) with real differences in environmental parameters (Table 4), there was a general lack of correlation between environmental parameters and known phenotypic traits under selective pressure, which is in line with previous findings that phenotype does not necessarily predict ecology. Thus, our results suggest that the use of phenotypic traits as a surrogate for ecology in studies dealing with phylogenetic niche conservatism may be problematic , but see . The complexity in phenotypic variation found here, encourages future studies that aim at partitioning phenotypic variation into independent contributions of ecology, phylogenetic inertia, and phylogenetically structured ecological variation, as proposed for higher taxon levels .
Our results indicate that divergence due to both historical geographic barriers and environmental differences may have led through time to the evolution of three P. hispanicus lineages, and that these processes are still acting now to prevent population mixing over the largest part of their allopatric distributions. Our preliminary results on the phylogeographic patterns observed in P. algirus suggest that similar patterns may also exist in other related lizard species, as previously suggested. Here, we highlight the importance of taking a multidisciplinary approach for disentangling the relative roles of vicariant and adaptive divergence in generating currently observed biological diversity. We found that a vicariant event was at the origin of the first splitting event, which was followed by a second splitting event (split between P. hispanicus hispanicus Central and P. hispanicus edwardsianus lineage) in which the role of ecologically based divergent selection (i.e. niche divergence) may have been more prominent. This indicates that diversity in a lineage is the result of different temporally separated evolutionary processes, which is concordant with patterns observed in other groups (e.g. cichlids, ; Gobiidae, ; Passerine birds, ; anoles, [52, 53]).
Samples and DNA Extraction
We conducted a systematic population sampling all over Spain (one population per intersection of a ± 250 km grid) and captured a total of 285 specimens between April and May 2006. The capture and handling of lizards was conducted under the licenses provided by Junta de Andalucía, Gobierno de Aragón, Junta de Castilla y León, Junta de Comunidades de Castilla - La Mancha, Generalitat de Catalunya, Junta de Extremadura, Xunta de Galicia, Comunidad de Madrid, Gobierno de Navarra, Generalitat Valenciana, Parque Natural de l'Albufera (Valencia), Parque Natural del Delta del Ebro (Cataluña), Parque Nacional de Doñana (Huelva), and Gobierno de España. Of the 285 captured specimens, 265 were identified as members of P. hispanicus, whereas 20 specimens were identified as members of P. algirus. We were able to collect individuals in 11 previously known populations (50%; see Additional File 1), failed to find any sample in three known populations, but found samples in adjacent unknown populations (13.6%; ). We also screened potential habitats in locations where no records existed previously, and were able to successfully collect specimens in eight yet unknown populations (36.4%). Sample locations are shown in Figure 4, and the location, sample size and collection numbers are given in Additional File 4.
For each captured individual, we collected a small piece of the tail tip, which was preserved in 70% ethanol at -20ºC. Genomic DNA was extracted from ethanol-preserved tissues using the ChargeSwitch gDNA Micro Tissue Kit (Invitrogen). All individuals were sequenced for mt cytochrome b (cytb) gene, and a subsample including 56 representative specimens of P. hispanicus and 16 specimens of P. algirus was also sequenced for mt NADH dehydrogenase subunit 4 (nad4) gene, and two nuclear loci (see Additional File 4 for specimen number and GenBank accession numbers).
List of newly developed primers used for PCR and DNA Sequencing
suppressor of SWI4 1
suppressor of SWI4 1
PCR amplifications were conducted in 25 µl reactions containing 67 mM Tris- HCl, pH 8.3, 1.5 mM MgCl2, 0.4 mM of each dNTP, 2.5 µM of each primer, template mtDNA (10-100 ng), and Taq DNA polymerase (1.5 U, Roche). For the PCR amplification of the cytb gene fragment, an initial 60 s denaturing step at 93º C was followed by 35 cycles of denaturing at 93º C for 60 s, annealing at 45-50º C for 60 s, and an extension phase at 72º C for 60 s. The final extension phase at 72º C lasted for 6 min. PCR cycling conditions for amplifying the nad4 fragment were: cycle 1 (94º C for 60 s, 72° C for 60 s), cycle 2 - 36 (94º C for 60 s, 50-58º C for 60 s, 72º C for 60 s), cycle 37 (72º C for 6 min). PCR cycling conditions for amplifying partial suppressor of SWI4 1 gene were: cycle 1 (94º C for 60 s, 72° C for 60 s), cycle 2 - 41 (94º C for 60 s, 53.5-59.5ºC for 60 s, 72º C for 60 s), and cycle 42 (72° C for 6 min). Those for amplifying locus 17 were: cycle 1 (94º C for 60 s, 72° C for 60 s), cycle 2 - 41 (94º C for 60 s, 53-66º C for 60 s, 72º C for 60 s), and cycle 42 (72° C for 6 min).
PCR products were checked in 1.5% agarose gels, purified by standard ethanol precipitation, and sequenced in an automated DNA sequencer (ABI PRISM 3700, Applied Biosystems) with the corresponding PCR primers using the BigDye Terminator v3.1 Cycle Sequencing Kit, and following manufacturer's instructions.
Three sequence data sets were analyzed: the first data set (mt data set) included cytb gene partial sequences of P. hispanicus (285 specimens), and of seven individuals of P. algirus (accession numbers: DQ150367, DQ150366, DQ150365, DQ150364, DQ150363, DQ150362, AF206535) obtained from GenBank. In addition, cytb sequences of three individuals of Gallotia caesaris caesaris (AY151843, AY154903, AF439948) and three individuals of Gallotia caesaris gomerae (AY151842, AY154902, AY154901) were used as outgroup taxa, and for molecular clock calibration. The second data set (nuclear data set) included partial sequences of two nuclear loci (suppressor of SWI4 1 and clone 17) of a subset of 56 individuals representing the major lineages of P. hispanicus (as per previous mt analyses), 16 individuals of P. algirus, and six specimens of G. caesaris. A total of 16 individuals were heterozygotes in 3.8% of their nucleotide positions, which were coded as N in all subsequent analyses. The third data set (combined data set) included cytb, nad4, and the two nuclear loci (suppressor of SWI4 1 and clone 17) for the subset of 54 individuals representing the major lineages of P. hispanicus.
Sequences were aligned using Clustal × version 1.83  with default penalties for gap opening and gap extension, and alignments were visually verified. For each molecular marker, independent alignments were prepared, and the best-fit models of nucleotide substitution were inferred using the Akaike information criterion (AIC;  as implemented in Modeltest version 3.7 ). The mt, nuclear, and combined data sets were analyzed using maximum likelihood (ML; , and Bayesian inference (BI; . ML analyses were performed with RAxML version 7.2.6  using the rapid hill-climbing algorithm  and starting from 100 distinct randomized maximum-parsimony starting trees. For BI analyses, we used MrBayes version 3.1.2 [63, 66]. We ran four simultaneous Markov chains for 20 million generations, sampling every 2000 generations (10,000 trees), and discarding the first 10% of generations (1,000 trees) as burn-in to prevent sampling before reaching stationarity. Adequate convergence of the Bayesian Markov chain Monte Carlo runs was assessed using Tracer version 1.5 (http://tree.bio.ed.ac.uk/software/tracer/). Two independent BI runs were performed to increase the chance of adequate mixing of the Markov chains, and to give some chance of spotting failure to converge. For both ML and BI analyses, three partitions were used for the mt data set (accounting for each codon position), two partitions were used for the nuclear data set (one per locus), and eight were used for the combined data set (accounting for each nuclear marker, plus each codon position of each mt marker). For BI analyses, independent best-fit models of nucleotide substitution (as selected by Modeltest) were used for each partition with model parameters unlinked and estimated separately among partitions. For ML analyses, the GTR + Γ model was used for all partitions due to software (RAxML) constraints, and model parameters were unlinked and estimated separately among partitions. Statistical support for internal branches in the ML analyses was evaluated by non-parametric bootstrapping  with 2,000 replicates and using posterior probabilities in the BI analyses.
Dating of Divergence Times
The combined dataset was used to date major cladogenetic events within the Psammodromus phylogeny using the Bayesian relaxed clock method  as implemented in BEAST version 1.6.1 . This widely used method for dating phylogenies  assumes a relaxed uncorrelated clock with rates drawn from a lognormal distribution across branches. The ML optimal topology was used as a starting tree, and the birth-death process  was used to describe diversification. The partitions used for BI of the combined data set and corresponding models (see above) were employed for the dating analysis. A first run of 10 million generations was first performed to optimize the scale factors of the prior function. The final Markov chain was run twice for 100 million generations, sampling every 10,000 generations, and burn-in and convergence of the chains were determined with Tracer. Effective Sample Size (ESS) values were over 350 for all parameters sampled.
Time estimates were calibrated using the formation of El Hierro 1.12 ± 0.02 Mya  as internal time constraint (maximum age) for the split between Gallotia lizards from El Hierro and La Gomera (Canary Islands). The time constraint was used as 'soft' bound : the mean and standard deviation of the normal distribution were chosen so that the mean is the arithmetical mean of the interval, and 95% of the probability lies within the lower and upper bound. Analyses using a uniform instead of a normal distribution yielded almost the same average time estimates and confidence interval width (not shown).
Population genetic analyses of P. hispanicus were performed using Arlequin 3.11. . Minimum-spanning networks based on cytb, nad4, suppressor of SWI4 1 and clone 17 haplotypes were inferred independently. To infer population dynamics from cytb haplotype data, each of the three major lineages of P. hispanicus (as recovered in the phylogenetic trees and the minimum-spanning networks) was split into northern and southern populations (approximately above and below 40ºN separating the three lineages into approximately 50% northern and 50% southern populations; see Additional File 1). Descriptive statistics including haplotype diversity (Hd; ) and nucleotide diversity (π; ) were determined for northern and southern populations, respectively. Inter-haplotype levels of divergence between northern and southern populations were estimated using the fixation index ΦST , which includes information on mitochondrial haplotype frequency , and genetic distances (TrN93; ) with gamma correction. Significance of pairwise population comparisons was tested by 20,000 permutations. An analysis of molecular variance (AMOVA) was used to examine the amount of genetic variability partitioned within and among populations . AMOVA tests were organized in a hierarchical manner so that population structure was studied at increments of increasing spatial scale, which range from structure within and among different populations to northern versus southern groups. Permutation procedures (N = 20,000) were used to construct null distributions, and test the significance of variance components for each hierarchical comparison . In all instances with multiple tests, p-values were adjusted using the sequential Bonferroni correction .
Immediately after capture, standardized digital photographs of each lizard belonging to P. hispanicus were taken according to the methods used by Fitze & Richner . In brief, living lizards were placed in a box covered with a photographic filter lens to immobilize the individual. This box was placed in an opaque camera box where two flashes were mounted. The settings of the camera and flashes were always identical and the distance between the objective and the object was fixed. Thus all photographs received a standard light exposure. Standard white chips (Kodak Colour Control Patches with R = 255, G = 255, B = 255) were fixed to each side of the filter for detecting potential errors in light exposure and allowing for calibrating the sizes. Photos were taken of the lizards' belly, back and flanks. Lizards were weighed to the nearest 0.001 g, snout-to-vent length (SVL) and total length were measured to the nearest 1 mm, and the number of femoral pores was counted.
List and brief description of the phenotypic measurements taken
snout to vent length
Total length (mm)
snout to tail tip length
Body mass (g)
Head length (mm)
distances between the tip of the snout and the occipital edge
Head width (mm)
distances between the borders of the outermost left and right
supraocular scales (located behind the eyes).
degree of head sharpness. Head width/head length 
Snout width (mm)
distance between the left and right foremost intersection point of
the first supraocular and the first supraciliar scale
Snout length (mm)
distance between the tip of the snout and the orthogonal
intersection with the snout width
degree of snout sharpness. Snout length/snout width
Anal scale width (mm)
distance between the posterior borders of the anal scale
Relative anal scale
anal scale width/SVL 
mean number of right and left femoral pores
number of longitudinal ventral scale rows
number supralabial scales below subocular scale 
number throat scales
number well-differentiated collar scales
Number of ocelli
mean number of left and right ocels
sum of presence/absence (1/0) of green coloration on neck, belly, subocular and supralabial scales + number of green coloured longitudinal lines + number of green coloured longitudinal lines that spread further than the middle of the body/2
Statistics used for the Analyses of Phenotypic Data
All statistical analyses were conducted using R 2.7.0 software (Free Software Foundation, GNU Project, Boston, MA, USA). A total of 211 adults were used for the multivariate analyses. We applied a permutational MANOVA (NP - MANOVA) based on distance measures  to investigate differences between lineages revealed by phylogenetic analyses. A total of 9,999 permutations were conducted, following Manly . We first standardized the data in order to avoid differential impact of unequally scaled variables on the posterior analysis . Thereafter, we calculated the dissimilarities between observations based on Euclidean distances, and applied a NP - MANOVA (Adonis function in Vegan package). Results from paired contrasts between lineages were corrected using Bonferroni procedures, indicated as P adj . The assumption of homogeneity of multivariate dispersion between lineages was fulfilled for all presented analyses .
Discriminant functions were derived from the linear combinations of the variables [88, 89], to assess the relative contribution of each variable to the differences between lineages and to visualize multivariate differences between lineages.
We also run univariate ANOVAs on each response variable separately in order to understand which traits differed among lineages. All assumptions were verified and, if necessary, transformations or non-parametric analyses were applied. Model assumptions were fulfilled in all cases.
Ecological Niche Modelling
To obtain predictive maps for the potential distributions of each P. hispanicus lineage and assess the degree of ecological divergence among lineages, we performed species distribution modelling based on the ecological niche concept . We built models based on environmental predictors and the lineage presences resulting from the phylogenetic analyses (9, 5, and 8 sampled populations for the edwardsianus, Central, and Western lineage respectively; Additional File 1).
We used biologically relevant environmental variables at a 1-km resolution. We initially considered 18 climatic, one topographic and two vegetation index variables, all of them typically being used in biogeographic models as direct and/or indirect predictors of species distributions (e.g. [91–93]). Climatic and topographic variables were obtained from the Worldclim source , whose raster maps were downloaded at the EDIT Geoplatform http://edit.csic.es/. As a measurement of primary productivity we used monthly maps of Enhanced Vegetation Indexes (EVI) generated from satellite MODIS images available at the NASA-LP DAAC web page (https://lpdaac.usgs.gov/lpdaac). EVI is an improvement of Normalized Difference Vegetation Index that minimizes canopy background variations and maintains sensitivity over dense vegetation conditions. Besides predicting primary productivity, EVI is also a measurement of shade availability, which in turn, may modulate direct climate effects in ectotherms (e.g. ). We generated year-averaged monthly values of EVI for each cell from the oldest period available to the year of taxon sampling (2000-2006). Thereafter, we calculated for each cell the minimum and maximum values over all months.
To meet model assumptions, environmental variables were Box-Cox transformed. Redundancy and colinearity between variables was analysed using Pearson correlations. Eight variables were excluded from the analyses because they were highly correlated with other variables (r > 0.90). A total of 13 variables were used for the subsequent analyses: seven temperature predictors (annual mean temperature, mean temperature of wettest quarter, mean temperature of driest quarter, minimum temperature of coldest month, maximum temperature of warmest month, annual temperature range and isothermality), four precipitation predictors (annual precipitation, precipitation of the warmest quarter, precipitation of the coldest quarter and precipitation seasonality), one topographic predictor (elevation) and one EVI index (minimum EVI index).
We built GIS-based models to estimate each lineage's multidimensional niche. We geo-referenced the sampled populations and used digital maps of environmental variables. To obtain habitat suitability (HS) maps for each lineage we used the Ecological Niche Factor Analysis (ENFA) implemented in the GIS-statistical tool Biomapper . ENFA allows the calculation of HS scores for each cell in a gridded map and it is especially suited if absence data are not available, unreliable, or meaningless. ENFA is analogous to principal component analysis with the difference that it is based on the niche concept of marginality (defined as the ecological distance between the lineages optimum and the average habitat within the study area) and specialization (the ratio of the ecological variance in average habitat to that associated to the focal lineage ). Environmental variables are compacted into a few factors where the first factor maximizes the marginality and the other factors maximize the specialization of the focal lineage. Finally, those factors that explained the biggest part of variance (i.e. those that best explained a lineage's ecological range) were used to obtain HS scores for each cell in the map ranging from 0 to 100. The distribution of the eigenvalues was compared with the MacArthur's broken stick distribution to decide which subset of factors was used for HS map computation . In addition to the predictive maps based on continuous HS scores, we obtained binary maps reclassified as predictions of suitability and non-suitability. The threshold applied to transform continuous into binary maps was the minimum training presence (the lowest HS scores associated to the populations of each lineage), which classified as suitable areas those corresponding to scores bigger than zero. This method is conservative since it may include cells where the HS scores of lineages are small.
Evaluation of models
We evaluated the predictive capacity of our niche models using an independent data set, i.e the P. hispanicus species distribution published in the official Spanish Atlas of Amphibians and Reptiles . The binary predictions derived for the three lineages were superimposed to obtain an overall prediction for P. hispanicus. We thereafter calculated the omission error with respect to the known P. hispanicus presences (% cells with predicted absences, with respect to the number of cells with Atlas presences).
A crucial assumption for assessing niche divergences and overlaps among lineages is that the predictive capacity should be similar for the compared lineages. We therefore investigated the robustness of the lineage specific model predictions, using a jack-knife methodology (e.g. ). In brief, we run an individual lineage model with one sampled population excluded (Nlineages - 1) and calculated the proportional prediction difference ((HS N lineages - HS N lineages-1)/HS N lineage). This procedure was repeated for each population. To evaluate potential biases between lineages we used a one-way ANOVA with the proportional prediction difference as the dependent variable and lineage as a factor.
Assessment of environmental predictor relevance
First, we analyzed differences in environmental predictors between lineages of P. hispanicus using one-way ANOVAs. These analyses allow linking phenotypic differences between lineages with differences in potential selective pressures. Second, we estimated for each lineage the relevance of each environmental predictor using ENFA. For each environmental predictor, we determined the highest predictor value among ENFA factors that explained an important part of the variation. Factor importance was determined using the broken stick distribution. Environmental predictor relevance increases with increasing predictor values (Figure 7; for further details see ).
Assessment of ecological niche divergence within P. hispanicus
We applied different principles to assess divergence of niches between lineage pairs, and examined whether there is concordance among procedures. First, we estimated interpredictivity among lineages [99, 100]. We used ecological niche models derived for a given lineage and predict all sampled populations (Figure 9). High predictability for the populations belonging to the lineage for which the niche model was constructed and low predictability for populations of the other lineages (i.e. low interpredictivity) indicates divergent ecological niches. We calculated interpredictivity among lineages using continuous HS scores and binary predictions (see Model building) and used one-way ANOVAs based on angularly transformed scores  for HS scores to compare predictability between lineages. Post-hoc comparisons were performed using Tukey´s range test.
Second, we compared niche models based on split and lumped taxonomic groupings following Raxworthy et al.'s rational . We calculated the excess of prediction (area predicted by lumped model - area predicted by superposed split models) for each lineage pair and investigated whether the combined ecological niche model predicts more suitable area than the superposed split models. We also calculated the false-positive rates of lumped and superposed split models by determining the percentage of predicted niche space, where now presences have been cited  (wrongly predicted presences/total of absences). Since the Atlas does not discriminate among lineages, we estimated the false-positive rates excluding Atlas presence data within the minimum-convex polygon of the not considered lineage. To estimate niche divergence we calculated differences in false-positive rates for each lineage pair (false positive rate of lumped - split models). Bigger false-positive rate differences indicate higher niche divergence.
Finally, the extent of geographic overlap of habitat suitability maps between lineage pairs can reveal information about niche divergence and about the availability and spatial structure of the shared environmental conditions. To estimate overlap, we used binary maps (see model building) and determined areas suitable for two lineages (Figures 8c, d, and 8e). To visualize the spatial structure and the connectivity of shared environmental conditions, we first localized potential lineage contact zones using lineage-specific minimum-convex polygons (Additional File 3 - Figure S3) and considered areas as potential contact zones where two polygons overlapped. Zoomed contact zones of lineage specific predictions were provided (Figures 8a, b, and 8c), that show the extent to which lineages are isolated/connected through continuous HS scores.
Assessing inter-model variability with respect to ENFA
To understand whether the main conclusions derived using ENFA were robust, we assessed inter-model variability using MaxEnt, Multidimensional Envelope (MDE). We examined to which extent niche divergence and geographic overlaps may depend on the used algorithms, and obtained ensemble predictions , which join the results of the three modelling techniques (ENFA, Multidimensional Envelopes and MaxEnt). For more detailed information about the applied techniques and analyses, see Additional File 2.
LMSJ was supported by a PhD fellowship (I3P 060501) from the Consejo Superior de Investigaciones Científicas (CSIC) of Spain co-financed by the European Social Fund, VGJ was supported by a PhD fellowship from the Ministry of Science and Innovation of Spain (FPU AP2006-01678), DSM was supported by a postdoctoral fellowship of the Ministry of Science and Innovation of Spain (MEC/Fulbright 2007-0448), and PSF by a grant from the Ministry of Science and Innovation of Spain (Programa Ramón y Cajal, RYC-2003-006136). The project costs were financed by a grant from the Comunidad de Madrid, Spain (200530M090 to PSF and RZ).
We thank Borja Milá, John Wiens, and two anonymous reviewers for insightful comments on a previous version of the manuscript. We also thank José María Delgado, Benet Pera Gresely, Guillem Pérez i de Lanuza, and Tomás Ponce for field assistance, Mario García-París, Jesús Mellado, Gustavo A. Llorente, Fernando Palacios, and Juan Manuel Pleguezuelos for indicating locations where they observed specimens of P. hispanicus. The Asociación Herpetológica Española for providing us with the data gathered for the Spanish Atlas of Reptiles and Amphibians . Ana Cristina Andreu Rubio, Adolfo Marco Llorente, and Isidro Román Maudo for helping with permits and logistics at the Doñana National park.
- Losos JB: Integrative approaches to evolutionary ecology - Anolis lizards as model systems. Annu Rev Ecol Syst. 1994, 25: 467-493. 10.1146/annurev.es.25.110194.002343.Google Scholar
- Wiens JJ, Graham CH: Niche conservatism: Integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol Syst. 2005, 36: 519-539. 10.1146/annurev.ecolsys.36.102803.095431.Google Scholar
- McCormack JE, Zellmer AJ, Knowles LL: Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation?: Insights from tests with niche models. Evolution. 2010, 64: 1231-1244.PubMedGoogle Scholar
- Pease KM, Freedman AH, Pollinger JP, McCormack JE, Buermann W, Rodzen J, Banks J, Meredith E, Bleich VC, Schaefer RJ, et al: Landscape genetics of California mule deer (Odocoileus hemionus): the roles of ecological and historical factors in generating differentiation. Mol Ecol. 2009, 18 (9): 1848-1862. 10.1111/j.1365-294X.2009.04112.x.PubMedGoogle Scholar
- Mayr E: Systematics and the Origin of Species from the Viewpoint of a Zoologist. 1942, New York: Columbia University PressGoogle Scholar
- Coyne JA, Orr HA: Speciation. 2004, Cambridge, MA: Harvard Univ. PressGoogle Scholar
- Graham CH, Ron SR, Santos JC, Schneider CJ, Moritz C: Integrating phylogenetics and environmental niche models to explore speciation mechanisms in dendrobatid frogs. Evolution. 2004, 58: 1781-1793.PubMedGoogle Scholar
- Chaves JA, Pollinger JP, Smith TB, LeBuhn G: The role of geography and ecology in shaping the phylogeography of the speckled hummingbird (Adelomyia melanogenys) in Ecuador. Mol Phylogenet Evol. 2007, 43: 795-807. 10.1016/j.ympev.2006.11.006.PubMedGoogle Scholar
- Wiens JJ, Donoghue MJ: Historical biogeography, ecology and species richness. Trends in Ecology & Evolution. 2004, 19: 639-644. 10.1016/j.tree.2004.09.011.Google Scholar
- Kozak KH, Wiens JJ: Does niche conservatism promote speciation? A case study in North American salamanders. Evolution. 2006, 60: 2604-2621.PubMedGoogle Scholar
- Losos JB: Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. Ecol Lett. 2008, 11: 995-1003. 10.1111/j.1461-0248.2008.01229.x.PubMedGoogle Scholar
- Thomassen HA, Buermann W, Mila B, Graham CH, Cameron SE, Schneider CJ, Pollinger JP, Saatchi S, Wayne RK, Smith TB: Modeling environmentally associated morphological and genetic variation in a rainforest bird, and its application to conservation prioritization. Evol Appl. 2010, 3: 1-16. 10.1111/j.1752-4571.2009.00093.x.PubMedPubMed CentralGoogle Scholar
- Orr MR, Smith TB: Ecology and speciation. Trends in Ecology & Evolution. 1998, 13: 502-506. 10.1016/S0169-5347(98)01511-0.Google Scholar
- Schluter D: The Ecology of Adaptive Radiation. 2000, Oxford: Oxford University PressGoogle Scholar
- Pleguezuelos JM, Marquez R, Lizana M, (eds): Atlas y Libro Rojo de los Anfibios y Reptiles de España. 2002, Madrid: Ministerio de Medio AmbienteGoogle Scholar
- Verdú-Ricoy J, Carranza S, Salvador A, Busack SD, Díaz JA: Phylogeography of Psammodromus algirus (Lacertida) revisited: systematic implications. Amphib Rept. 2010, 31: 576-582. 10.1163/017353710X521555.Google Scholar
- Busack SD, Lawson R: Historical biogeography, mitochondrial DNA, and allozymes of Psammodromus algirus (Lacertidae): a preliminary hypothesis. Amphib Rept. 2006, 27: 181-193. 10.1163/156853806777239968.Google Scholar
- Busack SD, Salvador A, Lawson R: Two new species in the genus Psammodromus (Reptilia: lacertidae) from the Iberian Peninsula. Ann Carnegie Mus. 2006, 75: 1-10. 10.2992/0097-4463(2006)75[1:TNSITG]2.0.CO;2.Google Scholar
- Raxworthy CJ, Ingram CM, Rabibisoa N, Pearson RG: Applications of ecological niche modeling for species delimitation: A review and empirical evaluation using day geckos (Phelsuma) from Madagascar. Syst Biol. 2007, 56: 907-923. 10.1080/10635150701775111.PubMedGoogle Scholar
- Pérez-Mellado V: Psammodromus hispanicus Fitzinger, 1826. Fauna Ibérica. Edited by: Salvador A. 1998, Madrid: Museo Nacional de Ciencias Naturales-CSIC, 318-326.Google Scholar
- Hewitt GM: Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996, 58: 247-276.Google Scholar
- Hewitt G: The genetic legacy of the Quaternary ice ages. Nature. 2000, 405: 907-913. 10.1038/35016000.PubMedGoogle Scholar
- Carranza S, Harris DJ, Arnold EN, Batista V, de la Vega JPG: Phylogeography of the lacertid lizard, Psammodromus algirus, in Iberia and across the Strait of Gibraltar. J Biogeogr. 2006, 33: 1279-1288. 10.1111/j.1365-2699.2006.01491.x.Google Scholar
- Weijermars R: Geology and tectonics of the betic zone, SE Spain. Earth-Science Reviews. 1991, 31: 153-236. 10.1016/0012-8252(91)90019-C.Google Scholar
- Martinez-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. Mol Ecol. 2006, 15: 3375-3388. 10.1111/j.1365-294X.2006.03013.x.PubMedGoogle Scholar
- Martinez-Solano I, Goncalves HA, Arntzen JW, Garcia-Paris M: Phylogenetic relationships and biogeography of midwife toads (Discoglossidae Alytes). J Biogeogr. 2004, 31: 603-618. 10.1046/j.1365-2699.2003.01033.x.Google Scholar
- Albert EM, Zardoya R, Garcia-Paris M: Phylogeographical and speciation patterns in subterranean worm lizards of the genus Blanus (Amphisbaenia: Blanidae). Mol Ecol. 2007, 16: 1519-1531. 10.1111/j.1365-294X.2007.03248.x.PubMedGoogle Scholar
- Gomez A, Lunt DH: Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. Phylogeography of Southern European Refugia - Evolutionary Perspective on the Origins and Conservation of European Biodiversity. Edited by: Weiss SFN. 2007, 155-188.Google Scholar
- Zardoya R, Doadrio I: Molecular evidence on the evolutionary and biogeographical patterns of European cyprinids. J Mol Evol. 1999, 49: 227-237. 10.1007/PL00006545.PubMedGoogle Scholar
- Cunha C, Coelho MM, Carmona JA, Doadrio I: Phylogeographical insights into the origins of the Squalius alburnoides complex via multiple hybridization events. Mol Ecol. 2004, 13: 2807-2817. 10.1111/j.1365-294X.2004.02283.x.PubMedGoogle Scholar
- Doadrio I, Carmona JA: Phylogenetic relationships and biogeography of the genus Chondrostoma inferred from mitochondrial DNA sequences. Mol Phylogenet Evol. 2004, 33: 802-815. 10.1016/j.ympev.2004.07.008.PubMedGoogle Scholar
- Garcia-Paris M, Alcobendas M, Buckley D, Wake DB: Dispersal of viviparity across contact zones in Iberian populations of fire salamanders (Salamandra) inferred from discordance of genetic and morphological traits. Evolution. 2003, 57: 129-143.PubMedGoogle Scholar
- Veith M, Kosuch J, Vences M: Climatic oscillations triggered post-Messinian speciation of Western Palearctic brown frogs (Amphibia, Ranidae). Mol Phylogenet Evol. 2003, 26: 310-327. 10.1016/S1055-7903(02)00324-X.PubMedGoogle Scholar
- Veith M, Mayer C, Samraoui B, Barroso DD, Bogaerts S: From Europe to Africa and vice versa: evidence for multiple intercontinental dispersal in ribbed salamanders (Genus Pleurodeles). J Biogeogr. 2004, 31: 159-171. 10.1111/j.1365-2699.2004.00957.x.Google Scholar
- Pérez Albert A, Valcárcel Díaz M, Blanco Chao R: Pleistocene glaciation in Spain. Quartenary Glaciations - Extent and Chronology Volume 2. Edited by: Ehlers J, Gibbard PL. 2004, London: Elsevier Ltd, 390-394.Google Scholar
- Kukla G: Saalian supercycle, Mindel/Riss interglacial and Milankovitch's dating. Quaternary Science Reviews. 2005, 24: 1573-1583. 10.1016/j.quascirev.2004.08.023.Google Scholar
- 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. Proc R Soc Lond Ser B-Biol Sci. 2001, 268: 1625-1630. 10.1098/rspb.2001.1706.Google Scholar
- Surget-Groba Y, Heulin B, Guillaume GP, Thorpe RS, Kupriyanova L, Vogrin N, Maslak R, Mazzotti S, Venczel M, Ghira I, et al: Intraspecific phylogeography of Lacerta vivipara and the evolution of viviparity. Mol Phylogenet Evol. 2001, 18: 449-459. 10.1006/mpev.2000.0896.PubMedGoogle Scholar
- Hewitt GM: Post-glacial re-colonization of European biota. Natural-Environment-Research-Council-Royal-Society-of-Edinburgh Symposium on Molecular Genetics in Animal Ecology: Sep 1998; Edinburgh, Scotland. 1998, 87-112.Google Scholar
- Quinteiro J, Rodriguez-Castro J, Castillejo J, Iglesias-Pineiro J, Rey-Mendez M: Phylogeny of slug species of the genus Arion: evidence of of Iberian endemics and of the existence of relict species in Pyrenean refuges. J Zool Syst Evol Res. 2005, 43: 139-148. 10.1111/j.1439-0469.2005.00307.x.Google Scholar
- Garcia-Paris M, Jockusch EL: A mitochondrial DNA perspective on the evolution of Iberian Discoglossus (Amphibia: Anura). J Zool. 1999, 248: 209-218.Google Scholar
- de Bruijne CH, Andriessen PAM: Far field effects of Alpine plate tectonism in the Iberian microplate recorded by fault-related denudation in the Spanish Central System. 9th International Conference on Fission Track Dating and Thermochronology: Feb 06-11 2002; Lome, Australia. 2002, 161-184.Google Scholar
- Wiens JJ: Speciation and ecology revisited: phylogenetic niche conservatism and the origin of species. Evolution. 2004, 58: 193-197.PubMedGoogle Scholar
- Endler JA: A predator's view of animal color patterns. Evol Biol. 1978, 11: 319-364.Google Scholar
- Cuthill IC, Stevens M, Sheppard J, Maddocks T, Parraga CA, Troscianko TS: Disruptive coloration and background pattern matching. Nature. 2005, 434: 72-74. 10.1038/nature03312.PubMedGoogle Scholar
- Soule M, Kerfoot W: On the climatic determination of scale size in a lizard. Syst Zool. 1971, 21: 97-105.Google Scholar
- Swenson NG, Enquist BJ, Thompson J, Zimmerman JK: The influence of spatial and size scale on phylogenetic relatedness in tropical forest communities. Ecology. 2007, 88: 1770-1780. 10.1890/06-1499.1.PubMedGoogle Scholar
- Desdevises Y, Legendre P, Azouzi L, Morand S: Quantifying phylogenetically structured environmental variation. Evolution. 2003, 57: 2647-2652.PubMedGoogle Scholar
- Albertson RC, Markert JA, Danley PD, Kocher TD: Phylogeny of a rapidly evolving clade: The cichlid fishes of Lake Malawi, East Africa. Proc Natl Acad Sci USA. 1999, 96: 5107-5110. 10.1073/pnas.96.9.5107.PubMedPubMed CentralGoogle Scholar
- Rüber L, Van Tassel JL, Zardoya R: Rapid speciation and ecological divergence in the American seven-spinned gobies (Gobiidae, Gobiosomatini) inferred from a molecular phylogeny. Evolution. 2003, 57: 1584-1598.PubMedGoogle Scholar
- Ricklefs RE, Birmingham E: The causes of evolutionary radiations in archipelagoes: Passerine birds in the Lesser Antilles. Am Nat. 2007, 169: 285-297. 10.1086/510730.PubMedGoogle Scholar
- Thorpe RS, Surget-Groba Y, Johansson H: Genetic Tests for Ecological and Allopatric Speciation in Anoles on an Island Archipelago. PLoS Genet. 2010, 6: e1000929-10.1371/journal.pgen.1000929.PubMedPubMed CentralGoogle Scholar
- Thorpe RS, Surget-Groba Y, Johansson H: The relative importance of ecology and geographic isolation for speciation in anoles. Philos Trans R Soc B-Biol Sci. 2008, 363: 3071-3081. 10.1098/rstb.2008.0077.Google Scholar
- 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. Proc Natl Acad Sci USA. 1989, 86: 6196-6200. 10.1073/pnas.86.16.6196.PubMedPubMed CentralGoogle Scholar
- San Mauro D, Gower DJ, Oommen OV, Wilkinson M, Zardoya R: Phylogeny of caecilian amphibians (Gymnophiona) based on complete mitochondrial genomes and nuclear RAG1. Mol Phylogenet Evol. 2004, 33: 413-427. 10.1016/j.ympev.2004.05.014.PubMedGoogle Scholar
- Arevalo E, Davis SK, Sites JW: Mitochondrial-DNA Sequence Divergence and Phylogenetic-Relationships among 8 Chromosome Races of the Sceloporus-Grammicus Complex (Phrynosomatidae) in Central Mexico. Syst Biol. 1994, 43: 387-418.Google Scholar
- Murphy R, Sites J, Buth D, Haufler C: Molecular Techniques. Molecular Systematics. Edited by: Hillis DM MC, Mable BK. 1996, Sunderland: Sinauer Associates, 51-381.Google Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990, 215: 403-410.PubMedGoogle Scholar
- Thompson JD, Gibson TJ, Plewniak F, Jeanmougin J, Higgins DG: The CLUSTAL × windows interface: Flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Research. 1997, 25: 4876-4882. 10.1093/nar/25.24.4876.PubMedPubMed CentralGoogle Scholar
- Akaike H: Information theory as an extension of the maximum likelihood principle. Second international symposium of information theory. Edited by: Petrov BN, Csaki F. 1973, Budapest, Hungary: Akademiai KiadoGoogle Scholar
- Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.PubMedGoogle Scholar
- Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Evol. 1981, 17: 368-376. 10.1007/BF01734359.PubMedGoogle Scholar
- Huelsenbeck JP, Ronquist FR: MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics. 2001, 17: 754-755. 10.1093/bioinformatics/17.8.754.PubMedGoogle Scholar
- Stamatakis A: RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006, 22: 2688-2690. 10.1093/bioinformatics/btl446.PubMedGoogle Scholar
- Stamatakis A, Blagojevic F, Nikolopoulos D, Antonopoulos C: Exploring new search algorithms and hardware for phylogenetics: RAxML meets the IBM Cell. The Journal of VLSI Signal Processing. 2007, 48: 271-286. 10.1007/s11265-007-0067-4.Google Scholar
- Ronquist F, Huelsenbeck JP: MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.PubMedGoogle Scholar
- Felsenstein J: Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985, 39: 783-791. 10.2307/2408678.Google Scholar
- Drummond AJ, Ho SYW, Phillips MJ, Rambaut A: Relaxed phylogenetics and dating with confidence. PLoS Biology. 2006, 4: 699-710.Google Scholar
- Drummond AJ, Rambaut A: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.PubMedPubMed CentralGoogle Scholar
- San Mauro D, Agorreta A: Molecular systematics: a synthesis of the common methods and the state of knowledge. Cell Mol Biol Lett. 2010, 15: 311-341. 10.2478/s11658-010-0010-8.PubMedGoogle Scholar
- Gernhard T: The conditioned reconstructed process. J Theor Biol. 2008, 253: 769-778. 10.1016/j.jtbi.2008.04.005.PubMedGoogle Scholar
- Guillou H, Carracedo JC, Pérez Torrado F, Rodriguez Badiola E: K-Ar ages and magnetic stratigraphy of a hotspot-induced, fast grown oceanic island: El Hierro, Canary Islands. J Volcanol Geoth Res. 1996, 73: 141-155. 10.1016/0377-0273(96)00021-2.Google Scholar
- Yang Z, Rannala B: Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds. Mol Biol Evol. 2006, 23: 212-226.PubMedGoogle Scholar
- Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol Bioinf online. 2005, 1: 47-50.Google Scholar
- Nei M: Molecular evolutionary genetics. 1987, New York: Columbia University PressGoogle Scholar
- Excoffier L, Smouse PE, Quattro JM: Analysis of molecular variance inferred from metric distances among DNA haplotypes - application to human mitochondrial-DNA restriction data. Genetics. 1992, 131: 479-491.PubMedPubMed CentralGoogle Scholar
- Weir BS, Cockerham CC: Estimating F-Statistics for the analysis of population-structure. Evolution. 1984, 38: 1358-1370. 10.2307/2408641.Google Scholar
- Tamura K, Nei M: Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993, 10: 512-526.PubMedGoogle Scholar
- Guo SW, Thompson EA: Performing the exact test of Hardy-Weinberg proportion for multiple alleles. Biometrics. 1992, 48: 361-372. 10.2307/2532296.PubMedGoogle Scholar
- Rice WR: Analyzing tables of statistical tests. Evolution. 1989, 43: 223-225. 10.2307/2409177.Google Scholar
- Fitze PS, Richner H: Differential effects of a parasite on ornamental structures based on melanins and carotenoids. Behav Ecol. 2002, 13: 401-407. 10.1093/beheco/13.3.401.Google Scholar
- Harmon LJ, Gibson R: Multivariate phenotypic evolution among island and mainland populations of the ornate day gecko, Phelsuma ornata. Evolution. 2006, 60: 2622-2632.PubMedGoogle Scholar
- Pérez-Mellado V, Gosá A: Biometría y folidosis en Lacertidae (Sauria, Reptilia), algunos aspectos metodológicos. Rev Esp Herp. 1988, 3: 103-119.Google Scholar
- Blasco M: Contribución al conocimiento de los lacertidos de Andalucía. 1974, Granada: Universidad de GranadaGoogle Scholar
- Lessells CM, Boag PT: Unrepeatable repeatabilites: a common mistake. Auk. 1987, 104: 116-121.Google Scholar
- Anderson MJ: A new method for non-parametric multivariate analysis of variance. Austral Ecology. 2001, 26: 32-46.Google Scholar
- Manly BFJ: Randomization, Bootstrap and Monte Carlo Methods in Biology. 2001, Boca Raton: Chapman & Hall/CRCGoogle Scholar
- Quinn GP, Keough MJ: Experimental design and data analysis for biologists. 2002, Cambridge: Cambridge University PressGoogle Scholar
- Maindonald J, Braun WJ: Data Analysis and Graphics Using R - An example Approach. 2003, Cambridge: Cambridge University PressGoogle Scholar
- Hutchinson GE: Concluding remarks, Cold Spring Harbor Symposium. Quant Biol. 1957, 22: 415-427.Google Scholar
- Rissler LJ, Apodaca JJ: Adding more ecology into species delimitation: Ecological niche models and phylogeography help define cryptic species in the black salamander (Aneides flavipunctatus). Syst Biol. 2007, 56: 924-942. 10.1080/10635150701703063.PubMedGoogle Scholar
- Aragón P, Lobo JM, Olalla-Tárraga MÁ, Rodríguez MÁ: The contribution of contemporary climate to ectothermic and endothermic vertebrate distributions in a glacial refuge. Global Ecology and Biogeography. 2009, 19: 40-49.Google Scholar
- Guisan A, Thuiller W: Predicting species distribution: offering more than simple habitat models. Ecol Lett. 2005, 8: 993-1009. 10.1111/j.1461-0248.2005.00792.x.Google Scholar
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25: 1965-1978. 10.1002/joc.1276.Google Scholar
- Kearney M, Richard S, P PW: The potential for behavioral thermoregulation to buffer ''cold-blooded'' animals against climate warming. Proc Natl Acad Sci USA. 2009, 106: 3835-3840. 10.1073/pnas.0808913106.PubMedPubMed CentralGoogle Scholar
- Hirzel AH, Hausser J, Chessel D, Perrin N: Ecological-niche factor analysis: How to compute habitat-suitability maps without absence data?. Ecology. 2002, 83: 2027-2036. 10.1890/0012-9658(2002)083[2027:ENFAHT]2.0.CO;2.Google Scholar
- Vaughan IP, Ormerod SJ: The continuing challenges of testing species distribution models. J Appl Ecol. 2005, 42: 720-730. 10.1111/j.1365-2664.2005.01052.x.Google Scholar
- Aragón P, Baselga A, Lobo JM: Global estimation of invasion risk zones for the western corn rootworm Diabrotica virgifera virgifera: integrating distribution models and physiological thresholds to assess climatic favourability. J Appl Ecol. 2010, 47: 1026-1035. 10.1111/j.1365-2664.2010.01847.x.Google Scholar
- Peterson AT, Soberón J, Sánchez-Cordero V: Conservatism of Ecological Niches in Evolutionary Time. Science. 1999, 285: 1265-1267. 10.1126/science.285.5431.1265.PubMedGoogle Scholar
- Peterson AT, Holt RD: Niche differentiation in Mexican birds: using point occurences to detect ecological innovation. Ecol Lett. 2003, 6: 774-782. 10.1046/j.1461-0248.2003.00502.x.Google Scholar
- Araújo MB, M N: Ensemble forecasting of species distributions. Trends in Ecology & Evolution. 2007, 22: 42-47. 10.1016/j.tree.2006.09.010.Google Scholar
- Hartigan JA, Wong MA: Algorithm AS 136: A K-Means Clustering Algorithm. App Stat. 1979, 28: 100-108. 10.2307/2346830.Google Scholar
- Salvador A: Reptiles. Fauna Ibérica. Edited by: Ramos MAeaE. 1998, Madrid: Museo Nacional de Ciencias Naturales, CSIC, 10: 705-Google Scholar
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 (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.