The combined use of raw and phylogenetically independent methods of outlier detection uncovers genome‐wide dynamics of local adaptation in a lizard

Abstract Local adaptation is a dynamic process by which different allele combinations are selected in different populations at different times, and whose genetic signature can be inferred by genome‐wide outlier analyses. We combined gene flow estimates with two methods of outlier detection, one of them independent of population coancestry (CIOA) and the other one not (ROA), to identify genetic variants favored when ecology promotes phenotypic convergence. We analyzed genotyping‐by‐sequencing data from five populations of a lizard distributed over an environmentally heterogeneous range that has been changing since the split of eastern and western lineages ca. 3 mya. Overall, western lizards inhabit forest habitat and are unstriped, whereas eastern ones inhabit shrublands and are striped. However, one population (Lerma) has unstriped phenotype despite its eastern ancestry. The analysis of 73,291 SNPs confirmed the east–west division and identified nonoverlapping sets of outliers (12 identified by ROA and 9 by CIOA). ROA revealed ancestral adaptive variation in the uncovered outliers that were subject to divergent selection and differently fixed for eastern and western populations at the extremes of the environmental gradient. Interestingly, such variation was maintained in Lerma, where we found high levels of heterozygosity for ROA outliers, whereas CIOA uncovered innovative variants that were selected only there. Overall, it seems that both the maintenance of ancestral variation and asymmetric migration have counterbalanced adaptive lineage splitting in our model species. This scenario, which is likely promoted by a changing and heterogeneous environment, could hamper ecological speciation of locally adapted populations despite strong genetic structure between lineages.


| INTRODUC TI ON
Local adaptation, the process by which natural selection leads a population of organisms to be better suited to its particular environment than other populations of the same species, requires that different gene combinations are selected in different populations at different times during their evolutionary history (Endler, 1977;Hereford, 2009;Savolainen, Lascoux, & Merilä, 2013;Schluter, 2001;Wright, 1932). It is thus relevant to investigate the genetic footprint of past and present events that may shape adaptation as a dynamical process (Kokko et al., 2017). Such dynamism relies not only on how local adaptation interacts with current genetic dynamics, but also on how local adaptation has been favored, constrained, and/or blurred by the evolutionary history of populations (Laurent et al., 2016;McEntee, Tobias, Sheard, & Burleigh, 2018;Rosenblum, Hickerson, & Moritz, 2007;Rundle & Nosil, 2005). This temporal perspective can be provided by several genomic methods that in that respect supersede other approaches such as common garden or transplant experiments (Haasl & Payseur, 2016;Tiffin & Ross-Ibarra, 2014).
It is possible to screen genome-wide patterns of DNA polymorphism to detect locus-specific signals of divergent selection by using theoretical predictions about the effects of positive natural selection on allele frequencies (Kim & Stephan, 2002;Luikart, England, Tallmon, Jordan, & Taberlet, 2003;Przeworski, 2002;Schlötterer, 2003). Such screening may reveal loci reflecting high levels of genetic differentiation between groups of individuals putatively subject to divergent selection (Han et al., 2017;Wu, 2001). However, the genetic underpinnings of adaptation may depend on the particular genetic background of populations, which have the potential to develop consistent phylogeographic patterns depending on how long they have been reproductively isolated from one another.
Thus, the way in which populations locally adapt will depend on their phylogeographic history as realized in terms of past genetic change (either adaptive or stochastic, or influenced by gene flow; Barrett & Schluter, 2008;Han et al., 2017). This may often complicate the detection and interpretation of loci subject to selection (Rhode, Bester-Van Der Merwe, & Roodt-Wilding, 2017;Schield et al., 2017;Tigano, Shultz, Edwards, Robertson, & Friesen, 2017).
The most immediate approach to outlier analysis consists of uncovering which loci present traces of selection, and whether this selection is divergent or not (Tiffin & Ross-Ibarra, 2014). To this end, Bayesian approaches that use logistic regression to partition F ST coefficients into a population-specific term and a locus-specific term (e.g., Foll & Gaggiotti, 2008) have proved convenient. It is also possible to find outliers whose allele frequencies deviate from the expected under genetic drift alone, given the population tree (e.g., Bonhomme et al., 2010). For simplicity, we refer to these two approaches as raw outlier analysis (ROA) and coancestry-independent outlier analysis (CIOA), respectively. If the two sets of outliers detected by these alternative approaches do not overlap, ROA and CIOA, combined with gene flow estimates, can shed light on the relative roles of ancestral genetic variation, admixture, and genetic innovations, as drivers of local adaptation dynamics (Barrett & Schluter, 2008;Haasl & Payseur, 2016). The lack of overlap between the outliers detected by ROA and CIOA could be meaningful in ecology-driven phylogenetic structuring scenarios. This is because the ability of CIOA to detect SNPs related to the adaptive processes should be low when the relationship between ecology and genetic structure in the sample is strong (Bonhomme et al., 2010). In the limit, coancestry could be able to explain entirely the pattern of allele fixation. In these cases, such outliers should be detected more easily by ROA, while the ones underlying particular scenarios of local adaptation within the general branching pattern (independent to coancestry) would be more easily detected by CIOA (Haasl & Payseur, 2016).
We used genotyping by sequencing (GBS; Elshire et al., 2011) and combined ROA and CIOA to explore how the fixation of ancestral genetic variation and genetic innovation underlie convergent local adaptation. To do so, we used as model species the most abundant and widespread lizard in the Mediterranean region of the Iberian Peninsula: the large Psammodromus Psammodromus algirus  Figure 1). This species shows consistent genetic and phenotypic divergence associated with environmental heterogeneity across its distribution range. The northwest part of its range is more humid and temperate than the southeast, and this is mirrored by broad changes in the vegetation patterns, with forests dominating in the west and more open spaces in the east . The genetic diversity of P. algirus is structured in two mtDNA lineages, eastern and western, separated ca. 3-3.5 Ma (Carranza, Harris, Arnold, Batista, & Gonzalez De La Vega, 2006). Each mtDNA lineage typically shows a distinct phenotype for a heritable trait that could be adaptively linked to crypsis: lizards that inhabit the eastern, more open and grassier regions tend to display a striped dorsal pattern absent among western lizards, which inhabit more vegetated forests . However, striped and unstriped phenotypes do not perfectly match the eastern-western lineage dichotomy. For instance, the population in Lerma that has been the focus of several previous studies (see Díaz, Pérez-Tris, Tellería, Carbonell, & Santos, 2005;Santos, Díaz, Pérez-tris, Carbonell, & Tellería, 2008; we would expect if they were locally adapted to the forests they live in. However, the analysis of their mtDNA places Lerma within the E2 clade of the eastern lineage, which is almost exclusively composed of striped lizards . Therefore, the Lerma population provides an excellent opportunity to investigate the sources of variation that can be locally selected to generate phenotypic convergence after long periods of genetic divergence. First, ancestral genetic variation standing within lineages since before their divergence could allow populations to face high levels of environmental heterogeneity (i.e., adaptation by ancestral genetic diversity; Barrett & Schluter, 2008). This source of ancestral genetic variation should be detected by ROA in the form of divergent outliers with high levels of homozygosity in populations at the extremes of the phenotypic (and environmental) gradient, but with higher than average heterozygosity in the Lerma population. If this was the case, adaptive genetic variation would predate lineage splitting, when it could be either neutral, adaptive to an intermediate environment or maintained by adaptation to a geographically unstructured environmental heterogeneity (Huang, Tran, & Agrawal, 2016;Mackay, 1981;McDonald & Ayala, 1974;Paccard et al., 2018;Parter, Kashtan, & Alon, 2008;Welch & Jiggins, 2014). This variation would be divergently selected in the different ecological scenarios occupied by where eastern and western lineages . Therefore, we would expect an overrepresentation of alternative alleles in lizards inhabiting each extreme of the environmental gradient (Nosil, Funk, & Ortiz-Barrientos, 2009). Such overrepresentation should be consistent with the branching pattern of populations, since lineage splitting seemed to be ecologically driven . However, these outliers should be less detectable by CIOA, because they would not show allele frequencies with enough deviation from the expected given the overall ecologically driven phylogenetic structure of populations. Second, local adaptation of populations faced with novel environments could have been attained through genetic innovation. Genetic innovation might be provided in Lerma by immigrants bringing genetic variants from populations already adapted to similar environments (i.e., adaptation by immigration; Peniston, Barfield, & Holt, 2019). Both eastern and western lineages would show highly divergent homozygosity not only in neutral loci but also in the ones that are under selection, because of both the maintenance of ancestral lineage splitting (genetic drift) and the adaptation to their respective extremes of an environmental gradient, respectively Räsänen & Hendry, 2008). Thus, if genetic innovation is the result of admixture, we also expect ROA to detect outliers with high heterozygosity, but not higher than the average heterozygosity of the entire genome.
Finally, genetic innovation may also evolve from natural selection of locally arisen variants (Huber et al., 2015). In Lerma, such innovations should be uncovered by CIOA in the form of outliers with higher levels of homozygosity than the average for the genome (ROA would have more difficulty to detect such outliers if the phylogenetic inertia of adaptation is large). These scenarios are not mutually exclusive, and their contribution to lizards' local adaptation can be inferred with the outlier analyses described above.
Our goal is not to identify the genetic basis of existing adaptations, but to unravel how evolution supplies the genetic variation that ecology demands for the convergence of phylogeographically divergent populations that are locally adapted to similar habitat. 1,230 m asl) is a more productive montane location covered by deciduous forest where 17% of lizards are striped. Vegetation structure and species composition are very similar at Lerma and the western localities, especially Navacerrada; the three sites share a high cover of Quercus forests (deciduous Q. pyrenaica at Navacerrada and 66% of the Lerma forest fragments, and Q. ilex at El Pardo and 33% of the Lerma fragments) and a well-developed layer of Cistus shrubs (C. laurifolius at Lerma and Navacerrada, and C. ladanifer at El Pardo).

| Phenotypic variation
During the breeding season of 2012, we took tissue samples from 18 adult lizards (>1 year old) noosed or captured by hand in Lerma. During the breeding season of 2015, we captured and sampled with the same methods another 77 lizards (20, 18, 19, and 20 from Aranjuez, Brihuega, El Pardo, and Navacerrada, respectively). These 95 lizards were used for genomic analyses. In 2016, we captured 18 lizards in Lerma that were used for the analysis of phenotypes along with the 77 lizards captured in the other populations during 2015. Therefore, different individuals from Lerma were used for genetic and phenotypic analyses. In order to characterize lizards' phenotype, we brought them to the laboratory to record their snout-vent length (SVL) and body mass and to measure their dorsal coloration from images (these are phenotypic measurements that have been suggested to be adaptive in previously published studies; Díaz et al., 2005Díaz et al., , 2017Iraeta, Salvador, & Díaz, 2013). For this later purpose, we held the individuals stretched under an antireflective glass and immobilized them against a soft sponge to avoid damaging them. We took pictures of the lizards' back in a dark room, using a fix setup with two white light sources placed on either side of the lizard at 25 cm from the subject, with an angle of 45°. The camera was set at 35 cm from the subject to obtain an overhead image, with automatic exposure compensation and shutter speed of 0.5 ms. We used Adobe Photoshop CS6 We standardized the analyzed area using the "magnetic lasso" tool to delimit a 5-cm-long surface from the shoulder joint (which set the width of the analyzed surface) toward the posterior end of the animal. We measured the size of the striped surface with the "magic wand" tool (at 10% tolerance) after clicking at the middle of the stripe. We subsequently used the "similar" option of the "magic wand" tool (at 10% tolerance) to select pigmented scales and used the percentage of colored pixels in the analyzed surface as a measure of the size of the stripe. We measured the darkness of the stripe using the inverse of brightness of the colored layer.
Note that we did not use a color standard to reduce the error in the measurement of brightness, but the fact that all lizards were photographed under the same conditions of illumination and exposure made our images directly comparable.

| DNA extraction, sequencing, and variant calling
We obtained tissue samples by clipping 2 cm of the tail tip of lizards, which were released unharmed at their site of capture. We kept the samples in absolute ethanol at 4°C until DNA extraction. We purified DNA for library preparation using the Speedtools Tissue DNA Extraction Kit (Biotools) with a cell lysis step of 24 hr and resuspension in DNase-free water at 60°C.
We used the restriction enzyme Pst1 for GBS library preparation. Sequencing was done in an Illumina HiSeq2500 sequencer.
The pipeline used for SNPs discovery was UNEAK, implemented in TASSEL v.3.0 (Bradbury et al., 2007), and specifically designed for samples with no reference genome. Sequence tags were aligned to each other to form "networks" of tags, where each node is a single tag sequence, and each edge represents a single base pair difference between two tags. The networks were pruned to remove putative sequencing errors (low frequency alleles) using the error rate threshold parameter. The resulting dataset had 83,648 SNPs with a site depth of 5.96 ± 6.56 (mean ± SD) and a site missingness of 0.49 ± 0.33. We discarded loci with minor allele frequencies <0.01 or that could be successfully sequenced in less than 10% of individuals. The resulting dataset had 73,291 loci, with only one SNP per locus, and with a site depth of 6.60 ± 6.75 and a site missingness of 0.42 ± 0.31.

| Admixture and gene flow
Prior to performing clustering analyses, we used the software PLINK v.1.9. (Purcell et al., 2007) to prune the SNP database for linkage disequilibrium (LD), according to observed sample correlation coefficients. This was necessary because our clustering model (described below) did not take into account LD, and therefore, linked SNPs could bias the grouping. We performed clustering analyses to find genetic structure among populations and lineages (eastern or western) and to obtain the assignment probabilities of each individual to each one of the resulting clusters. We used the program ADMIXTURE v.1.3 (Alexander, Novembre, & Lange, 2009) for maximum-likelihood estimation of individual ancestries. We used cross-validation errors to identify the number of clusters (K-value) for which the model had highest predictive accuracy. However, we also explored higher K values to detect patterns of genetic structure among the five populations. To complement such clustering analysis, we performed a multidimensional scaling analysis (MDS) on TASSEL 5.2. (Bradbury et al., 2007).
In order to improve the interpretation of our outlier analyses, we estimated gene flow among populations to test whether there could be less opposition to immigration of genetically western lizards in Lerma than in other populations of the eastern lineage.
To obtain estimates of gene flow among populations, we used a maximum-likelihood method based on coalescence implemented in MIGRATE 3.6 (Beerli & Felsenstein, 2001). MIGRATE uses an equilibrium model that estimates migration rates averaged across the coalescent history. Including undifferentiated populations in the analysis of migration rates would hamper the estimation of migration rates among less connected populations (Pfenninger & Posada, 2002). Therefore, we pooled individuals from populations that remained undifferentiated in an admixture analysis with K = 5 (equal to the number of sampled populations). Although grouping will make the estimates of local population size not interpretable for pooled populations, this strategy does provide robust estimates of gene flow (Beerli & Felsenstein, 2001). We used MIGRATE to estimate the effective population size scaled by mutation rate Θ = N e μ, together with the effective number of migrants N e m, where N e is the effective population size, μ is the mutation rate per generation, and m is the migration rate per generation. Migration rates were estimated as M = m/μ and presented with 95% credibility intervals. Starting values for all parameter estimates were initially obtained using F ST (Beerli & Felsenstein, 1999); all other searching parameters were set to default values.
We performed a Markov chain Monte Carlo (MCMC) search running 10 heated chains with 100,000 recorded genealogies each, and three cold chains with 50,000 recorded genealogies, sampling each chain every 100 steps with a burn-in period of 10,000 steps.

| Outlier analyses
After estimating overall genetic structure and migration, we applied another filtering step to our SNP dataset to minimize false positives in the subsequent outlier analyses. To this end, we discarded loci with minor allele frequencies <0.05 in each population (thereby excluding all private alleles from the dataset), or loci that could not be successfully sequenced from at least 75% of individuals in each population. The resulting SNP dataset included 6,421 loci.
We used a Bayesian approach to perform a ROA as implemented in Bayescan v.2.1, a very useful approach when the number of populations is low (Foll & Gaggiotti, 2008). Bayescan uses a logistic regression model to partition F ST coefficients into a population-specific term (β) and a locus-specific term (α). We selected loci with α > 0, suggesting positive selection and a conservative false discovery rate (corrected for multiple testing) of q ≤ 0.05. To obtain these parameters, we ran the reversible-jump MCMC algorithm implemented in the program with a prior odd value of 10, and using 20 pilot runs of 5,000 iterations each, followed by 100,000 iterations with a burn-in of 50,000. In order to search for outliers taking into account the variation in the degree of relatedness among populations, we per-  proportion of striped lizards in Lerma (44%) was higher than in western populations, but their stripes were so tenuous that their overall phenotype was clearly western-like (Figure 2a,b).

| Admixture and gene flow
Using 73,291 SNPs, the best model produced by ADMIXTURE had K = 2, recovering two clusters corresponding to the two mtDNA lineages (cv-error = 0.586; Figure 3a). When we forced the model toward higher K values (with higher values of cv-error), Aranjuez was the first population to separate from the rest of its lineage (K = 3, cv-error = 0.665), followed by Brihuega and Lerma (K = 4, cv-error = 0.744). However, when we set K equal to the number of populations in our sample, we did not detect any genetic structure between the two western populations (K = 5, cv-error = 0.810). The MDS using the same dataset produced an axis with eigenvalue = 0.463, which defined an east-west gradient that discriminated between both mitochondrial lineages with p = 1.0 (F 1,93 = 1,369.19, R 2 = .936, p < 10 -6 ; mean ± SD scores of eastern and western lizards = −0.056 ± 0.017 and 0.081 ± 0.019, respectively; Figure 4).
Our maximum-likelihood estimates of gene flow revealed highly asymmetrical patterns of migration among lizard populations (ANOVA for among-group variation in the estimates of immigration:

| Outlier analyses
Our ROA with Bayescan detected 12 outlier loci with α > 0 (0.97 < α < 1.35) and q ≤ 0.05. The CIOA identified nine additional outlier loci with p < .001, none of which was detected by Bayescan.  (Figure 3b). The population of Lerma was composed of individuals assigned either to one of the eastern clusters (the one including Brihuega, with a mean ± SD probability of assignment of 0.393 ± 0.147) or to one of the two western groups (probability of assignment of 0.547 ± 0.155). When we forced the model toward higher cv-errors, the next result was for K = 2 (cv-error = 0.686; Figure 3c), with all lizards being correctly assigned to their mitochondrial lineage except for those from Lerma that were assigned to the western cluster with mean ± SD probability of 0.613 ± 0.140.
Aranjuez and the two western populations, located at opposite extremes of the genetic gradient identified by MDS and admixture analyses, scored high homozygosity (with different alleles approaching fixation) for the outliers detected by ROA. Lerma, and to a lesser extent also Brihuega, scored significantly higher heterozygosity of ROA outliers than expected from overall genome heterozygosity (ANCOVA with population as the factor, H ROA as the dependent variable and H ALL as the covariate: F 4,89 = 10.121, p < .001; Table 1). However, for the outliers detected by CIOA, Lerma scored significantly higher homozygosity than the other populations (H CIOA as the dependent variable: F 4,89 = 5.812, p < .001; Table 1). To further support the hypothesis that genetic variants subject to selection have been fixed along the evolutionary history of highly divergent populations (Aranjuez and western), we examined the probabilities of correct assignment of each individual to its lineage using SNP data with either whole genome data (PG, Figure 3a) or outliers (PO, Figure 3c). In the populations located at the extremes of the genetic gradient, we found a strong correlation between PO and PG (Aranjuez: r = .823, p < .001; western populations: r = .924, p < .001). However, the correlation was not significant in Brihuega (r = .294, p = .236) or in Lerma (r = .154, p = .541).

| D ISCUSS I ON
The large Psammodromus provides an example of local adaptation in which the phenotype of individuals fits better to the expected from habitat characteristics than from population ancestry. Our combined analysis of phenotypic variation, population genetic structure and genetic polymorphisms subject to selection, allowed us to explore the sources of genetic variants favored in different ecological contexts in this system.
The analysis of >73,000 SNPs recovered the east-west division of genetic lineages previously described on the basis of mitochondrial DNA sequences, thereby confirming the ancestral division of this species on a genome-wide basis (Carranza et al., 2006).
Compared to previous studies, our approach boosted the resolution with which the genetic background of lizard populations (and individuals) could be characterized. We found different degrees of population structure within genetic lineages, with hierarchical levels of between-population differentiation within the eastern lineage (in which Aranjuez arose as an evolutionarily distinct population), and absence of differentiation between the two western populations studied (cf. Díaz et al., 2017). Within this genetic framework, our analysis of the dorsal striped versus unstriped pattern, a heritable trait with adaptive value in this species , showed a phenotypic gradient with Aranjuez at one end, and the western populations at the other. The lizards from Lerma are phenotypically very similar to western lizards, a pattern that contradicts their F I G U R E 3 Genetic structure according to 73,291 SNPs (a) or 21 outlier SNPs (b, K = 4; c, K = 2). Bar colors represent posterior probabilities of cluster membership unmistakably eastern genetic background, but is consistent with their western-like habitat. This link between habitat and phenotype was in harmony with the pattern of genetic structuring of the outliers identified in the genome, an observation which substantiates the adaptive value of these genetic variants (Rhode et al., 2017;Rosenblum, Rompler, Schoneberg, & Hoekstra, 2010;Schield et al., 2017). Highly divergent populations (Aranjuez and western) have genetic variants subject to selection that seem to have been fixed along their evolutionary history. The high correlation between PG and PO in these populations further indicates that the position of individuals along the genetic gradient defined by selected alleles (outliers) corresponds to their location along the axis of genomic backgrounds opposing eastern and western lineages. In other words, it seems that the genetic variability that is adaptive in these populations (detected by ROA) is given by lineage sorting, and therefore, it could have been fixed along the long branches of the phylogenetic structure of our sample. However, the same correlation was not significant in F I G U R E 4 MDS analysis according to 73,291 SNPs that produced two major axes; the first of them consisted of a west-east gradient that discriminated between the two lineages, confirming the ancestral east-west division on a genome-wide basis Therefore, our results provide compelling evidence of ancestral genetic variation contributing to lizards' adaptability through its fixation along eastern and western lineages (Barrett & Schluter, 2008).
The process of local adaptation in Lerma seems to be favored not only by a genetic background consolidated through a long history of local adaptation to changing conditions, but also by the incorporation of genetic novelty (locally or imported from other populations; Barton & Etheridge, 2018). Gene flow among lizard populations was apparently influenced by ecological processes and not just by geographic proximity (Orsini, Vanoverbeke, Swillen, Mergeay, & De Meester, 2013;Rosenblum, 2006). For example, gene flow into ecologically extreme populations was asymmetric: Aranjuez was almost completely isolated while western populations scored the greatest levels of immigration. Aranjuez seems thus to be the place where incoming genotypes are more strongly purged out, probably because lizards occupy singular habitat of perhaps reduced quality. In forested habitat typical of western populations, the large Psammodromus prefers open forests where sunlit patches are readily available (Santos et al., 2008). In this habitat, genetically eastern lizards would not be selected against as strongly as the western ones may be purged out in shrubby habitat typical of eastern populations.
In addition, eastern lizards may find opportunities in western habitat, where higher vegetation cover results into lower detectability by visually oriented predators (even with their striped phenotype; Díaz et al., 2017), and they can benefit from higher food abundance in forests (Consuegra, Verspoor, Knox, & García De Leániz, 2005;Tenan, Fasola, Volponi, & Tavecchia, 2017). The opposite is not true: Western lizards find poor habitat in open scrubland (such as Aranjuez in our study; Díaz et al., 2017). Also, eastern (striped) lizards could be better adapted to dispersal through unforested, grassy habitat than western (unstriped) ones (Duckworth, 2008;LaRue, Holland, & Emery, 2018). Whatever the causes of the asymmetry in migration rates found in our study, unrestricted gene flow from western populations, along with ancestral genetic variation, would allow eastern populations (especially Lerma) to approach the western lizard phenotype in forested habitat (Alleaume-Benharira, Pen, & Ronce, 2006;Phillips, 1996).
Eastern and western lineages of the large Psammodromus have most likely remained locally adapted during the glacial cycles, maintaining their profound divergence but also great adaptive potential favored by the maintenance of ancestral genetic variation in the face of habitat heterogeneity (Barrett & Schluter, 2008;Gallet, Froissart, & Ravigné, 2018;Svardal, Rueffler, & Hermisson, 2015). Our outlier analysis (ROA) supports this idea, by documenting both the divergence of populations at the extremes of the environmental gradient, and the maintenance of high heterozygosity in eastern populations whose habitat has seemingly acquired western-like conditions. In addition, admixture may also have been important for maintaining populations adapted to changing local conditions. Thus, a variable permeability to immigration conditioned by a dynamic landscape configuration may have prevented speciation while maintaining strong genetic structure (Räsänen & Hendry, 2008;Whitney, Bowen, & Karl, 2018). This is remarkable for a terrestrial vertebrate with limited dispersal ability (Santos et al., 2008) after 3 millions of years of lineage divergence (Carranza et al., 2006). At the same time, the combination of divergence subject to selection, maintenance of ancestral adaptive variation, and permeability to gene flow may have rendered the large Psammodromus a species able to adapt to a wide variety of habitats (Lande & Shannon, 2006;Mackay, 1981;Paccard et al., 2018;Parter et al., 2008;Peniston et al., 2019;Welch & Jiggins, 2014), a trait which may underlie its success as the most abundant and widespread Iberian lizard (Huang et al., 2016;Pleguezuelos, Marquez, & Lizana, 2004 Levinst, 1962). After all, if only one dimension of environmental variation is considered, moving away from one end of the gradient inevitably leads to approaching the other end (Lahti et al., 2009).
All in all, the combined use of ROA and CIOA provided compelling evidence of different sources of genetic variation contributing to local adaptation dynamics in a lizard species. We are aware that this was feasible thanks to various singularities of our study system, with populations distributed over a wide, heterogeneous range that has been changing since the split of eastern and western lineages in the late Pliocene. In this biogeographical scenario, adaptive lineage splitting was counterbalanced by ancestral genetic variation and asymmetrical patterns of gene flow (likely promoted by historical or current environmental instability), which could hamper speciation despite strong genetic structure. Complex systems like this one may thus exemplify how to shed light on a variety of important questions at the interface between local adaptation and the first steps toward speciation (or the lack of it; Rosenblum, 2006;Rundle & Nosil, 2005;Schluter, 2001;Via, 2009). Overall, the combined use of ROA and CIOA may offer a powerful tool to elucidate how evolution provides the genetic variation that ecology demands throughout the history of a species. and an anonymous referee improved the paper with their comments.

CO N FLI C T O F I NTE R E S T
None declared.