Discussion
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, 2010, 2011, 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 individuals (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, Llanos-Garrido et al. 2017). 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, 2010, 2011, 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, the background genetic differentiation between these two populations is very low (Díaz et al. 2017, Verdú-Ricoy et al. 2010), even when analytical power is increased with genomic scans. Irrespective of phenotypic differenciation, other populations of Psammodromus algirus showed marked genetic isolation despite being much closer (both geographycally and environmentally) than the ones presented here (Díaz et al. 2017, Llanos-Garrido et al. 2019). 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 (Santos et al. 2008, Pérez-Tris et al. 2019). 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 and Keitt 2000, Takahashi et al. 2016). Such decrease in genetic diversity and/or effective population size could have led to a subsequent reduction of migration rate and increased genetic structure in this marginal population (Holt and Keitt 2000, Zalewski et al. 2016, Li et al. 2018; see Deng et al. 2020), 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 any hint of genetic structure, it is nearly impossible to properly estimate gene flow among undefined populations (Pfenninger and 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.
However, the lack 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 (Philips et al. 2008, Moody et al. 2015, Poelstra et al. 2014, Shaner et al. 2015). In our case, it seems that the divergence between populations is so recent 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 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 in our system, because background divergence is very low across all the genome, which suggests that SNPs with almost any degree of divergence could have been detected as outliers.
Finally, there could be no significant divergence at all on a genome-wide scale. It is true that previous studies point to a genetic basis, especially those that included reciprocal transplants or common garden experiments, but it could be the case that the observed differences were the result of a heritable epigenetic response to the environment (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 y 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 any detectable genotypic isolation between them. Moreover, and for reasons that remain unveiled, genome-wide differenciation 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 2009). Highlighting natural systems where divergent adaptation occurs with different degrees of genomic-wide differenciation (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 2009, Shafer and Wolf 2013, Sexton et al. 2014, Krohn et al. 2018).