Low genome‐wide divergence between two lizard populations with high adaptive phenotypic differentiation

Abstract Usually, adaptive phenotypic differentiation is paralleled by genetic divergence between locally adapted populations. However, adaptation can also happen in a scenario of nonsignificant genetic divergence due to intense gene flow and/or recent differentiation. While this phenomenon is rarely published, findings on incipient ecologically driven divergence or isolation by adaptation are relatively common, which could confound our understanding about the frequency at which they actually occur in nature. Here, we explore genome‐wide traces of divergence between two populations of the lacertid lizard Psammodromus algirus separated by a 600 m elevational gradient. These populations seem to be differentially adapted to their environments despite showing low levels of genetic differentiation (according to previously studies of mtDNA and microsatellite data). We performed a search for outliers (i.e., loci subject to selection) trying to identify specific loci with FST statistics significantly higher than those expected on the basis of overall, genome‐wide estimates of genetic divergence. We find that local phenotypic adaptation (in terms of a wide diversity of characters) was not accompanied by genome‐wide differentiation, even when we maximized the chances of unveiling such differentiation at particular loci with FST‐based outlier detection tests. Instead, our analyses confirmed the lack of genome‐wide differentiation on the basis of more than 70,000 SNPs, which is concordant with a scenario of local adaptation without isolation by environment. Our results add evidence to previous studies in which local adaptation does not lead to any kind of isolation (or early stages of ecological speciation), but maintains phenotypic divergence despite the lack of a differentiated genomic background.

cause local adaptation (Funk et al., 2011;Nosil et al., 2008;Schluter, 2000). In theory, isolation by adaptation should be theoretically detectable in the form of divergence peaks among locally selected loci (or among loci genetically correlated with them), but if populations remain isolated long enough, selectively neutral loci can also become differentiated, which makes the detection of outlier loci empirically challenging (Krohn et al., 2019;Llanos-Garrido et al., 2019a).
However, discordances may also appear, whereby two populations may be genetically undifferentiated while showing evident phenotypic divergence (Moody et al., 2015;Palmer & Kronforst, 2015;Shaner et al., 2015). Such discordances can happen, despite the existence of processes that blur overall genomic divergence (e.g., gene flow, recent divergence), only when natural (or sexual; Yang et al., 2018) selection is strong enough to overcome such processes, favoring the presence of locally divergent regions with variants under selection within an otherwise undifferentiated overall genomic background (Burri, 2017;Wang et al., 2019).
While numerous studies have dealt with locally adapted populations occupying different environments or isolated by ecological barriers (Orsini et al., 2013;Rosenblum, 2006;Sexton et al., 2014;Zhao et al., 2020), only a few have focused on the lack of correlation between isolation and adaptation (Feder et al., 2014). Moreover, it has been suggested that there is a significant publication bias against such studies (Krohn et al., 2019). This possible under-representation of apparent negative results may lead to biased estimates of how frequently local adaptation occurs without isolation, and this is precisely the reason why replicated studies on species or populations with different degrees of isolation are needed (Feder et al., 2014;Sendell-Price et al., 2020;Talla et al., 2017). In fact, if studies that do not find any effect of environmental gradients on genetic differentiation are rarely published, examples of incipient ecological speciation and/or isolation by environment may artificially be deemed frequent or even widespread (Hendry et al., 2009;Sexton et al., 2014;Shafer & Wolf, 2013).
The aim of this study is to elucidate the patterns of genetic differentiation that underlie phenotypic divergence between two populations of the lacertid lizard Psammodromus algirus separated by a 600-700 m altitudinal gradient. This gradient is associated with a large number of habitat differences, including forest type (deciduous vs. perennial) or average annual rainfall (1,170 vs. 438 mm; see the Methods section for a detailed explanation of these habitat differences). The analysis of mitochondrial DNA sequences has shown that these populations present very little genetic differentiation Verdú-Ricoy et al., 2010), even though they differ in a wide variety of adaptive phenotypic characteristics, including many life history traits (Iraeta et al., 2006(Iraeta et al., , 2011(Iraeta et al., , 2013, sexual ornamentation (Iraeta et al., 2011;Llanos-Garrido et al., 2017), and prevalence of ticks , among others (Table 1).
The evidence for such adaptations is based on previous studies that have shown, through reciprocal transplant and common garden TA B L E 1 Phenotypic differences found between the populations of El Pardo and Navacerrada in previous studies

Phenotypic trait Reference
Shorter incubation times in the montane population Iraeta et al. (2006) Larger hatchlings in the lowland population Iraeta et al. (2006) Faster growth in lowland juveniles (reciprocal transplant experiment) Iraeta et al. (2006) Faster growth in lowland juveniles (common garden experiment) Iraeta et al. (2013) More plastic activity levels in response to food availability in the lowland population Iraeta et al. (2008) Larger clutches of smaller eggs in the lowland population Iraeta et al. (2008Iraeta et al. ( , 2013 Longer flight distance for pregnant females in the montane population Iraeta et al. (2010) Relatively longer hind limbs in the lowland population Iraeta et al. (2011) Larger adult females in the montane population Iraeta et al. (2013) Countergradient variation in body size: the genotypes that presumably control the key adaptations of the lowland population (larger eggs and hatchlings, and faster growing juveniles) occur in a low-productivity environment in which lizards grow more slowly and reach a smaller adult size Iraeta et al. (2006), Iraeta et al. (2013) More and relatively larger femoral pores in males from the lowland population Iraeta et al. (2011) Greater development of the sexual coloration of the head (i.e., larger colored surface) in males from the lowland population Iraeta et al. (2011 Increased saturation of the sexual coloration of the head in males from the montane population Llanos-Garrido et al. (2017) Males respond to the activation of the immune system by reducing the extent of the sexual coloration of the head in the lowland population, and its saturation in the montane population Llanos-Garrido et al. (2017) Higher rates of infestation by tick nymphs in the montane population (no ticks in the lowland population) Llanos-  experiments, that these adaptive phenotypic differences are not sustained solely by environmental effects (Iraeta et al., 2006(Iraeta et al., , 2013. Therefore, there must be a genetic basis to determine such phenotypic differences between these apparently undifferentiated populations, even if such genetic basis does not lead to an environmentally based isolation. To define the degree of genetic differentiation between the two populations, we used a GBS genomic scan based on 73,291 SNPs that allowed us to analyze the genetic structure and distance between them. In addition, we used a Bayesian method of detection of SNPs with genetic distances between populations greater than expected given the degree of background genomic differentiation (i.e., F ST -based outlier test; Bayescan: Foll and Gaggiotti, 2008). With this approach, we tried to define polymorphisms possibly associated with the patterns of divergent selection that promote the observed adaptive differences (Bonhomme et al., 2010). Thus, if local adaptation has led to isolated populations, we expect to find a genome-wide differentiation standard that could hamper the detection of divergence peaks at adaptive loci. Alternatively, we might find an undifferentiated genomic background where such peaks should be easy to detect, at least in theory.
We used an approach in which all genetic variants are used to infer the basal level of genomic differentiation, defining outliers as SNPs with a greater degree of divergence and with allelic frequencies deviated from the expected under neutral selection (Lewontin & Krakauer, 1975). The reason why comparing differentiation peaks with adjacent regions should facilitate this task is that the degree of differentiation is heterogeneously distributed throughout the genome (Campagna et al., 2015). Therefore, detecting a divergence peak in a specially conserved region is challenging using delocalized genetic variation, which may find such degree of divergence even below the background genomic differentiation, but which may actually be very divergent within its genomic location (Lawson and Petren, 2017). This scenario is especially common in the coding regions where the genetic basis of phenotypic diversity is located, so that approaches with delocalized SNPs do not usually respond to what is the genetic basis of a given phenotype, but describe the general patterns of genetic differentiation that lay behind the process of local adaptation (e.g., Llanos-Garrido et al., 2019a;Tigano et al., 2017). Thus, the aim of this study is not to uncover the specific genetic basis behind already published adaptations, but to elucidate whether such local adaptation is accompanied by a significant degree of genome-wide differentiation or not (Krohn et al., 2019).

| Study system
Psammodromus algirus (Linneus, 1758) is the most abundant lacertid lizard in the Mediterranean scrubland and forest habitats of the Iberian Peninsula (it is only absent from the northern Eurosiberian area). It is a medium-sized species (snout-vent length: 60-90 mm; body mass: 6-16 g) in which variation in body size is an important indicator of fitness, as it determines survival, male attractiveness, and female fecundity (Díaz, 1997;Díaz et al., 2005;Iraeta et al., 2013;Martıń & Forsman, 1997). Psammodromus algirus males present a reddish nuptial coloration during the months of highest activity (April-June) that can range from small sublabial marks in subdominant males to occupy practically the entire head in the largest males. This visual signaling is complemented by the existence of femoral pores as a means of chemical signaling. The number of femoral pores is higher in males and their activity increases as the breeding season advances in response to increased blood androgen levels (Chiu & Maderson, 1975;Cole, 1966;Diaz et al., 1994).  Table 1).
The montane population is located at Navacerrada ( in both locations, but reaches higher densities at higher altitude (Díaz, 1997). the area of analysis using the "magnetic loop" tool, we measured the red-colored surface with the "magic wand" tool (with 30% tolerance) after selecting a random red-colored point, and we then used the "similar" option of the magic wand, with the same tolerance, to select all areas with a similar coloration. The colored surface was measured as the percentage of colored pixels in the area of analysis.

| Field sampling and phenotypic analysis
To calculate the red saturation, we use the proportion of the red channel within the RGB channel, that is: R/(R + G + B), where R, G and B are the red, green, and blue channels of the graphics card. All these measures were taken blindly with respect to the population of origin. Once all this was done, the lizards were released on their site of capture.
All phenotypic analyses were performed with Statistica software (Statsoft). In order to analyze between populations differences, we used two-way ANOVAs with sex and population as factors and with different morphological measures as the dependent variables. When necessary, we used ANCOVAs to control for the effect of body size. To assess the relationship between tick load and different explanatory variables (e.g., sex, size, or number of femoral pores), analyses were restricted to the Navacerrada population. To summarize variation in male size in this population, we ran a Principal Component Analysis (PCA) that produced a single factor explaining 84.2% of the variation in the data matrix and giving high factor loadings to all size variables (snout-vent length, head length and width, mean leg length, and body mass; Table 2).

| DNA extraction and sequencing
We obtained tissue samples by clipping 2 cm of the tail tip of lizards, and we purified DNA for library preparation using the Biotools Speedtools Tissue DNA Extraction kit. One of the samples (from the El Pardo lowland population) lacked sufficient DNA to be included in our analyses, leaving a final sample size of 39 genotyped lizards.
We used the restriction enzyme Pst1 to prepare the libraries for genome-wide genotyping with GBS (Genotyping by Sequencing).
The sequencing was done in an Illumina HiSeq2500 sequencer, and once the sequences were obtained, we characterized single nucleotide polymorphisms (SNPs). For this purpose, we used the UNEAK pipeline, implemented in TASSEL v.5.0 (Bradbury et al., 2007), and specifically designed for samples of species without a reference genome. The sequences were pruned to eliminate sequencing errors using the error rate threshold parameter. The resulting database consisted of 83,648 SNPs with a coverage of 5.68 ± 6.56 (mean ± SD) and missing site rate of 0.49 ± 0.33. We then discarded loci with minor allele frequencies <0.05, and those that were sequenced in less than 10% of individuals. After this pruning process, the final database consisted of 73,291 loci, with a single SNP per locus and a coverage of 6.6 ± 6.75 and a missingness of 0.42 ± 0.31.
This dataset was the one used to analyze the neutral variation of our samples.

| Characterization of neutral variation and SNPs under selection
Before performing the genetic structure analysis, we used PLINK v1.9 software (Purcell et al., 2007) to prune the SNPs according to We also applied another filtering step before doing the outlier analysis to minimize the false positive rate. We discarded loci whose minor allele frequencies were <0.05 in any population (thus excluding privative alleles) and loci that could not be sequenced in at least 75% of the individuals in each population. The resulting database consisted of 6,421 SNPs. To identify SNPs putatively under selection, we performed an outlier analysis with Bayescan v.2.1. (Foll and Gaggiotti, 2008), a very conservative method which is not prone to false positives, and is very useful when the number of populations is low (Foll and Gaggiotti, 2008). This program uses a logistic regression to split the F ST coefficients into a population-specific effect (β) and a locus-specific effect (α). We selected loci with α > 0, suggesting positive selection, and a false discovery rate (corrected by multiple testing) of q < 0.05.
Once this was done, we tried to annotate SNPs with the high-

| Environmental analysis
To evaluate environmental variation between El Pardo and Navacerrada (Figure 1)

| Phenotypic variability
We found significant differences between populations that confirmed the results of previous studies: lizards from the lowland population were smaller and had relatively longer legs than those from the montane population; in addition, they never presented tick infestation, despite adult ticks being present in their habitat (Table 3). Sexual coloration data were especially scarce from males sampled in this study (few colored males with very weak colorations), and as a consequence differences between populations could not be analyzed. When we only considered lizards from Navacerrada, we found that males had more ticks (F 1,18 = 10.85, p = .004) and that only the largest males could develop a sexual coloration (F 1,8 = 42.17, p < .001; Figure 2).

| Neutral genetic variability and SNPs under selection
The analysis of genetic structure with admixture yielded a single more plausible model that grouped all samples into a single cluster (K = 1, cv-error = 0.382), indicating no apparent genetic structure between both populations. Even when forcing the model toward higher values of K (i.e., less plausible results), we still failed to find any signs of between-populations genetic structure (K = 2, cverror = 0.405; Figure 3). A PCA on the covariance matrix obtained from the same data set returned a scree plot that showed a sharp decrease in eigenvalues between PC1 and PC2 ( Figure S1), indicating that PC1 was the single axis that best captured overall genetic variation. One-way analyses of variance showed that both populations had overlapping scores ( Figure S2) with means that were nearly identical for PC1 (F 1,37 = 0.000018, p = .997) and did not differ significantly for PC2 (F 1,37 = 1.08, p = .305).
Regarding outliers analyses, Bayescan did not detect any locus subject to selection (q-values range: 0.883-0.910; mean F st ± SD = 0.023 ± 0.0003). However, we explored the loci with largest values of α, the magnitude of the locus-specific effect of selection, in case they were biologically meaningful. The largest increase in α was concentrated on four maximum (although low) values, two of which could not be annotated at all (Table 4) (Table 4).

| Environmental variability
By using the scree plot criterion, the environmental PCA yielded a  TA B L E 3 Phenotypic differences between the populations studied in this work: results of one-way ANOVAs with population as factor, and of two-way ANCOVAs (for the last two variables, marked with an asterisk) with population and sex as factors and size (SVL) as the covariate (mean = −1.563, SD = 0.725; Figure 1). In fact, such variation was much greater than the one found in the geographical cells within any of the two populations, either lowland (mean = −1.069, SD = 0.048; Levene test: F 1,28 = 12.048, p = .002) or montane (mean = −2.121, SD = 0.043; Levene test: F 1,42 = 25.807, p < .001).

| D ISCUSS I ON
Our results point to a discordance between the absence of genomic divergence and consistent phenotypic differentiation in traits whose adaptive value seems obvious (Iraeta et al., 2006(Iraeta et al., , 2010(Iraeta et al., , 2011(Iraeta et al., , 2013.
The fact that, as previously found, in the montane population the load of ectoparasites is related to secondary sexual characteristics of males, and that there are no infestations in the lowland population, shows that the presence of ticks is associated with different selective pressures in both populations (Llanos-Garrido et al.,

2017). Apart from their negative effect on the survival of individu-
als (Salvador et al., 1996), ticks appear to have other significant impacts on fitness, influencing the expression of sexual ornaments that are important for reproduction Salvador et al., 1996). In addition, the altitudinal gradient involves many other selective pressures that, while less easy to identify than the one imposed by ticks (but still evident given the environmental differences shown by our PCA scores), generate an undeniable impact on fitness-related traits such as growth rate, egg size, or clutch size (Iraeta et al., 2006(Iraeta et al., , 2010(Iraeta et al., , 2011(Iraeta et al., , 2013. In addition, the adaptive phenotypic differences found in this and earlier studies span the whole life cycle of lizards (e.g., incubation time, hatchling size, juvenile growth rate, and sexual ornamentation and reproductive investment of adults), which should involve a considerable variety of divergent selective pressures underlying their phenotypic differentiation (Roff, 2007). Moreover, the fact that these traits maintain their differentiation in common garden conditions (Iraeta et al., 2013) or in reciprocal transplant experiments (Iraeta et al., 2006) reveals a divergent genetic basis between the two populations (De Kort et al., 2014). However, background genetic differentiation between these two populations is very low Verdú-Ricoy et al., 2010), even when analytical power is increased with genomic scans.
Irrespective of phenotypic differentiation, other populations of P.
algirus showed marked genetic isolation despite being much closer (both geographically and environmentally) than the ones presented here Llanos-Garrido et al., 2019a). This is true even at a fragmented metapopulation where suitable habitat patches are as close as hundreds of meters apart, and where subpopulations have diverged less than 100 years ago (Pérez-Tris et al., 2019;Santos et al., 2008). However, this latter scenario occurs at the northern edge of the species' range, where genetic diversity might be lower than at El Pardo or Navacerrada, which occupy the center of the Iberian distribution range (Holt & Keitt, 2000;Takahashi et al., 2016). Such decrease in genetic diversity and/or effective population size could F I G U R E 3 Probability of assignment to groups generated by Admixture for the model with two groups (K = 2). On the left, CV-errors for each of the models analyzed; the most credible model is the one that includes only one group (K = 1) have led to a subsequent reduction of migration rate and increased genetic structure in this marginal population (Deng et al., 2020;Holt & Keitt, 2000;Zalewski et al., 2016), whereas the opposite might be true at the center of the species' distribution range. However, as migration rate between El Pardo and Navacerrada is high enough (or divergence time recent enough) to blur overall genetic structure, it is nearly impossible to properly estimate gene flow among undefined populations (Pfenninger & Posada, 2002). Thus, elucidating how gene flow affects local adaptation in this system constitutes a distinct challenge, really difficult to circumvent even by increasing our analytical power or by getting longer assemblages.

TA B L E 4 SNPs with greater divergence between populations detected by Bayescan
Nevertheless, a low degree of overall genetic structure should not hinder the genetic differentiation of specific variants on which the basis of adaptive divergence is located. There are other systems in which the background genetic differentiation between phenotypically well-differentiated groups is low, and yet islands of selection whose local divergence outstands above the rest of the genome can still be found (Aguillon et al., 2018). These genome fractions may persist differentiated even in continuous gene flow scenarios between the study groups (Moody et al., 2015;Poelstra et al., 2014;Shaner et al., 2015). In our case, it seems that the divergence between populations is so low that it has not translated into a genomic imprint to a scale large enough to be detected with the approach we have used (Aguillon et al., 2018). However, there is compelling evidence that points to the existence of at least a part of the genome that has diverged, giving rise to the observed phenotypic differences. One possible explanation is that the only divergent regions are those expressed as phenotypes (Safran et al., 2016). These regions would be located in very specific areas of the genome and could be really scarce (Campagna et al., 2015), to the extent that low coverage genotyping approaches could miss them.
Another possibility is that the SNPs responsible for the observed adaptive divergence could be located in regions of the genome with low genetic diversity, which would express a peak of divergence at a local genomic scale, but which would escape the outliers detection analysis used in this study (Campbell et al., 2018) if the amount of genetic variability is distributed heterogeneously throughout the genome (Feulner et al., 2015;Nosil et al., 2009;Poelstra et al., 2014). However, such heterogeneity does not seem to be an issue  (Groot et al., 2018;Trerotola et al., 2015). Such epigenetic marks could be inherited only up to a specific number of generations, after which the signal would be lost (Trerotola et al., 2015). To discard this hypothesis, we could confirm previous findings by extending common garden experiments a few generations more (Groot et al., 2018), or we could perform a transcriptomic differentiation study (McGirr & Martin, 2018).
In summary, we succeeded to demonstrate that genome-wide differentiation was not associated with local adaptation in an altitudinal gradient with marked environment differences between two lizard populations. The phenotypic differentiation pattern of these populations was obviously concordant with a scenario of ecologically driven divergence, but it did not lead to overall genotypic isolation between them. Moreover, and for reasons that remain unveiled, genome-wide differentiation seems much lower than expected by geographic distance. Thus, we uncovered another example of locally (and divergently) adapted populations, which remain undifferentiated on a genome-wide scale, a phenomenon that could be more frequent than previously thought due to publication bias against this apparent lack of positive results (Hendry et al., 2009). Highlighting natural systems where divergent adaptation occurs with different degrees of genomic-wide differentiation (including none) is of paramount importance to understand how frequently isolation by adaptation occurs, and to what extent the first incipient stages of ecological speciation are widespread (Feder et al., 2012;Hendry et al., 2009;Krohn et al., 2019;Sexton et al., 2014;Shafer & Wolf, 2013).

This research was funded by the Spanish Ministry of Science and
Technology (grants CGL2013-41642-P and CGL2017-82117-P to JPT, and PID2019-108341GB-I00 to JAD).

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