Deforestation limits evolutionary rescue under climate change in Amazonian lizards

The impact of climate change on biodiversity is often analysed under a stable evolutionary perspective focused on whether species can currently tolerate warmer climates. However, species may adapt to changes, and particularly under conditions of low habitat fragmentation, standing adaptive genetic variation can spread across populations tracking changing climates, increasing the potential for evolutionary rescue. Here, our aim is to integrate genomic data, niche modelling and landscape ecology to predict range shifts and the potential for evolutionary rescue.


| INTRODUC TI ON
Accelerated habitat degradation and climate change can elevate extinction rates in natural populations to unprecedentedly high levels (Andermann et al., 2020;Neubauer et al., 2021).Most predictions across different ecosystems and biological groups indicate extreme loss of species ranges and species diversity, especially when accounting for the interactions between climate change and habitat degradation (Feeley & Rehm, 2012;Newbold et al., 2019).However, some species might have populations adapted to different climate conditions across their ranges (Carvalho et al., 2011;Millien et al., 2006), including demographic and genetic factors able to potentially avoid local extinction (Forester et al., 2022) and enable their survival in future climates (Diniz-Filho & Bini, 2019).Therefore, climate change forecasts for the future of biodiversity are expected to be improved and more accurate by integrating intraspecific genetic information with the eco-evolutionary dynamics among populations (e.g.adaptation and dispersal) (Bothwell et al., 2021).
Traditionally, species ranges are the primary source of data used for predicting biodiversity responses under changing climates (Bellard et al., 2012).Most attempts to predict range shifts have used ecological niche modelling (ENM) based on species records and environmental data alone (Wiens et al., 2009).However, these methods generally assume that species will move across less suitable environments to reach newly suitable areas or become locally extinct in the parts of the current range predicted to become unsuitable (Guisan et al., 2017).In this regard, population genomics can help estimate the tolerance and evolutionary potential of species in the face of climate change or the potential for evolutionary rescue (Balkenhol et al., 2016;Gotelli & Stanton-Geddes, 2015;Waldvogel et al., 2020).Evolutionary rescue refers to the process by through which populations avoid extinction through rapid adaptation to a changing environment (Bell, 2017).This process occurs not only through novel adaptations but also when standing adaptive genetic variation enables the persistence of populations, especially when individuals with climate-adapted genotypes can disperse across the landscape and reach maladapted populations (Bell, 2017;Razgour et al., 2019).The latter case is often also called genetic rescue; however, this term is generally related to the deliberate introduction of new genetic variants via human intervention to improve the fitness of a population (Hoffmann et al., 2021).Genotype-environment association analyses can be used to identify climate-driven adaptive genetic variation (Rellstab et al., 2015).Information resulting from these approaches can then be integrated with ENMs for future habitat suitability and landscape structure models that include more biologically realistic parameters (e.g.filters or barriers to dispersal) to improve predictions of range shifts under changing climates and habitats in species and natural populations (Razgour et al., 2019).
The effects of climate change and habitat loss on biodiversity are predicted to be highly detrimental in South America (Sales et al., 2020), particularly for ectothermic organisms such as lizards (Sinervo et al., 2010).Temperatures in the tropical parts of the continent may increase by 2.5-4.5°C by the end of this century (IPCC, 2021), resulting in intensified droughts and, possibly, the conversion of forest to savannas, and savannas to xeric scrublands (Cooper et al., 2020;Zemp et al., 2017).These changes may also interact with the ever-increasing deforestation rates (Staal et al., 2020), leading to increased extinction rates and a lower potential for evolutionary rescue (Razgour et al., 2019).Habitat degradation may lead to the amplification of climate change (e.g.increasing edge effects and fires occurrence and intensity) and is already affecting the world's largest and most biodiverse rainforest, the Amazonia, and most biodiverse savanna, the Cerrado (Souza et al., 2020).However, we still do not know to what extent natural populations and species in such tropical biomes might be resilient to anthropogenic changes.
Reptiles are suitable models for testing the effects of climate change and habitat degradation on tropical biodiversity.Most species have relatively lower dispersal abilities than birds or large mammals (Azevedo et al., 2021;Saladin et al., 2019), meaning that they are more likely to respond to local changes.Also, reptiles are ectothermic, which means that climate changes are more likely to affect local populations and cause range shifts since their metabolism is directly related to the environmental temperature (Azevedo et al., 2021;Huey, 1982).Lizards are conspicuous elements of the South American biological communities (Gasnier & Magnusson, 1994).Nevertheless, they face elevated thread levels due to climate change (Sinervo et al., 2010) and land use (Palmeirim et al., 2017).
Many Amazonian lizards are predicted to experience local extinction or demographic reduction due to climate change (Diele-Viegas et al., 2019).For instance, populations of the whiptail lizard Kentropyx calcarata (Squamata: Teiidae) present elevated risks of local extinction when considering data on thermal physiology (Pontes-da-Silva et al., 2018).This vulnerability is high at the ecotone between the Amazonian rainforest and the Cerrado savanna, coincident with the so-called 'Deforestation Arc' most impacted by human activities in Amazonia (Albert et al., 2023;Costa & Pires, 2010; until 2070 could foster evolutionary rescue of ectothermic organisms.These actions could prevent substantial biodiversity loss in Amazonia, emphasizing the importance of understanding species adaptability in maintaining biodiversity.

K E Y W O R D S
dispersal, ecotones, genome-environment association analysis, neotropics, niche modelling, range shifts et al., 2020;Rehm et al., 2015).The geographically structured thermal responses identified so far suggest local adaptation to the climate in certain populations, even though Amazonian Kentropyx species seem to have a high degree of climate niche conservatism across deeper timescales (Sheu et al., 2020).Kentropyx calcarata is distributed across Amazonian lowlands and the northern Atlantic Forest, restricted to forested habitats (Ribeiro-Junior & Amaral, 2016).Understanding the widespread distribution and habitat specificity to forests in K. calcarata may provide general insights into the impacts of climate change and deforestation on biodiversity.

Marques
Here, we investigate range shifts and the potential for evolutionary rescue in future climate scenarios (Figure 1) after determining the neutral genetic structure and climate-driven adaptive genomic variation in K. calcarata.We expect adaptive genetic variation to be structured across Amazonia, as in thermal traits of the species (Pontes-da-Silva et al., 2018).We also expect that local extinction forecasts will be higher across the Amazonia-Cerrado ecotone due to the higher deforestation rates in that region relative to more central parts of Amazonia (Silva Junior et al., 2021).We predict current and future distributions of climate-adapted individuals by directly integrating ENMs and dispersal constraints with deforestation predictions in the region.Our approach provides a spatial assessment of the potential for evolutionary rescue in Amazonian ectothermic organisms, yielding insights on local climate adaptation, with implications for biodiversity management and conservation of this keystone tropical forest in the face of climate change.

| Summary
To investigate adaptive selection to climate and the potential for evolutionary rescue in K. calcarata across Amazonia, we used genomic data derived from RAD-sequencing.After controlling for neutral genomic variation, we identified distinct genotypes adapted to different parts of the climate gradient.Using environmental niche modelling, we predicted distributions of each genotype until 2100, considering deforestation as a potential barrier to dispersal and a cause of habitat loss.Finally, we estimated the potential for evolutionary rescue of each genotype by analysing the overlap between the areas with range loss of one genotype and the areas of permanence or expansion of the other genotype.

| Genetic sampling, processing and population structure
We sampled 112 individuals of all Kentropyx species within the known range of K. calcarata in Amazonia, the Atlantic Forest and the Cerrado savannas (information on specimens, museums and localities are available in Table S1).
We extracted genomic DNA from each sample (liver and muscles) using the Macherey-Nagel® Mini Kit with a high salinity protocol F I G U R E 1 Graphical representation of our approach for inferring areas with potential for evolutionary rescue.Potential for evolutionary rescue in individuals with genotypes adapted to distinct environmental conditions in a hypothetical species.Note that due to dispersal limitations, not all suitable areas are occupied by the respective genotypes (represented as asterisks and crosses) as represented in the environmental gradient in the top figures.Below, the current and future (2100) distributions are overlapped.Then, areas in which Genotype 2 becomes extinct but which are occupied by individuals with Genotype 1 (either by dispersal or permanence) are considered areas with potential for evolutionary rescue through the spread of climate-adapted genomic variation.The drawing represents a Kentropyx calcarata.
with proteinase K.We then visualized fragment sizes in agarose gels and measured DNA concentration and quality with Qubit™ 3.0 Fluorometer and Nanodrop, selecting only samples with a low degree of DNA fragmentation and normalizing the concentration to 20 ± 2 ng DNA for each 50 μL solution (TE buffer).RAD-sequence library preparation (as detailed in Etter et al., 2011) was then performed by Floragenex, Inc. (Eugene, OR, USA).Individual samples were digested with the SbfI enzyme, linked with barcoded RAD ligators, multiplexed, sonicated and then size selected to a range within 300-500 bp.After PCR amplification, DNA sequencing was performed on a 2 × 100 bp Illumina HiSeq platform.This procedure resulted in ~8,400,000 reads per individual with a length of 100 base pairs.We deposited the demultiplexed raw-sequencing data in the Sequence Read Archive (PRJNA1111372).
We demultiplexed and assigned the reads for each sample and performed a de novo assembly and SNP calling using iPyrad 0.9.55 (Eaton & Overcast, 2020), checking the quality of demultiplexed samples using MultiQC (Ewels et al., 2016).For the iPyrad pipeline, we allowed one mismatch from individual barcodes, clustering the reads of each sample (minimum size of 35 base pairs) and across the samples (de novo assembly), sequence coverage of 6× and minimum Phred score of 33.We selected a minimum clustering threshold of 0.90 after testing values from 0.85 to 0.95 for missingness and number of recovered SNPs (McCartney-Melstad et al., 2019).We allowed a maximum of two alleles and a maximum proportion of 0.5 for heterozygous sites per locus, including only loci present in 75% of the individuals.
From the pipeline above, we obtained 139,205 SNPs used in a principal component analysis (PCA) in iPyrad to identify samples clustered with other K. calcarata samples or with the remaining species.To check whether these samples would be similarly grouped regardless of the approach, we additionally inferred a species tree using SVDquartets (Chifman & Kubatko, 2014) with the tetrads function of iPyrad in 100 bootstraps.tetrads uses a multi-species coalescent method from quartet trees inferred from unlinked SNPs.After these steps, 66 samples across 33 localities with 348,109 SNPs (~20% missing sites) were consistently assigned to K. calcarata using both methodologies and used in the following steps (Figure S1).We then used VCFtools v. 0.1.16(Danecek et al., 2011) to filter SNPs with minor allele frequency (MAF) higher than 0.05 (Ahrens et al., 2018).
We also selected only one SNP per RAD stack (cluster of loci derived from an RAD de novo assembly) to decrease the probability of sampling linked SNPs, thus minimizing linkage disequilibrium (scripts provided by Prates et al., 2018).We kept only SNPs with less than 20% of missing data, obtaining a final data set of 30,589 SNPs for downstream analyses of population structure and genome-environment association.
We inferred the number of genetic clusters (k), population structure and degree of admixture for the 66 K. calcarata individuals using sparse non-negative matrix factorization (sNMF) (Frichot et al., 2014).We calculated least-squares estimates of ancestry proportions for k ancestral populations varying from one to 10 and evaluated model fit with an entropy criterion based on 100 crossvalidation repetitions.Next, we chose the value of k with the lowest entropy, with a 5% tolerance error (Frichot et al., 2014).Additionally, we produced a PCA from the same SNP data and used a scree plot to identify the number of significant components.These analyses were implemented in the R-package LEA 2.0 (Frichot & François, 2015).
We selected the following bioclimatic variables for the subse- dispersal constraints (cellular automata).We used the available moderate (AB1) and extreme (A2) scenarios, which have previously been used to model dispersal and range loss in South America (Sales et al., 2020) and closely match the patterns of habitat loss in Amazonia.We transformed the five available landscape classes (forest, grassland, farmland, urban and barren) into the two relevant classes for our study: forest (habitat of K. calcarata) and non-forest (all remaining categories-Figure S2), the last considered as barriers for dispersal in the downstream analyses.These non-forest areas also include some of the largest Amazon Basin rivers as well (when wider than 1 km-the resolution of the layer).The forest/barrier layer for the present time was used for delimiting current ranges by cropping the ENMs predictions for current climates (see below).We used the future projections to assess corridors and barriers for dispersal in range projections until 2100.
We conducted genome-environment association analyses to identify signatures of climate-driven genetic variation based on associations between allele frequencies and local climates (Rellstab et al., 2015), using the four selected bioclimatic variables, the filtered SNPs and accounting for the underlying neutral genetic structure (sNMF and PCA, as described above).We employed two approaches to identify candidate SNPs with a signature of climate-driven genetic variation, selecting those recovered in both methods for downstream analyses (Forester et al., 2018;Razgour et al., 2019).First, we used redundancy analysis (RDA) for detecting outlier loci, that is, candidate SNPs under climate selection.RDA enables the simultaneous detection of candidate SNPs across multiple environmental predictors, showing how groups of SNPs covary in the multivariate space, even for weak associations from each locus (Rellstab et al., 2015).We used the first three axes of the PCA (>80% of the variation) calculated from the K. calcarata SNPs (Forester et al., 2018) as covariates to control for the neutral genetic structure.We followed Capblancq and Forester (2021) The second genome-environment association method we implemented was a latent factor mixed model with regularized leastsquares using ridge penalty (LFMM 2- Caye et al., 2019).LFMM screens each SNP for signatures of local adaptation using the neutral genetic structure (k value derived from sNMF) as latent unobserved variables (Frichot et al., 2013).Unlike RDA, LFMM is implemented against only one environmental variable each time.
Therefore, we ran LFMM for all SNPs against the first three PCA axes, which explained >80% of the variation of the same bioclimatic variables used for RDA (Capblancq et al., 2018).We performed additional runs with k−1 and k+1, selecting the run with the lowest genomic inflation factor (λ).We then corrected the pvalues by the genomic inflation factor while ensuring a uniform distribution of p-values along the interval from 0 to 1, with a peak of values close to 0 (candidate loci), visualized through histograms (François et al., 2016).As for the RDA procedure, we considered loci with p-values inferior to .01 divided by the number of SNPs as candidate SNPs.

| Climate-adapted genotypes, modelling ranges and evolutionary potential
Each genome-environment association method may detect a set of candidate SNPs that do not entirely match, potentially leading to incongruences when allele mapping frequencies (Ahrens et al., 2018).We used the most conservative approach of retaining for the downstream analyses only the SNPs selected in both RDA and LFMM (Forester et al., 2018).However, we also checked for the spatial congruence of the allele frequencies for all possible sets of candidate SNPs recovered (RDA only, LFMM only, RDA-LFMM intersection or RDA-LFMM union), highlighting any potential differences in downstream analyses.
We reran the RDA including only each set of candidate SNPs to maximize the spread of the adaptive genetic variation across the environmental space (adaptively enriched RDA) (Capblancq & Forester, 2021).To aid in visualizing the degree of spatial congruence among the distinct sets of candidate SNPs, we predicted the genetic similarity for every pixel in the study area using the enriched adaptively RDA and environmental predictors (Capblancq & Forester, 2021;as described in Steane et al., 2014).
We used a k-means classification procedure as suggested in Considering the potentially small number of localities representing each class of climate-adapted genotypes, we complemented the geographical distribution of each class with verified occurrence records for the species retrieved from Sheu et al. (2020).We used the predicted values of genetic similarity for every pixel described above to classify the occurrence records into the respective climateadapted genotypes.Only records in which the predicted values of genetic similarity felt within the observed values for each genotype were included.We then filtered these additional records according to the geographical proximity to the observed climate-adapted genotypes (i.e.sampled individuals with genetic information), with increasingly larger buffers for each genotype up to the maximum extension in which the buffers do not overlap.
To model the distribution of the species and each climateadapted genotype, we used three different classes of algorithms for an ensemble forecasting (Araújo & New, 2007), including an envelope method (DOMAIN), logistic regression (GLM) and maximum entropy (Maxent) (Carpenter et al., 1993;Guisan et al., 2002;Phillips & Dudík, 2008).We performed 10 replications for each combination of variables for each method, using 10-fold cross-validation to evaluate the results.We kept the same proportion of randomly sampled pseudo-absences for each modelled class (2× the number of unique occurrences) to make the results comparable (Gomes et al., 2020).
We obtained the final ensemble models using true skill statistics (TSS) for weighting, where higher TSS values indicate more reliable model predictions, maximizing both model sensitivity and specificity (Allouche et al., 2006).Additionally to TSS, we also checked the area under the curve (AUC) for model performance, considering only models with AUC > 0.75 (Guisan et al., 2017).Finally, we projected the models into the current climate conditions and averaged the predictions for the five distinct GCMs for the SSP 3-7.0 (moderate/middle-range) and SSP 5-8.5 (extreme) emission scenarios We also noticed that commonly used thresholds such as maximum sensitivity + maximum specificity (max se + sp), although maximizing evaluation metrics such as max TSS and max Kappa (Liu et al., 2016), tended to crop out even areas with the presence of samples with genomic data.Finally, we cropped the resulting binary predictions for the present by the current forest cover.
We used dispersal simulations to predict whether each class of climate-adapted individuals would be able to track shifts in suitable climates and habitats until 2100, become extinct (low suitability or forest loss) or persist in each location.For this, we used a cellular automata model of dispersal based on kernel densities and barrier constraints, as implemented in the R package MigClim (Engler & Guisan, 2009).We ran MigClim until 2100 with intermediate steps in 2040 and 2070, using the ENMs predictions of habitat suitability for these years to indicate the potential to dispersal into unoccupied cells, where values lower than the minimum training presence were considered unsuitable (e.g. the site is considered unoccupied or extinct).We used the 2050 (for the above-mentioned intermediate steps in 2040 and 2070) and 2100 projections of land cover as strong barrier constraints, that is, dispersal was not enabled between two diagonally adjacent non-forest barrier cells.We aligned the land cover projection scenarios AB1 (moderate) and A2 (extreme deforestation) with corresponding moderate and extreme climate change scenarios.Complementary, to assess the relative influence of deforestation alone, we also investigated the effects of contrasting combinations of deforestation and climate scenarios (e.g.moderate deforestation with extreme climate change).Lastly, although no estimates of dispersal for K. calcarata exists, general estimates of range shifts due to climate change across multiple taxonomy groups range from a few metres to ~2 km per year (Chen et al., 2011).Therefore, we used maximum distance steps of 30 cells as the maximum distance possible to move between timesteps on each modelled interval (current to 2040, 2040-2070 and 2070-2100), which is equivalent to ~1 km per year, within the range estimated for other taxa (Chen et al., 2011) and at the same resolution of the land cover layers available.
The results from the ENM-based dispersal simulations were then used to indicate which parts of the ranges of distinct climateadapted genotypes will remain stable until 2070 and 2100 (climatically suitable and no deforestation), which new areas will become available and reachable through dispersal, and the range loss (climatically unsuitable or deforested).The potential for evolutionary rescue was estimated as the areas where the predicted distribution of one of the genotypes in 2100 overlaps with areas where the other genotype has lost its range (Figure 1).

| Results summary
We found considerable potential for evolutionary rescue in the moderate/middle-range scenarios of climate change, with most of the range of the species being rescued or remaining suitable until 2070.The genome-environment analyses revealed two genotypes uniquely adapted to different climates in Amazonia, after considering the effect of population structure.By integrating species distribution modelling with dispersal filters (forest cover), we were able to estimate range shifts (range loss in central and eastern Amazonia) and potential for evolutionary rescue in central Amazonia.

| Population structure
We found six geographically structured genetic clusters (k = 6) in K. calcarata (Figure 2 and Figure S3), with most individuals presenting some level of admixture, especially in eastern Amazonia (Figure 2a,b).Significant dispersal barriers roughly delimited genetic clusters, including river drainages, interfluves and the diagonal of dry biomes (Figure 2b,c).

| Genome-environment analyses
Our RDA results indicated a significant association of SNPs to the bioclimatic variables (adjusted r 2 = .046)when accounting for population structure (PCA), with the first three RDA axes explaining more than 80% of the variation (p < .01)(Figure S4).After correcting for the genomic inflation factor (GIF = 1.41), checking the uniform distribution of p-values and correcting for multiple tests (Figure S5), we found a set of 82 candidate SNPs with the RDA-only approach.The lowest value for the genomic inflation factor of our LFMM analyses for each of the 30,589 SNPs against three PCs axis of the environmental variables was reached for k = 6 (GIF = 1.36).After correcting for the genomic inflation factor (Figure S5) and for multiple tests, we found a set of 860 candidate SNPs with the LFMM-only approach.

| Climate-adapted individuals, modelling ranges and evolutionary potential
Of the SNPs found with RDA-only and LFMM-only approaches, We obtained good performances from the ENMs for each class of genotypes (dry/seasonal and wet/non-seasonal models) and for all individuals (species model), AUC > 0.9 and TSS > 0.6 (Figures S7-S9).The predicted range shift into future climates for all individuals (SSP 3-7.0 and SSP 5-8.5 scenarios), accounting for dispersal constraints and habitat availability, indicated large extensions of stable, suitable habitats until 2070 (Figure S10), but considerable losses in central Amazonia and south of the Amazon River until 2100 in the extreme scenario (Figure 3a; Figure S10).A similar pattern was found for dry/seasonal genotypes, with a high potential for range expansion only in the SSP 3-7.0 scenario (Figure 3c; Figure S11).Individuals with wet/non-seasonal genotypes will experience the highest degree of habitat loss and local extinction across the entire range, retaining just a few patches of suitable habitats until 2100 in either scenario (Figure 3b; Figure S12).When comparing the results considering dispersal constraints and habitat availability (Figure 3;  S7), where darker tones indicate a higher prevalence of each genomic cluster, and dots indicate all known species records used in the SDMs.The species presents a north-south population structuring separated by the Amazon River (number 1 in the map), a west to east differentiation roughly following some interfluves (2-Madeira and 3-Tapajós in green, 3-Tapajós and 4-Xingu in orange, 5-Tocantins/Araguaia in red, east Tocantins in blue) and an isolated cluster in the Atlantic Forest (purple).non-seasonal ones in the middle-range scenario (Figure 4a,b), with most areas being rescued or remaining suitable until 2070.For 2100, only a few patches of habitats were predicted to remain suitable but with considerable potential for evolutionary rescue.In the extreme scenario (Figure 4c,d), there is potential evolutionary rescue until 2070, but it is drastically reduced and confined to a small portion of the northern range in 2100, coincident with predicted remnants of forest cover (Figure S2e).No significant evolutionary rescue is predicted from the wet/non-seasonal genotypes to the dry-seasonal ones (not shown).When contrasting the effects of moderate and extreme scenarios of deforestation against both climate change scenarios, we consistently observed that moderate deforestation always results in greater areas with potential for evolutionary rescue, 76%-191.6%higher than extreme deforestation scenarios in 2070 and 2100 (Table S2).

| DISCUSS ION
We find a substantial potential for evolutionary rescue across natural populations of lizards in Amazonia, especially in moderate/middlerange emission scenarios until 2070, followed by a decreasing trend from 2070 to 2100 (Figure 4).We show that the persistence and dispersal of individuals with standing adaptive genetic adaptation may be an essential mechanism for biodiversity maintenance potentially preventing demographic processes that could lead to local extinction (Forester et al., 2022).This mechanism holds even without considering plasticity or the appearance of novel adaptive variation (Bay et al., 2017;Diniz-Filho & Bini, 2019).However, in extreme climate change scenarios, and if dispersal and persistence are restricted due to high pressures of deforestation, considerable parts of the species range will be lost together with local genetic variation across eastern, southern and central Amazonia.Considering that genetic lineages of K. calcarata (Figure 2), as well as populations of other organisms and even entire species are endemic to these areas (Ribeiro-Junior & Amaral, 2016), our results indicate a dire scenario for Amazonian biodiversity if extreme climate scenarios happen.
The distribution of the two climate-adapted genotypes identified for K. calcarata match important patterns of environmental conditions in the region.In Amazonia, a west-to-east climate gradient is marked by an increase in climate seasonality eastwards and in the direction of the Cerrado (Cheng et al., 2013).These gradients influence biodiversity patterns of multiple taxonomic groups at various levels, with biodiversity generally decreasing towards the east (Mesquita et al., 2015;Ter Steege et al., 2015).The distribution of K. calcarata mostly in eastern Amazonia may explain the dominance of the dry/ .Individuals with WnS genotypes will lose considerable extensions of suitable areas in both scenarios, whereas, DS and all individuals, considerable areas will remain stable mostly in the middle-range scenario.The extreme scenario indicates extensive habitat loss, especially to the south of the Amazon River, mainly related to the loss of forest cover (Figure S2e).Biome borders (grey lines) correspond to Am: Amazonia; Ce: Cerrado; Ca: Caatinga; At: Atlantic Forest.
seasonal genotype relative to the wet/non-seasonal one.Previous studies integrating mtDNA phylogeographic structure, thermal physiology and mechanistic distribution modelling of K. calcarata have shown a degree of thermal tolerance at the Amazonia-Cerrado ecotone, suggesting adaptation to intense climatic selective pressures (Avila-Pires et al., 2012;Cronemberger et al., 2022;Pontes-da-Silva et al., 2018), coincident with the distribution of the dominant dry/seasonal genotype.Although the current rate of environmental changes may require species to adapt or speciate faster than documented in evolutionary studies (Román-Palacios & Wiens, 2020), our results indicate that, at least for some species, the presence of standing genetic variation and local adaptation related to such climate gradient that could facilitate adaptation processes.
Beyond current climate gradients, eastern Amazonia and the ecotones have been climatically unstable since the Pleistocene (Cheng et al., 2013;Oliveras & Malhi, 2016) leading to habitat fragmentation and resulting in local extinctions, but also potentially providing opportunities for genetic divergence, speciation and adaptation in the long term (Baker et al., 2020).Changes in herbivory may have further strengthened the vegetation fluctuations, as until the Late Pleistocene, all of South America had a highly diverse megafauna (Faurby & Svenning, 2015).These may have substantially opened the vegetation (Doughty et al., 2016) which may have created warmer microhabitats than seen under the same climate today.Species distributed in these regions may have evolved adaptive responses, providing a repository of genetic diversity against future climate change across climatically unstable regions (Killeen & Solorzano, 2008).In that case, we could expect lower levels of extinction than predicted by our models due to climate change alone.However, deforestation is also predicted to be strong in this region, limiting evolutionary rescue, and some species do not seem to present adaptive genetic variation across F I G U R E 4 Potential for evolutionary rescue.Potential for evolutionary rescue from individuals with dry/seasonal (DS) genotypes to wet/ non-seasonal genotypes (WnS): the distribution of DS genotypes that will overlap with areas that will become unsuitable for WnS genotypes (compare to Figure 3), potentially allowing the permanence of the species in such areas.These results are shown in two climate change and forest cover scenarios, SSP 3-7.0/AB1 (moderate/middle-range) and SSP 5-8.5/A2 (extreme scenarios).The potential for evolutionary rescue is considerable until 2070, but sharply decreases afterwards.Biome borders (grey lines) correspond to Am: Amazonia; Ce: Cerrado; Ca: Caatinga; At: Atlantic Forest.
the same environmental gradient, as in the case of some anole lizards in Amazonia (Prates et al., 2018).
We show that climate change will mainly impact populations of K. calcarata currently adapted to the milder climates of centralsouthwestern Amazonia, contrary to the expectation of more intense anthropogenic climate change in eastern and southern Amazonia (Blois et al., 2013;Parsons, 2020).Most of the direction of change in suitable habitats for the wet/non-seasonal genotype was towards the southwest, where other species of the genus already occur (Sheu et al., 2020).Thus, tracking climate change in the case of this genotype could potentially increase interspecific competition or faunal turnover (Sales et al., 2020), although sympatric species of Kentropyx generally use distinct microhabitats (Avila-Pires, 1995).The direction of climate change in our models slightly favoured stability and expansion of habitats for populations with the dry-seasonal genotype in the direction of range loss of the wet/nonseasonal ones, thus increasing the potential for evolutionary rescue (Figures 2 and 3).However, this tendency was mostly restricted to the northernmost areas and a few additional small patches in the south, suggesting that the central-southern Amazonia might not have climates that are analogous to the current ones experienced by the species, a climate change tendency predicted to affect biodiversity hotspots profoundly (Williams et al., 2007).
Like previous studies using mtDNA data (Avila-Pires et al., 2012;Cronemberger et al., 2022), our genomic-level data indicate that K. calcarata presents geographically structured populations, associated with major interfluves in the transition zones between Amazonia and the Cerrado, as well as an isolated population in the Atlantic Forest (Figure 1).In addition, each main population occupies a distinct climate space (Figure S4a,b).Except for the one north of the Amazon River, all remaining populations are at great risk due to climate change per our future predictions, partially reflecting forecasts based on thermal physiology (Pontes-da-Silva et al., 2018).In the extreme deforestation scenario, all populations south of the Amazon River are predicted to lose most of their ranges (Figure 2).
The observed structured genetic clusters within the species suggest limited gene flow among populations, which could further diminish the potential for evolutionary rescue as indicated by our models.
Therefore, both climate-driven adaptive genetic variation and overall genetic diversity are in danger of being lost in large portions of southern and eastern Amazonia.
Modelling the species distribution into the future without considering the adaptive genetic variation and potential dispersal constraints resulted in slight differences in the predicted range shifts compared to the individual predictions for each main locally adapted genotype (Figure 3).We confirm a high total range loss tendency in modelling individual genotypes (Razgour et al., 2019).However, irrespective of the inclusion of adaptive genetic variation, we predicted substantial range loss in central-south Amazonia in moderate/ middle-range scenarios and across the entire southern Amazonian region in extreme scenarios.We also identified which part of the climate-related adaptive variation is at risk of disappearing by the end of the century (i.e.populations carrying genotypes adapted to less seasonal climates).Although vast areas were predicted to be unsuitable for the species, this does not necessarily translates into immediate extinction.Some areas may keep populations in scattered forest patches and small habitat patches (e.g.microrefugia), but they still might be affected by stochastic events, pathogens or inbreeding (extinction debt) (Kuussaari et al., 2009).
The potential for evolutionary rescue in buffering extinction risks for an ecologically important, forest-associated Amazonian lizard is considerably high until 2070.This potential remains high even in pessimistic climate change scenarios, but only as long as forest cover is retained so that it allows for spatial connectivity (at least 76.2% higher with moderate deforestation in 2100).However, in extreme deforestation scenarios-which could become a reality if Amazonian deforestation and other human impacts continues at the same levels seen over the past few years (Albert et al., 2023;Maeda et al., 2021;Silva Junior et al., 2021) -and at the emission levels predicted towards the end of the century, neither evolutionary rescue nor persistence of forest specialist species will be possible.Instead, range loss and genetic diversity erosion will occur, likely causing range-wide extinction of lizards and potentially other species.Our results suggest that if the goal is to reduce biodiversity loss in Amazonia, changes in land-use practices and actions to mitigate climate change should be considered.These actions could potentially provide sufficient time for demographic processes, adaptation and evolutionary rescue to occur, which are crucial for biodiversity recovery and maintenance.

ACK N OWLED G EM ENTS
quent analyses: Bio4 (temperature seasonality), Bio5 (max temperature of warmest month), Bio15 (precipitation seasonality) and Bio18 (precipitation of warmest quarter).First, we selected bioclimatic variables representing temperature and precipitation extremes and climate seasonality instead of annual averages.Such extreme temperatures are more strongly associated with local extinctions that have already occurred due to climate change worldwide (Román-Palacios & Wiens, 2020).We included only variables with considerable variation across the sampled localities, which after visual inspection of the environmental gradients resulted in 10 out of the 19 variables.We calculated the variation inflation factors (VIFs), with the usdm R package (Naimi & Araújo, 2016), keeping only variables with VIF lower than five to minimize multicollinearity among predictors.For mapping the current and future estimates of forest cover and land use in the study area, we used a high-resolution (1 km 2 ) global database for the years 2010, 2050 and 2100 derived from Li et al. (2017), which are based on remote sensing data, GCMs, future land use simulations (human development predictions) and , estimating Mahanolabis distances from the RDA scores (number of axes explaining at least 80% of the variation), which we then corrected for the inflation factor and transformed into p-values.Loci with p-values lower than .01divided by the number of SNPs (3.3 × 10 −7 -Bonferroni correction) were considered candidate SNPs.
Carvalho et al. (2021) on the adaptively enriched RDA scores (weighted by RDA axis importance) of each set of candidate SNPs to classify individuals as potentially adapted to distinct portions of the environmental space, hereafter climate-adapted genotypes.We used 30 indexes included in the R package NbClust to select the best number of k-means clusters (Charrad et al., 2014), choosing the value indicated by the majority of the indexes.From this classification, we used ENMs constrained by dispersal barriers (non-forested areas, see below) to predict current and future ranges of each climateadapted genotype and compared these results to the predicted range of the species using all records.
for 2011-2040 (hereafter 2040), 2041-2070 (hereafter 2070) and 2071-2100 (hereafter 2100).To calculate the current ranges (binary maps) of each climate-adapted class from the ENM outputs, we applied a threshold of the minimum training presence (Pearson et al., 2007) calculated considering only the records of individuals with genomic data confirmed for each genotype class, allowing to retrieve the minimal conditions for the presence of each genotype.
only 25 were common to both analyses and used in downstream analyses.The distribution of the individuals across the adaptively enriched multivariate space (i.e.RDA scores containing only candidate SNPs) resulted in two main groups (k-means classification): one from central-southwestern Amazonia, which we classified as a genotype climatically adapted to wet/non-seasonal climates, and another one from eastern Amazonia and the ecotones with the Cerrado and northern savannas, classified as the dry-seasonal genotype (Figure S6).Alternative selections of candidate SNPs (RDA only, LFMM only and RDA + LFMM union) yielded quite similar results (Figure S6).

Figures
Figures S10-S12) with the suitability values of the raw SDM results (Figures S7-S9), it is noticeable that climate change alone is responsible for a great extent of areas in central Amazonia and south of the Amazon River being considered areas of local extinction for the species in 2100.We found considerable potential for evolutionary rescue from dry-seasonal genotypes to the areas dominated by the wet/