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 .