# Future climate will drive changes of suitable habitats and sympatric areas for two green lizards in Western Europe

Igor Boyer, Romain Bertrand, Olivier Lourdais, Francis Isselin-Nondedeu

► To cite this version:

Igor Boyer, Romain Bertrand, Olivier Lourdais, Francis Isselin-Nondedeu. Future climate will drive changes of suitable habitats and sympatric areas for two green lizards in Western Europe. Frontiers of Biogeography, 2023, 15 (4), pp.e59543. 10.21425/f5fbg59543. hal-04318021

HAL Id: hal-04318021

https://hal.science/hal-04318021v1

Submitted on 1 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Future climate will drive changes of suitable habitats and sympatric areas for two green lizards in Western Europe

Igor Boyer<fnref n="1" />,<fnref n="2" />* ID, Romain Bertrand<fnref n="3" /> ID, Olivier Lourdais<fnref n="4" />,<fnref n="5" /> ID, and Francis Isselin-Nondedeu<fnref n="6" />,<fnref n="7" /> ID

**Abstract**

Species may respond to climate change by redistributing their distribution areas, but because they do not share the same climatic affinities, they should not respond in the same way. Consequently, distribution shifts of species that are currently found in sympatry may change the extent of the area of sympatry and therefore interspecific interactions at the local scale. In Western Europe, the green lizards *Lacerta agilis* and *Lacerta bilineata* live in partial sympatry, share morphological similarities, and can locally compete for resources. In this study, we used a correlative species distribution model (SDM), Maxent, to explore the effects of climate change on the distribution of suitable areas for each species and also within their sympatric area under future scenarios. Our simulations showed that all *L. agilis* subspecies are more likely to lose suitable habitats throughout their distribution areas whatever the scenario. Conversely, even if *L. bilineata* should lose less suitable areas, gains may be more restricted. In addition, when looking within the areas both reachable for *L. a. agilis* and *L. bilineata* to identify where sympatric areas will occur, we find that they should be spatially restricted. Climatic refugia potentially shared for both lizards may therefore be limited. Consequently, competitive exclusion expected in areas suitable for both species (i.e., the potential area of sympatry) could decrease in the future as climate changes likely lead *L. agilis* subspecies and *L. bilineata* to follow different species range shifts trajectories over the 21st century.

**Highlights**

* Closely related species living in partial sympatry should respond differently to future climate change

* Green lizards are not expected to gain suitable areas in their leading edges and therefore are unlikely to colonize newly suitable habitats in their northern range margins

* Climate change should prevail over competition in the reduction of suitable areas within areas that are reachable for both species

* Adding scenarios of competitive interaction to the responses to climate change strengthens predictions of species distribution

* Our results highlight the transverse effect of climate changes for biodiversity, promoting its redistribution and mismatch on species interactions.

**Keywords:** Lacertids, Maxent, range shift, species competition, species distribution models, species lineages

<footnote n="1">¹Laboratoire Cltés, TERritoires, Environnement et Sociétés, UMR 7324, Université de Tours, CNRS, Tours, France;</footnote>

<footnote n="2">²Laboratoire Biodiversité dans l'Anthropocène : Dynamique, Fonction et Gestion (BiodivAG), Université d'Angers, Angers, France;</footnote>

<footnote n="3">³Laboratoire Évolution et Diversité Biologique, UMR 5174, Université de Toulouse III Paul Sabatier, CNRS, IRD, Toulouse, France;</footnote>

<footnote n="4">⁴Centre d'Etudes Biologiques de Chizé, CEBC UMR 7372, CNRS-La Rochelle Université, 79360, Villiers en Bois, France;</footnote>

<footnote n="5">⁵School of Life Sciences, Arizona State University, Tempe, AZ, 85287-4501, USA;</footnote>

<footnote n="6">⁶Institut Méditerranéen de Biodiversité et Écologie (IMBE), UMR Université Aix-Marseille Avignon, 7263-CNRS, 237-IRD IRPNC (Ingénierie de la Restauration des Patrimoines Naturels et Culturels), Avignon, France;</footnote>

<footnote n="7">⁷Département d'Aménagement et d'Environnement, Ecole Polytechnique de l'Université de Tours, Tours, France;</footnote>

<footnote>*Correspondence: Igor Boyer, pro.igorboyer@gmail.com</footnote>
Introduction

Terrestrial species are facing marked spatial
redistribution along latitudinal or elevational gradients
in response to climate change (Lenoir et al. 2020).
Species' abilities to cope with these changes depend
on the adequacy between their response capacities
(dispersal, plasticity and adaptation) and the pace of
climate change (Kafash et al. 2018, Loarie et al. 2009).
Because species differ in their ecological niches their
future range-shifts may result either in expansion or
contraction of their distribution areas. The species
redistribution may also reshape the location and
extent of current sympatric and syntopic zones
(i.e., areas where two [or more] related species exist
together and live in the same habitats) (Urban et al.
2012). Despite the large body of scientific literature
that stressed future changes in suitable habitats,
few have considered the responses of co-existing
species (but see Boukal et al. 2019, Jones et al. 2016).
However, since related species may differ in their
dispersal abilities, sympatric or syntopic species may
not be able to track future environmental changes
in their habitats (in terms of pace and direction).
Consequently, some species could lose habitats,
while others could gain habitats within the sympatric
zones (Préau et al. 2020, Urban et al. 2012). This
eventually leads to the disappearance of sympatric
zones (Supplementary Fig. S1) and alters community
compositions. On the other hand, species range-shifts
could induce no-analog communities with newly
sympatric species that have similar ecological niches
and consequently will likely increase interspecific
competition (Bombi et al. 2011, Urban et al. 2012).

Among ectotherms, reptiles have received little attention in conservation-prioritization studies even if it has been shown that they are experiencing both decline and distribution changes (Cox et al. 2022, Diele-Viegas & Rocha 2018). Indeed, Moreno-Rueda et al. (2012) evidenced a current northward distribution shift for 22 of 30 species of reptiles in the Iberian Peninsula, in response to the current climate change. Reptiles that spread poleward or upward access to newly suitable habitats and could be serious competitors for resident species. For example, in England, range-expanding populations increased by 40% for *Podarcis muralis* and 36% for *Lacerta bilineata*, whereas the population of common lizards (*Zootoca vivipara*), which are resident, decreased by approximately 75% between 2002 and 2007 (Mole 2010). Competition may not be the only factor that explains these changes in population trends, indeed climate change is likely to favor non-native species (Mole 2010).

Reptiles are closely related to the abiotic condition
of their environment. Temperature and humidity
strongly influence their physiological functions, such
as locomotion, growth, and reproduction (Rozen-
Rechels et al. 2019) which possibly impacts their
fitness (Dezetter et al. 2021, Reedy et al. 2013). Both
thermal preferences and limits are specific to each
species and relate to their life histories, especially
their foraging modes (Corn 1971, Trochet et al. 2018).

However, contrasted thermal tolerances should result
in different responses to climate change (Catry et al.
2015, Díaz-Ricaurte et al. 2020). Species that coexist
in sympatry and use similar habitats may nonetheless
differ in their thermal and hydric requirements, partly
because of their life history (Guillon et al. 2014,
Lelèvre et al. 2010, Lourdais et al. 2013). These
mechanistic differences between sympatric species
are likely to influence the response to the ongoing
climate change. Many studies have already investigated
the impact of climate change on reptile distributions,
however, studies on the consequences of species'
redistributions on the sympatric zones (geographic
position, extent) are still lacking.

Locally, syntopic lizards may compete for both trophic resources and basking or retreat sites that allow them to maintain an optimal body temperature (Bergeron & Blouin-Demers 2020, Chukwuka et al. 2021). Competition for basking sites occurs between both conspecifics and other lizards (Langkilde and Shine 2004). Competitive outcomes generally go in the favor of larger species, regardless of their age classes (Langkilde and Shine 2004, 2007). Although aggressive interactions can occur between species, new aggressive interactions can also occur when species expand into the range of a resident species that has a similar Eltonian niche because these separate species are competing for resources (Grether et al. 2013, Shine and Langkilde 2005).

Correlative SDMs are one of the most used tools to
assess the vulnerability of species to climate change,
however the use of abiotic variables alone (such as
climate variables) cannot fully explain the distribution
of species which also rely on biotic interactions. Indeed,
at the local scale, species with similar ecological
requirements and partially overlapping climate
niche may compete for resources in their habitats
(Jones et al. 2016). Moreover, areas that are highly
climatically suitable for two or more species that have
similar ecological requirements are supposed to be
those where the outcomes of competition will be the
higher (because species will show high performance
and are therefore good competitors). Thus, we can
suppose that areas predicted highly suitable for two
species that have similar ecological requirements
would not become as suitable as predicted because
one species may outcompete the other one, resulting
in a reduced number of predicted suitable areas
(Jones et al. 2016).

In this study, we focus on two species of Western
European Lacertidae lizards belonging to the genus
*Lacerta*: the sand lizard *Lacerta agilis* (Linnaeus 1758) and
the western green lizard *Lacerta bilineata* (Daudin 1802).
*L. agilis* is observed mainly in open habitats with
compact bushes and is distributed over the temperate
Palearctic (Heltai et al. 2015, Nemes et al. 2006,
Roitberg et al. 2015), whereas *L. bilineata* is considered
both as a sub-Mediterranean species with thermophile
preferences and a generalist species that inhabits open
areas (Pikalik et al. 2020). They differ substantially
in their morphology and life-history (Arnold 1998,
Mezzasalma et al. 2014) but share a partial sympatric
distribution, mainly in the center of France (Fig. 1).
Escoriza and Amat (2021), who focused on niche segregation in lizards, found that *L. bilineata* has a wider niche than *L. agilis* and that they have great probability to be found in the same habitats. *L. agilis* numbers are declining in Western Europe, and an action plan has been developed to take appropriate measures to conserve this species (Edgar and Bird 2005).

Given their inherent differences we expect these species to respond differently to climate change and these species appear to be good candidates for testing whether climate change should redistribute suitable habitats in Europe and reshape their sympatric areas (Supplementary Fig. S1). Therefore, we asked the following questions:

1. Are the bioclimatic niches of *L. agilis* (sensu lato) and *L. bilineata* divergent enough to induce contrasting responses to climate change? Given their current global distributions, we expect that the bioclimatic niches of *L. agilis* and *L. bilineata* strongly differ. Moreover, in response to climate warming and accounting for potential dispersal abilities, we expect that *L. bilineata* will gain more habitats than *L. agilis*.

2. How will climate change reshape the current area of sympatry and competitive interactions affect the extent of potential presence of the two lizards? We predict that climatically suitable areas for both species will decrease, and that competitive exclusion within the remaining sympatric areas will exacerbate these decreases.

To address these questions, we used the species distribution modelling (SDM) software MaxEnt (Phillips et al. 2017) to predict future changes in habitat suitability for the periods 2041–2060 and 2061–2080. We consider the species' dispersal abilities by integrating reachable habitats (i.e., habitats that could theoretically be reached by individual species according to their dispersal abilities) into the models. Because *L. agilis* is divided into several distinct subspecies, we first tested for niche similarity between subspecies.

## Materials & Methods

### Study species

Sand lizards are divided into subspecies that are grouped into two major clades: the western clade, which comprises the *L. a. agilis*, *L. a. argus*, *L. a. chersonensis*, and *L. a. garzoni*, and the eastern clade, which comprises the *L. a. brevicaudata*, *L. a. exigua*, *L. a. grusinica*, and *L. a. tauridica* (Andres et al., 2014). *L. bilineata* belongs to the *Lacerta viridis* complex gathering European green lizards that belong to four deeply divergent lineages, such as *L. bilineata*, *L. viridis*, the Adriatic lineage of *L. viridis*, and a *L. viridis* lineage confined to the Black Sea coast of Anatolia (Marzahn et al. 2016). With isolated populations in Germany (Fig. 1), *L. bilineata* is distributed in the north of Spain, Italy, and France. This species is divided into two subspecies, with *L. b. chloronota* being endemic to Sicily and Calabria and *L. b. bilineata* on the rest of the distribution area (Marzahn et al. 2016).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Current distribution areas and locations of the sympatric areas of *Lacerta bilineata* and *Lacerta agilis* subspecies in Western Europe (source: IUCN, GBIF, INPN [Inventaire National du Patrimoine Naturel]). Maps were projected under the Lambert Azimuthal Equal Area (LAEA) Europe coordinate system (EPSG:3035). The colonization limits (dotted lines) represent the geographical areas that the species could potentially reach between 2041 and 2060 based on their respective dispersal capacities only.</figcaption></figure>
## Occurrence data

Occurrence data for *Lacerta agilis* and *Lacerta bilineata* were obtained from GBIF.org (accessed September 2021) GBIF Occurrence Download (DOI: https://doi.org/10.15468/dl.etdvus for *L. bilineata* and https://doi.org/10.15468/dl.znucb8 for *L. agilis*). We selected all occurrences from the Western European distribution areas of *Lacerta agilis agilis*, *Lacerta agilis argus*, *Lacerta agilis garzoni*, and *Lacerta bilineata* that have a maximal coordinate uncertainty of 5.000 meters for the period 1979-2016. Occurrence data were projected to the Lambert Azimuthal Equal Area (LAEA) Europe coordinate system (EPSG:3035) and firstly scaled to correspond to the centroid of a 10 x 10-km grid. This data uniformization allowed us to reduce the relative part of oversampled locations and, thus, increase the model's predictive abilities by reducing spatial autocorrelation (Boria et al. 2014).

Secondly, we removed occurrence data located outside the species' current ranges obtained through the International Union for Conservation of Nature Red List (UICN redlist) website (*L. agilis*: https://dx.doi.org/10.2305/IUCN.UK-4.RLTS.T15788A5071439.en accessed on 2021-04; *L. bilineata*: https://dx.doi.org/10.2305/IUCN.UK.2009.RLTS.T61519A12501065.en accessed on 2020-02) (Fig. 1), the sample size of occurrence data used for modeling was 3,356 for *L. bilineata*, 1,673 for *L. a. agilis*, 371 for *L. a. argus*, and 21 for *L. a. garzoni*. This last subspecies only inhabits the eastern part of the Pyrenean mountains, and according to the IUCN, our occurrence data covered most of this subspecies' current distribution area. Of the 3,356 occurrences of *L. bilineata*, only seven could be attributed to *L. b. chloronota* based on location. Of the 3,356 and 1,673 occurrences of *L. bilineata* and *L. a. agilis*, respectively, 423 were sympatric (i.e., the same coordinates).

## Climate data for species distribution modeling

Initially, we selected 11 bioclimatic variables that influence lizards' behavioral thermohydroregulation and, therefore, possibly their distribution (Rozen-Rechels et al. 2019, Žagar et al. 2015). We used isothermality (Bio3), maximum temperature of the warmest month (Bio5), minimum temperature of the coldest month (Bio6), mean temperature of the driest quarter (Bio9), precipitation seasonality (Bio15), precipitation of the warmest quarter (Bio18), precipitation of the coldest quarter (Bio19), radiation of the first trimester (RadTri1), radiation of the second trimester (RadTri2), radiation of the third trimester (RadTri3), and radiation of the fourth trimester (RadTri4). Isothermality, as a proxy of thermal fluctuation through seasons, is an important driver of reptile persistence (Taylor et al. 2021). In summer, maximum temperature could raise reptiles' body temperatures above their maximum critical temperatures, leading to their deaths. Conversely, minimal temperature in winter may lower body temperature below their minimum critical temperature but also greatly impact their body condition and ultimately their persistence (Brischoux et al. 2016).

In addition, performance of reptiles is closely related to the hydration state of individuals, climatic variables related to precipitation are therefore greatly important for reptile persistence (Rozen-Rechel et al. 2019). Also, changes in precipitation regime during the winter may facilitate the expansion of reptiles in new habitats (Araújo et al., 2006). Finally, heliothermic ectotherms, such as reptiles, highly depend on solar radiation to thermoregulate (Bogert 1949, Žagar et al. 2015). Quarterly radiation data were used to account for seasonality changes.

We calculated the current bioclimatic variables using the R package dismo (Hijmans et al. 2021) for the period 1981–2010 at a 10-km grid spatial resolution based on monthly precipitation, temperature, and radiation data obtained from the TerraClimate database (Abatzoglou et al. 2018). We conducted all analysis with R version 4.2.1 (R Core Team 2021). This resolution was selected to match with our occurrence data. We calculated future climate conditions for the periods 2041–2060 and 2061–2080 under the two representative concentration pathway (RCP) scenarios RCP2.6 and RCP8.5 (optimistic and pessimistic scenarios, respectively) using Tabor and Williams's (2010) delta method. Thus, we calculated climate anomalies between the current and future periods from the full dataset of CMIP5 GCMs used in the IPCC WG1 AR5 (van Oldenborgh et al. 2013). These anomalies were calculated as absolute anomalies for temperature and relative anomalies for precipitations and radiations (Trouet and Van Oldenborgh 2013, van Oldenborgh et al. 2013). Then, we resampled future anomaly grids through bilinear interpolation to obtain 10-km grids covering most of Europe. We finally added the 10-km grids to a second set of high-resolution climate grids covering the temperature, precipitation, and radiation conditions averaged during the period 1981–2010 (Abatzoglou et al. 2018).

We calculated the Pearson correlation coefficients of the bioclimatic variables, as well as their variance inflation factors (VIFs) in the SDM, using the R packages usdm and ecospat (Broennimann et al. 2012, Naimi et al. 2014). We retained only variables with a VIF of less than 10 and an absolute correlation coefficient of less than 0.75 (Naimi et al. 2014). We employed this method to limit collinearity issues in the model, which could have biased the SDM's forecasts. Predictor variables must be selected according to the biology of the studied species to better capture the species' ecological niche (Araújo et al. 2019). We retained six variables for this study: the maximum temperature of the warmest month (Bio5), minimum temperature of the coldest month (Bio6), precipitation seasonality (Bio15), precipitation of the warmest quarter (Bio18), precipitation of the coldest quarter (Bio19), and the radiations of the first trimester (RadTri1). Hence, despite their Pearson correlation coefficients that were slightly higher than 0.75, RadTri1 and Bio5 were retained because their VIFs indicated inconsequential collinearity and they are important for temperate reptile species (O'Brien 2007; Table S1, Table S2).
## Dispersal capacities and reachable habitats

To account for species' dispersal capacities, we confined our projections for future climatic suitability onto dispersal buffers. These buffers represent the geographical space that the species could potentially reach based on their maximum dispersal capacities. We identified the reachable habitats by multiplying the maximum yearly dispersal distance by the number of years between the current period and the future period. We considered the dispersal capacities based on a study by Berglind (2000), who monitored a population of *L. a. agilis* in Sweden. In Berglind's study, a subadult who moved approximately 500 m was responsible for the highest recorded dispersal event. Based on this result, we considered a dispersal distance of 500 m per year. Because no study has assessed reliable dispersal distance for *L. bilineata*, we used the dispersal distance of its sister species, *Lacerta viridis* (Kolora et al. 2018, Tvrtkovi et al. 1998). *L. viridis* can disperse over approximately 5,000 m per year. The buffer widths were therefore calculated as 5,000 x 41 for *L. bilineata* (i.e., 205 km) and 500 x 41 for *L. agilis* (*sensu lato*; i.e., 20.5 km) for the period 2041–2060.

## Species distribution modeling

We used MaxEnt v3.4.1 (using a "cloglog" transformation) to model the distribution of each lizard species based on the relationship between species' occurrences and bioclimatic conditions. Maxent is commonly used for presence-only data and has also demonstrated satisfactory performance in situations with limited occurrence data (Elith et al. 2006, Hernandez et al. 2006, Phillips et al. 2017, Tan et al. 2022).

Although *L. bilineata* is divided into subspecies, occurrence data within the distribution area of *L. b. chloronota* were too scarce (seven) to warrant individual consideration. Therefore, these occurrence data were considered as part of the 3,356 occurrences attributed to the nominal species *Lacerta bilineata*. Our occurrence dataset gathered occurrences of both *L. b. bilineata* and *L. b. chloronota* (the only two subspecies of *L. bilineata*), which we hereinafter categorize under *L. bilineata*. Conversely, based on the results that the niche similarity tests provided (see Appendix S1), we ran separate models for the subspecies of *L. agilis* (i.e., *L. a. agilis*, *L. a. argus*, and *L. a. garzoni*), which allowed us to partially account for intraspecific diversity that could lead to differences in species' range shifts (Boyer et al. 2021). Also, even close species with apparently similar ecological requirements, should respond differently to climate change because their current distributions may be likely driven by specific climatic conditions (Labadessa and Ancillotto 2022, Struck et al. 2018).

We calibrated models with current bioclimatic data within a buffer of 200 km to have an extended calibration area with enough cells to generate background points without duplicates. In addition, background points are more likely to be true absences when they are distant from known presence points (Grimmett et al. 2020). Also, this calibration area is not oversized, limiting the inclusion of climatically suitable habitats that are not reachable, due to the presence of dispersal barrier.

We randomly selected 75% of data for training (random seed and bootstrap) and retained the remaining 25% for testing. MaxEnt's bootstrap resampling method allowed for the selection of a random testing dataset for each run. We carried out the model performance evaluation with the testing dataset. Although testing datasets can be partly composed of occurrences also used within the training dataset, this partitioning method performs satisfactorily, especially with species with few occurrence data (Hassangavyar et al. 2022).

Furthermore, to avoid overcomplexity and overfitting issues in our models, we tested a range of regularization multipliers and feature classes (Blair et al. 2022, Radosavljevic and Anderson 2014) using the R package ENMeval (Muscarella et al. 2014). We tested all possible feature classes combinations between the Linear, Hinge, Quadratic, Product and Threshold features with regularization multipliers ranging from 1 to 3. We selected the most adapted combinations according to the lowest value of AICc.

The best regularization multipliers and feature classes were 1LQHPT (hinge plus quadratic plus hinge plus product plus threshold features) for *L. a. agilis*, *L. a. argus* and *L. a. garzoni* and 1LQ (linear plus quadratic features) for *L. bilineata*. We set a prevalence of 0.1 and randomly generated 3346, 3339, 189 and 6712 background points within the calibration area for *L. a. agilis*, *L. a. argus*, *L. a. garzoni* and *L. bilineata* respectively. We set a total of 100 runs for model building and set the other parameters by default for each lizard.

We calculated the area under the receiver operating characteristic curve (AUC), which is a metric of the predictive performance of a model that ranges from 0 to 1, where 0.5 and below indicate that the model does not have a better predictive ability than a random model (Hanley and McNeil 1982). Because AUC received criticism, we also calculated the continuous Boyce index, Kappa statistic and true skill statistic (TSS). Contrary to the AUC, Kappa statistics and TSS are threshold-dependent measures of accuracy (Allouche et al. 2006). Also, CBI assesses a model's ability to consistently predict several levels of suitability and is particularly relevant for presence-only data (Hirzel et al. 2006). These metrics range from -1 to 1. For the CBI and TSS negative values indicate an incorrect model, while values close to 1 indicate that predictions perfectly infer suitability indexes from occurrence data. Kappa metric ranges from -1 to 1 where $K < 0.40$ indicates poor reliability, $0.40 < K > 0.8$ good reliability and $K > 0.80$ almost perfect reliability (Landis and Koch 1977).

We projected future suitability, within the previously described dispersal buffers, with the future climatic data calculated following the delta method at the 10-km square resolution. We binarized the continuous predicted habitat suitability into suitable and unsuitable data using Maxent's "maximum training sensitivity plus specificity" threshold rule (Liu et al. 2013). This threshold maximizes the fraction of correctly predicted occurrences (true presence) and the fraction of correctly predicted absences (true absence; Barbet-Massin et al. 2012).
We therefore calculated the percentage of losses and gains in climatically suitable habitats. We hereafter combined changes in suitability for each lineage on averaged maps.

To identify the direction in which the global distribution of lizards will shift in the coming decades, we assessed the location of the centroids for current and future suitable areas and calculated the differences in latitudinal and longitudinal coordinates between these centroids.

We assessed the uncertainties in our projections using two methods. Multivariate environmental similarity surface (MESS) analyses made it possible to identify areas where climatic conditions were outside the range of those used for calibration, thereby inducing extrapolations in predictions. MESS analyses return values ranging from -100 to 100 where "values < -10 indicate high extrapolation, where values of -10 to 0 indicate low extrapolation, values of 0 to 10 low interpolation, and where values > 10 indicate high interpolation" (Ngo et al. 2022). We also calculated the standard deviation across forecasts to distinguish areas where the predictions converged (the more reliable predictions) from those where predictions were divergent and, therefore, less reliable.

Habitat redistributions within the potential sympatric zone in competitive interactions

We investigated the possible shifts in, as well as the contractions, expansions, and proportions of, suitable habitats within the future sympatric zones in three steps. We defined the future sympatric zones as the geographical areas that will be both reachable (based on their dispersal capacities) and abiotically (here climatically) suitable for both species.

We created buffers taking into account the intrinsic dispersal abilities of *L. a. agilis* and *L. bilineata* in order to precisely identify reachable areas for 2041-2060 and 2061-2080. We used precise occurrences obtained from the "Faune-France" network of associations for the period 1981-2010 as a starting point for the creation of buffers. We then intersected the two buffers (i.e., *L. a. agilis* buffer and *L. bilineata* buffer) to identify areas reachable for both species.

Second, we identified the climate suitability within these reachable areas. We used the model outputs previously obtained to forecast species' climate suitability within their distribution areas, at the resolution of 10 x 10 km. We cropped these forecasts with the buffer of reachable areas for both species to focus on the areas of potential sympatry.

Third, we investigated the consequences of competitive outcomes within the potential future sympatric areas for the total surface of habitats for *L. a. agilis* and *L. bilineata*. These species coexist at the coarse scale but rarely at the local scale. Indeed, when looking for accurate occurrences within the 423 sympatric occurrences at the coarse scale (10-km square), we find 167 precise occurrences of *L. agilis*, 516 for *L. bilineata* and only 40 syntopic occurrences (the precise occurrences come from high resolution dataset described above).

Therefore, competitive exclusion may occur between these species even if it was not proven (no study assessed competitive interaction between these species). Because no study assessed competitive exclusion between these species, we do not know in which competition may occur. However competitive exclusion should occur between two similar species when habitat is suitable for both (Real et al. 2022, Wang et al. 2005). We therefore explored three a priori theoretical conditions within cells, which were predicted as both climatically suitable and reachable for the two species (i.e., potential sympatric areas):

Case 1: Neither of the two species excludes the other when they coexist (coexistence scenario).

Case 2: *L. agilis* outcompetes and excludes *L. bilineata* (competitive exclusion).

Case 3: *L. bilineata* outcompetes and excludes *L. agilis* (competitive exclusion).

Initially, we tested the coexistence scenario (Case 1) wherein the areas that are reachable and climatically suitable for the species remain suitable for both species. To investigate Cases 2 and 3, we verified where areas will be abiotically suitable for both *L. bilineata* and *L. a. agilis* within reachable areas. When a 10-km square cell was predicted to be climatically suitable for both species, the cell was attributed to *L. agilis* in Case 2 and *L. bilineata* in Case 3.

## Results

### Current climate suitability

According to mean AUC<sub>test</sub>, CBI and Kappa statistic, the bootstrapped maxent models' performance is acceptable, and predicted suitability values are reliable (AUC<sub>mean</sub> = 0.826 ± 0.008; 0.892 ± 0.017; 0.951 ± 0.034; 0.806 ± 0.005 | CBI<sub>mean</sub> = 0.996 ± 0.003; 0.985 ± 0.01; 0.874 ± 0.057; 0.998 ± 0.002 | Kappa<sub>mean</sub> = 0.489 ± 0.017; 0.423 ± 0.049; 0.608 ± 0.122; 0.45 ± 0.009 for *L. a. agilis*, *L. a. argus*, *L. a. garzoni* and *L. bilineata* respectively [see Table S3 for details relative to other metrics]). According to model forecasts (Supplementary Fig. S2), an important part of the distribution area (defined using IUCN data) is climatically suitable for each lizard.

The niche overlap analysis does not show that niches between lizards are more similar than random, they should therefore respond differently to changes in their environment (Appendix S1). According to the niche analysis, the bioclimatic variables that contributed the most for western green lizards were Bio6 for *L. a. agilis* and *L.a. argus*, Bio5 for *L. a. garzoni* and RadTri1 for *L.bilineata* (Appendix S1).

For *L. a. agilis*, areas located in the range margins appeared less suitable. Most areas within the distribution area of *L. bilineata* appeared suitable, and only some areas located in its southern range margins and in altitude (e.g., in the Alps). Areas that are climatically unsuitable for *L. a. agilis* and *L. bilineata* covered about 40 to 50% of their current distribution areas (Supplementary Fig. S2).
Changes in suitability within distribution areas

MESS analyses indicated that most reachable
areas had low degrees of extrapolation for most
lizards, regardless of scenario (Supplementary Fig. S3).
Consequently, we did not predict suitability values from
environmental values outside those used in calibration,
and our predictions were reliable. However, suitability
in elevated areas tended to be more extrapolated
(Supplementary Fig. S3).

The standard deviation across each run followed a
similar pattern with more variation in elevated areas
(Supplementary Fig. S4). We noted that areas with
high standard deviations were often predicted to be
unsuitable (Fig. 2).

According to the location of the distribution
centroids, we noted that suitable areas for the complex
of *L. agilis* subspecies were always more displaced to
those for *L. bilineata* whatever the period or scenario.
Also, suitable areas for *L. agilis* shifted to the northeast
(Fig. 3). With the notable exception of models projected
in 2061-2080 under the pessimistic scenario, suitable
areas for *L. bilineata* shifted to the north.

Under the RCP2.6, forecasts indicated less losses of
currently suitable habitats for western green lizards.
Most of the time, gains in newly suitable areas are
greater under this scenario but never offset losses,
thereby resulting in a net decrease in the surface
of suitable areas of approximately -9.97% to -0.56%
(Table 1).

Losses of currently suitable habitats for *L. agilis* and *L. bilineata* were more important under RCP8.5 than under RCP2.6, ranging from -86.96% to -27.64% and -27.90% to -8.51% for RCP8.5 and RCP2.6, respectively, regardless of period (Table 1). Gains in newly suitable habitats were always more important under RCP8.5 than under RCP2.6 for *L. agilis* but not for *L. bilineata* (respectively +27.39% to +45.36% vs +3.03% to 5.89%). However, gains never compensated for losses of suitable habitats, which resulted in a global net reduction in the surface of suitable area of approximately -71.29% to -0.56% compared to current (Table 1).

Under RCP2.6, the losses of suitable habitats for *L. agilis* in the western range margins were partially offset by gains in the eastern and northern range margins (Fig. 4). Therefore, the net surface area of suitable future habitats ranged from -9.79% to -2.98% of currently suitable areas (Table 1). In addition, findings relating to RCP8.5 indicated that losses of currently suitable habitats were more than one and a half greater than the surface area of newly suitable habitats. The remaining suitable areas were 59.57% fewer than currently suitable areas in the period 2061-2080 (Table 1). Consequently, suitable habitats had almost disappeared in western Europe (Fig. 4). Nevertheless, areas in the northern and eastern parts of the distribution area remained and became climatically suitable, regardless of period or scenario, and populations inhabiting Denmark, Sweden, or the eastern part of Poland persisted over time (Fig. 4). Additionally, we noted that the Pyrenean mountains (i.e., where the subspecies *Lacerta agilis garzoni* currently lives) remained suitable, regardless of scenario.

<figure><img src="image_4.png" /></figure>

**Figure 2.** Current climatic suitabilities (a) of each *Lacerta agilis* subspecies in Western Europe within their distribution areas (performed separately; AGI = *Lacerta agilis agilis*; ARG = *Lacerta agilis argus*; and GAR = *Lacerta agilis garzoni*. The dotted line separates the distribution areas of *L. a. agilis* from those of *L. a. argus* and (b) of *Lacerta bilineata* within its distribution areas. We projected maps under the Lambert azimuthal equal-area coordinate system (EPSG:3035) at a 10 x 10-km spatial resolution.

<figure><img src="image_5.png" /></figure>

**Figure 3.** Mean shifts of distribution centroids for the complex of *Lacerta agilis* subspecies and *Lacerta bilineata* (see Materials and Methods).
<table><caption>**Table 1.** Gains in and losses of suitable habitats compared to numbers of currently suitable areas for green lizards in Western Europe. The column "Suitable areas loss (%)" represents the percentage of currently suitable habitats that were lost, and the column "Suitable areas gain (%)" represents the percentage of suitable habitats that were gained, with both sets of percentages being based on the period and the RCPs. For the complex of the *Lacerta agilis* subspecies, we trained separate models that were averaged and binarized according to the "maximum training sensitivity plus specificity" threshold rule for each subspecies. Subsequently, we pooled binarized predictions on one map per period and RCP.</caption><thead><tr><th>Species</th><th>Period</th><th>RCP</th><th>Suitable areas loss (%)</th><th>Suitable areas gain (%)</th><th>Gains/Losses balance</th></tr></thead><tbody><tr><td rowspan="4">L. agilis</td><td rowspan="2">2041-2060</td><td>2.6</td><td>27.9</td><td>17.92</td><td>-9.97</td></tr><tr><td>8.5</td><td>61.6</td><td>45.36</td><td>-16.24</td></tr><tr><td rowspan="2">2061-2080</td><td>2.6</td><td>24.14</td><td>21.16</td><td>-2.98</td></tr><tr><td>8.5</td><td>86.96</td><td>27.39</td><td>-59.57</td></tr><tr><td rowspan="4">L. bilineata</td><td rowspan="2">2041-2060</td><td>2.6</td><td>10.83</td><td>8.08</td><td>-2.74</td></tr><tr><td>8.5</td><td>27.64</td><td>5.89</td><td>-21.75</td></tr><tr><td rowspan="2">2061-2080</td><td>2.6</td><td>8.51</td><td>7.95</td><td>-0.56</td></tr><tr><td>8.5</td><td>74.32</td><td>3.03</td><td>-71.29</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Changes in future climatic suitability for the complex of *Lacerta agilis* subspecies in Western Europe. We ran separate models for each subspecies and pooled them on one map per period and RCP (see Methods). Red indicates future unsuitable areas (i.e., currently suitable areas); dark blue areas remain suitable; light blue areas become suitable (i.e., currently unsuitable areas); and gray areas remain unsuitable. We projected maps at a 10 x 10-km spatial resolution.</figcaption></figure>
Under RCP2.6, we noted that losses of suitable habitats for *L. bilineata* were partially counterbalanced by gains, which resulted in a net surface area ranging from -2.74% to -0.56% of suitable habitats compared to numbers of currently suitable areas (Table 1). Models did not forecast a northward expansion in response to climate change. More than 75% of the newly reachable areas were not predicted to become suitable in the coming decades; therefore, habitats reached in the northern range margin could be unsuitable (Fig. 5). Consequently, the location of suitable habitats throughout the distribution area of *L. bilineata* remained in its current state (Fig. 5).

Losses of currently suitable habitats under RCP8.5, occur throughout France for the period 2061-2080, represented -74.32% (Table 1; Fig. 4). Lastly, contrary to those for *L. a. agilis*, under RCP8.5, newly suitable habitats for *L. bilineata* were fewer in number than they were under RCP2.6 (i.e., +5.89% versus +8.08% and +3.03% versus +7.95% for the periods 2041-2060 and 2061-2080, respectively; Table 1).

*Suitability changes in competitive exclusion scenarios*

Within areas that are reachable for both species, areas suitable for both species (i.e., *L. agilis* and *L. bilineata*), which may represent the future area of syntopy, covered 2.65% to 21.94% of the total

surface (see Table S4). Among these areas, those that are located in mid-elevation areas or in altitude (e.g., in the Massif Central and the Alps) seems to persist through time even if others decreased (Fig. 6). Moreover, under the pessimist scenario, only about 10306 km² of areas remain suitable for both species left in 2061-2080. Conversely, unsuitable areas for both species increase through time and cover most of the reachable areas for both species (Table S4).

Under RCP2.6, habitats that are climatically suitable only for *L. agilis* exceed in number those for *L. bilineata* in the period 2041–2060 (Fig. 7). According to our forecast, under RCP2.6, suitable habitats for *L. agilis* were mainly located in the northern half of the accessible zone, unlike those for *L. bilineata*, which were always located in the southern part of zone, regardless of period or scenario (Fig. 6). Under RCP8.5, suitable habitats for *L. agilis* decreased over time (from 54,900 km² in the period 2041–2060 to 29,800 km² in the period 2061–2080) and only survived in the eastern and northern parts of the accessible zone (Fig. 6).

Case 1: Coexistence scenario. We observed that, under RCP2.6, the surface of suitable habitats was larger for *L. agilis* than they were for *L. bilineata*, regardless of period (e.g., representing about 48.1% to 70.4% of the reachable zone for both species versus 39.7% to 52.5% for *L. bilineata* in the period

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Changes in future habitat suitability for *Lacerta bilineata* (see details in Fig. 3)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Future suitability of reachable areas for both species (according to a starting point corresponding to current occurrences). We obtained maps by binarizing and averaging all runs. Purple and green indicate areas reachable for both species and only predicted suitable for *L. bilineata* and *L. agilis*, respectively. Yellow indicates areas that are reachable and suitable for both species (i.e., future sympatric zones), and gray indicates areas unsuitable but reachable.</figcaption></figure>

2041–2060; Table S5). Nevertheless, we noted that, under RCP8.5, suitable habitats were numerous for *L. bilineata*, regardless of period (Fig. 7).

Case 2: *L. agilis* outcompetes *L. bilineata*. Under RCP2.6, the median surface of suitable areas for *L. bilineata* may be approximatively halved in the period 2041–2060 compared to the “no competition” condition (Fig. 7; Table S5). Additionally, we noted that, under RCP8.5, numbers of suitable habitats for *L. bilineata* were comparable to those for *L. a. agilis*. For the period 2061–2080, we did not note consequent changes in the percentage of suitable habitats for *L. bilineata*.

Case 3: *L. bilineata* outcompetes *L. agilis*. Suitable habitats were always greater in number for *L. bilineata* than for *L. agilis*, regardless of period or scenario (Table S5). The most important reduction in the median

surface of suitable areas for *L. agilis* was observed under the RCP2.6 and represent a decrease about 36% whatever the period. Under RCP8.5, suitable habitats for *L. agilis* had almost disappeared in the period 2061–2080 and only represented approximately 3.5% to 18.8% of the total surface area of the potential accessible zone.

## Discussion

In the present study we find that western green lizards may be unequally threatened by climate change that agrees with the findings of studies that have shown that genetically close species should respond differently to climate change (Labadessa & Ancillotto 2022, Struck et al. 2018). We forecast that *L. bilineata* is likely to keep most of its suitable habitats over the 21st
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Percentage of suitable habitats regarding reachable areas for both species depending on the three a priori theoretical conditions for *Lacerta agilis* and *Lacerta bilineata*. We obtained values through 100 runs for each species. Cases 2 and 3 attributed habitats suitable for both species only to the dominant species (see Methods).</figcaption></figure>

century, whereas the subspecies complex of *L. agilis* is expected to suffer from the effects of global warming and face a severe loss of suitable habitats. However even if *L. bilineata* appears less threatened compared to *L. agilis* subspecies complex, it is unlikely to colonize newly suitable areas in its range margin.

Within the accessible area for both species, areas predicted to be suitable for both *L. a. agilis* and *L. bilineata*, wherein competitive interactions may occur, will cover a restricted surface area compared to those suitable for only one species. If competitive exclusion occurs, the surface of climatically suitable areas for the two species will not change overall.

In this study, we projected the future distribution changes for *Lacerta agilis* subspecies and *Lacerta bilineata* by considering dispersal limitations to determine whether one of these lizards could expand further than the other, and we explored the consequences for their sympatric zones. Moreover, our framework is reliable because it accounts for practices recommended in Blair et al. (2022), addressing model complexity, biases in input data,

and uncertainties in model projections. We obtained the maximum dispersal distances used in this study through studies that assessed maximum dispersal distances in homogenous landscapes. However, there is no guarantee that these distances are as high in fragmented landscapes. Consequently, our dispersal buffers had to overestimate the number of newly suitable areas in leading edges. We favored overestimating reachable areas over considering no dispersal or unlimited dispersal, which could have led to unrealistic predictions (Araújo et al. 2006, Ramirez-Villegas et al. 2014).

## Contrasted spatial responses

Our forecasts suggest that *Lacerta agilis* (*sensu lato*) will experience an unprecedented reduction in suitable areas throughout its current distribution area. Indeed, the distribution centroid shifts exceeds the dispersal abilities of this species, indicating a considerable contraction of its distribution areas to the northeast. Also, these discrepancies indicate a future lag between the ability of this species to cope with suitability
changes and its dispersal ability. In comparison, shifts of the distribution centroids for *L. bilineata* are less pronounced and do not exceed its dispersal abilities. Nevertheless, we noted this species may suffer from an important range contraction toward elevated areas in 2061-2080 under the pessimistic scenario.

Losses of suitable areas will be exacerbated under RCP8.5, leading to the spatial isolation of suitable habitats in the period 2061–2080 for both species. Suitable areas are mostly predicted to be found in mountainous areas (i.e., in Alps and Massif Central). This pattern has been observed in other reptiles that, despite climate change, could keep suitable habitats high in the mountains (Moreno-Rueda et al. 2012, Ngo et al. 2022). Nevertheless, successful migration and the establishment of populations at higher altitudes involve physiological adaptations relating to oxygen capacity and metabolism to counterbalance altitude hypoxia (Taylor et al. 2021). Moreover, although our forecasts indicate that large geographical areas are currently suitable regarding climate, the loss of heathland is known to pose a severe threat to sand lizards, whose populations are currently declining considerably and who are becoming increasingly isolated (Edgar and Bird 2005).

Our forecasts indicate future range contractions for *L. agilis*, although we did not include landscape variables because of methodological constraints (e.g., the irrelevance of the resolution of climatic variables for landscape). In addition, high resolution landscape data that account for heathland at the Europe scale are scarce, not uniformized among countries, obtained through distinct research programs, and therefore hardly comparable. Consequently, the range contraction predicted for *L. agilis* could be even more important if its habitats were to disappear, which could lead to the disappearance of most of the current populations. However, we observed that the Pyrenean mountains, where *L. a. garzoni* lives, will remain climatically suitable. This subspecies may therefore persist through time contrary to lowland subspecies.

We predicted that *L. bilineata* will lose fewer habitats than *L. agilis* (*sensu lato*), and, as our forecasts indicate, climate change should pose a lesser threat to *L. bilineata*. Our model identified that *Bio6* was an important variable for two of the three *L. agilis* subspecies. Contrary to *L. bilineata* that may accommodate from higher temperatures than *L. agilis* which seem infeduated from coldest temperatures according to response curves. Hence, a hypothetical global increase in temperature should not reduce the number of suitable future habitats for *L. bilineata*, which is consistent with the findings of a study that found that Mediterranean species were projected to lose fewer suitable habitats in Southern Europe than Euro-Siberian species were (Ruiz-Labourdette et al. 2012). Nevertheless, thermal variables seem not to be the main drivers of *L. bilineata* presence. Conversely, according to our forecasts, reachable habitats in the northern range margin are unlikely to become suitable soon, and therefore, higher dispersal capacities do not systematically result in more newly suitable habitats. Lastly, according to its responses to climate change, *L. bilineata* could, to

a certain extent, be considered the potential future winner, while *L. agilis* (*sensu lato*) could be the loser. Although we may have overestimated the dispersal distances used in this study, especially for *L. bilineata*, we did not observe that dispersal capacities could mitigate the impact of climate change. More than 75% of the areas potentially reachable through dispersal for this species will not become suitable, regardless of scenario (Fig. 5). Within its northern range margin *L. agilis* may locally benefit from climate change through phenology changes bringing a higher survival and fitness (Ljungström et al. 2015). Nonetheless, within its trailing edges or even its core distribution, this species may not take advantage of climate change.

## Effect of climate change on the sympatric zone

SDMs at the 10-km square resolution indicate that the distribution area of *L. bilineata* may not change much contrary to *L. a. agilis* which should suffer from a severe decrease in suitability. This finding follows the third case of our conceptual scheme (Supplementary Fig. S1), where Species B does not lose habitats in its range, while Species A does lose habitats, and the sympatric zone is reduced or even lost.

However, higher-resolution projections indicate a low number of suitable habitats for both species. We initially predicted that *L. bilineata* should gain newly suitable habitats where *L. agilis* currently lives because the species inhabits warmer areas. However, our results do not support this hypothesis, as climatic suitability is not expected to increase at all in the northern part of the accessible areas for both species over time and *L. bilineata* is not expected to colonize new areas in its leading edge. Also, these species appear strongly divergent in their climatic affinities, resulting in few areas that are suitable for both.

Because of these strong decreases in suitability, species must redistribute them. However, landscape connectivity may not be high enough to allow these species to redistribute themselves in higher altitudes (Guo et al. 2018). Consequently, because suitable areas will be reduced, those suitable for both species would be even more limited. Only one third of suitable habitats for *L. agilis* will also be suitable for *L. bilineata* in the period 2061–2080 within the reachable area for both species. Consequently, when we considered the case where *L. bilineata* outcompetes *L. agilis*, we found that the latter would almost disappear from Western Europe. We must keep in mind that even though these lost areas induced the extinction of *L. agilis* in its southern range margin, these areas only cover a small surface and are not of major importance compared to areas that remain suitable elsewhere in its distribution.

Another problem arises like our forecasts indicate that areas in altitude will remain climatically suitable for *L. agilis*, regardless of scenario. However, most of these areas are disconnected from each other, and sand lizard populations in elevated areas are likely to suffer from severe population isolation, which is one of the main drivers of species extinction (Heinrichs et al., 2016). Additionally, all these areas will be suitable for *L. bilineata*.
Therefore, nothing can guarantee that *L. bilineata* will not exclude sand lizards from these habitats.

Nevertheless, of the studies that have focused on niche segregation or competition between *L. agilis* and other lizard species, none have concluded that *L. agilis* could suffer from interspecific competition or exclude potential competitors from its habitats. However, Hamilton (2021) observed aggressive interactions between *L. agilis* and *P. muralis*, contrary to what had been observed in Germany (Heym et al. 2013). *L. agilis* may therefore exhibit similar behaviors to those of other potential competitors, such as *L. bilineata*. Langilde and Shine (2007), meanwhile, observed that competitive outcomes generally go in favor of the larger species, which, in our case, is *L. bilineata*. Moreover, this species is known to occasionally feed on reptiles, and consequently, the species is likely to exclude *L. agilis* from habitats suitable for both species. Lastly, we forecast that the potential future sympatric areas will be limited and regardless of the outcomes of competition, suitable areas for *L. agilis* decrease through time. Although *L. agilis* may be systematically excluded from habitats where *L. bilineata* lives, few habitats would be lost through exclusion.

## Conclusion

We found that climate change will negatively impact *L. agilis*. Because *L. bilineata* is projected to lose few of its suitable habitats, climate change could threaten this species less than the other one. However, newly suitable future habitats are unlikely to be plentiful for this species. Climatically suitable habitat for both species that could represent the future sympatric zones will often be located in altitude, decrease through time and only cover a small surface within areas reachable for both species. Considering the biology and associated competitive abilities of *L. bilineata*, the outcomes of competitive interactions may be an additional driver of range contraction for *L. agilis*, even if areas where competition should occur are spatially restricted. Although biotic interactions are often disregarded, their impact on the distribution of species is undeniable, and suitable areas could therefore be more limited than expected.

## Acknowledgements

We would like to thank Philipp Ginal and the two anonymous reviewers for their relevant comments which improved the manuscript and Darren McIntyre for English language corrections. The funding was supported by the Region Centre-Val de Loire (France) under the acronym MODELISE. RB's work was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41).

## Author Contributions

I.B., R.B., O.L. and F.I.N. conceived the ideas; I.B. performed the analysis, R.B., co-analysed some of the data; I.B. wrote, R.B., O.L. and F.I.N reviewed the paper. F.I.N. managed the whole project.

## Data Accessibility Statement

Occurrence data used in the present study were obtained from GBIF.org (accessed September 2021) GBIF Occurrence Download (DOI: https://doi.org/10.15468/dl.etdvus for *L. bilineata* and https://doi.org/10.15468/dl.znubc8 for *L. agilis*).

## Supplemental Material

The following materials are available as part of the online article at https://escholarship.org/uc/fb

**Figure S1.** Conceptual scheme of different expected range shifts scenarios for sympatric species in response to climate change.

**Table S1.** Pearson correlation matrix of the bioclimatic variables.

**Table S2.** Variation inflation factor scores of the climatic variables selected throughout modeling processes.

**Table S3.** Mean AUC and Boyce Index scores of the forecasts made for each western green lizard.

**Figure S2.** Percentage areas that are currently suitable within the distribution areas of each western green lizard.

**Figure S3.** Multivariate environmental similarity surface (MESS) analysis of future climate suitability.

**Figure S4.** Coefficient of variation across predictions.

**Table S4.** Surface and percentage of suitable areas within areas that are both reachable for *L. bilineata* and *L. agilis*.

**Table S5.** Minimal, maximal and median surface of suitable areas (in %) that are both reachable for *L. bilineata* and *L. agilis* across period and scenario.

**Appendix S1:** Niche overlap analysis

## References

Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., & Hegewisch, K. C. (2018) TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific Data, 5, 170191. https://doi.org/10.1038/sdata.2017.191

Allouche, O., Tsoar, A., & Kadmon, R. (2006) Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology, 43, Article 6. https://doi.org/10.1111/j.1365-2664.2006.01214.x

Andres, C., Franke, F., Bleidorn, C., Bernhard, D., & Schlegel, M. (2014) Phylogenetic analysis of the *Lacerta agilis* subspecies complex. Systematics and Biodiversity, 12, 43–54. https://doi.org/10.1080/14772000.2013.878000

Araújo, M. B., Thuiller, W., & Pearson, R. G. (2006) Climate warming and the decline of amphibians and reptiles in Europe. Journal of Biogeography, 33, Article 10. https://doi.org/10.1111/j.1365-2699.2006.01482.x
Araújo, M. B., Anderson, R. P., Márcia Barbosa, A., et al. (2019) Standards for distribution models in biodiversity assessments. Science Advances, 5, eaat4858. https://doi.org/10.1126/sciadvaat4858

Arnold, E. N. (1998) Cranial kinesis in lizards: variations, uses, and origins. Evolutionary Biology, 323–357.

Barbet-Massin, M., Jiguet, F., Albert, C. H., & Thuiller, W. (2012) Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution, 3, 327–338. https://doi.org/10.1111/j.2041-210X.2011.00172.x

Bergeron, R., & Blouin-Demers, G. (2020) Niche partitioning between two sympatric lizards in the Chiricahua Mountains of Arizona. Copeia, 108. https://doi.org/10.1643/CH-19-268

Berglind, S.-A. (2000). Demography and management of relict sand lizard *Lacerta agilis* populations on the edge of extinction. Ecological Bulletins, 48, 123-142.

Blair, M. E., Le, M. D., & Xu, M. (2022) Species distribution modeling to inform transboundary species conservation and management under climate change: promise and pitfalls. Frontiers of Biogeography, 14, e54662. https://doi.org/10.21425/F5FBG54662

Bogert, C. M. (1949) Thermoregulation in reptiles, a factor in evolution. Evolution, 3, 195–211.

Bombi, P., Akani, G. C., Ebere, N., & Luiselli, L. (2011) Potential effects of climate change on high- and low- abundance populations of the Gaboon viper (Bitis gabonica) and the nose-horned viper (B. nasicornis) in southern Nigeria. The Herpetological Journal, 21, 59–64.

Boria, R. A., Olson, L. E., Goodman, S. M., & Anderson, R. P. (2014) Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecological Modelling, 275, 73–77. https://doi.org/10.1016/j.ecolmodel.2013.12.012

Boukal, D. S., Bideault, A., Carreira, B. M., & Sentis, A. (2019) Species interactions under climate change: connecting kinetic effects of temperature on individuals to community dynamics. Current Opinion in Insect Science, 35, 88–95. https://doi.org/10.1016/j.cois.2019.06.014

Boyer, I., Cayuela, H., Bertrand, R., & Isselin-Nondedeu, F. (2021) Improving biological relevance of model projections in response to climate change by considering dispersal amongst lineages in an amphibian. Journal

of Biogeography, 48, 561–576. https://doi.org/10.1111/jbi.14019

Brischoux, F., Dupoué, A., Lourdais, O., & Angelier, F. (2016) Effects of mild wintering conditions on body mass and corticosterone levels in a temperate reptile, the aspic viper (Vipera aspis). Comparative Biochemistry and Physiology A, 192, 52–56. https://doi.org/10.1016/j.cbpa.2015.11.015

Broennimann, O., Fitzpatrick, M. C., Pearman, P. B., et al. (2012) Measuring ecological niche overlap from occurrence and spatial environmental data: measuring niche overlap. Global Ecology and Biogeography, 21, Article 4. https://doi.org/10.1111/j.1466-8238.2011.00698.x

Catry, I., Catry, T., Patto, P., Franco, A., & Moreira, F. (2015) Differential heat tolerance in nestlings suggests sympatric species may face different climate change risks. Climate Research, 66, 13–24. https://doi.org/10.3354/cr01329

Chukwuka, C. O., Mello, R. S. R., Cree, A., & Monks, J. M. (2021) Thermal heterogeneity of selected retreats in cool-temperate viviparous lizards suggests a potential benefit of future climate warming. Journal of Thermal Biology, 97, 102869. https://doi.org/10.1016/j.jtherbio.2021.102869

Corn, M. J. (1971) Upper thermal limits and thermal preferenda for three sympatric species of Anolis. Journal of Herpetology, 5, 17. https://doi.org/10.2307/1562838

Cox, N., Young, B. E., Bowles, P., et al. (2022) A global reptile assessment highlights shared conservation needs of tetrapods. Nature, 605, 285–290. https://doi.org/10.1038/s41586-022-04664-7

Dezetter, M., Le Galliard, J. F., Guiller, G., Guillon, M., Leroux-Coyau, M., Meylan, S., Brischoux, F., Angelier, F., & Lourdais, O. (2021) Water deprivation compromises maternal physiology and reproductive success in a cold and wet adapted snake *Vipera berus*. Conservation Physiology, 9, coab071. https://doi.org/10.1093/conphys/coab071

Díaz-Ricaurte, J. C., Serrano, F. C., Guevara-Molina, E. C., Araujo, C., & Martins, M. (2020) Does behavioral thermal tolerance predict distribution pattern and habitat use in two sympatric Neotropical frogs? PLoS ONE, 15,
e0239485. https://doi.org/10.1371/journal.pone.0239485

Diele-Viegas, L. M., & Rocha, C. F. D. (2018) Unraveling the influences of climate change in Lepidosauria (Reptilia). Journal of Thermal Biology, 78, 401–414. https://doi.org/10.1016/j.jtherbio.2018.11.005

Edgar, P., & Bird, D. R. (2005) Action plan for the conservation of the Sand Lizard (Lacerta agilis) in Northwest Europe. Convention on the Conservation of European Wildlife and Natural Habitats, Strousbourg.

Elith, J., H. Graham, C., P. Anderson, R., et al. (2006) Novel methods improve prediction of species' distributions from occurrence data. Ecography, 29, Article 2. https://doi.org/10.1111/j.2006.0906-7590.04596.x

Escoriza, D., & Amat, F. (2021) Habitat partitioning and overlap by large lacertid lizards in Southern Europe. Diversity, 13, 155. https://doi.org/10.3390/d13040155

Grether, G. F., Anderson, C. N., Drury, J. P., Kirschel, A. N. G., Losin, N., Okamoto, K., & Peiman, K. S. (2013) The evolutionary consequences of interspecific aggression: aggression between species. Annals of the New York Academy of Sciences, 1289, 48–68. https://doi.org/10.1111/nyas.12082

Grimmett, L., Whitsed, R., & Horta, A. (2020) Presence-only species distribution models are sensitive to sample prevalence: evaluating models using spatial prediction stability and accuracy metrics. Ecological Modelling, 431, 109194. https://doi.org/10.1016/j.ecolmodel.2020.109194

Guillon, M., Guiller, G., DeNardo, D. F., & Lourdais, O. (2014) Microclimate preferences correlate with contrasted evaporative water loss in parapatric vipers at their contact zone. Canadian Journal of Zoology, 92, 81–86. https://doi.org/10.1139/cjz-2013-0189

Guo, F., Lenoir, J., & Bonebrake, T. C. (2018) Land-use change interacts with climate to determine elevational species redistribution. Nature Communications, 9, 1315. https://doi.org/10.1038/s41467-018-03786-9

Hamilton, R. (2021) Aggressive interactions between a sand lizard *Lacerta agilis* and a common wall lizard *Podarcis muralis* in England. Herpetological Bulletin, 155, 28–29. https://doi.org/10.33256/hb155.2829

Hanley, J. A., & McNeil, B. J. (1982) The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology, 143, Article 1. https://doi.org/10.1148/radiology.143.1.7063747

Hassangavyar, M. B., Damaneh, H. E., Pham, Q. B., Linh, N. T. T., Tiefenbacher, J., & Bach, Q.-V. (2022) Evaluation of re-sampling methods on performance of machine learning models to predict landslide susceptibility. Geocarto International, 37, 2772–2794. https://doi.org/10.1080/10106049.2020.1837257

Heinrichs, J. A., Bender, D. J., & Schumaker, N. H. (2016) Habitat degradation and loss as key drivers of regional population extinction. Ecological Modelling, 335, 64–73. https://doi.org/10.1016/j.ecolmodel.2016.05.009

Heltai, B., Sály, P., Kovács, D., & Kiss, I. (2015) Niche segregation of sand lizard (Lacerta agilis) and green lizard (Lacerta viridis) in an urban semi-natural habitat. Amphibia-Reptilia, 36, 389–399. https://doi.org/10.1163/15685381-00003018

Hernandez, P. A., Graham, C. H., Master, L. L., & Albert, D. L. (2006) The effect of sample size and species characteristics on performance of different species distribution modeling methods. Ecography, 29, 773–785. https://doi.org/10.1111/j.0906-7590.2006.04700.x

Heym, A., Deichsel, G., Hochkirch, A., Veith, M., & Schulte, U. (2013) Do introduced wall lizards (Podarcis muralis) cause niche shifts in a native sand lizard (Lacerta agilis) population? A case study from south-western Germany. Salamandra, 9, 97–104.

Hijmans, R. J., Phillips, S., Leathwick, J., & Elith, J. (2021) dismo: species distribution modeling. Available at https://CRAN.R-project.org/package=dismo

Hirzel, A. H., Le Lay, G., Helfer, V., Randin, C., & Guisan, A. (2006) Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling, 199, 142–152. https://doi.org/10.1016/j.ecolmodel.2006.05.017

Jones, M., Bertola, L., & Razgour, O. (2016) Predicting the effect of interspecific competition on habitat suitability for the endangered African wild dog under future climate and land cover changes.
Hystrix, the Italian Journal of Mammalogy, 27.
https://doi.org/10.4404/hystrix-27.1-11678

Kafash, A., Ashrafi, S., Ohler, A., Yousefi, M., Malakoutikhah, S., Koehler, G., & Schmidt, B. R. (2018) Climate change produces winners and losers: differential responses of amphibians in mountain forests of the Near East. Global Ecology and Conservation, 16, e00471. https://doi.org/10.1016/j.gecco.2018.e00471

Kolora, S. R. R., Weigert, A., Saffari, A., et al. (2018) Divergent evolution in the genomes of closely-related lacertids, *Lacerta viridis* and *L. bilineata* and implications for speciation. *GigaScience*. https://doi.org/10.1093/gigascience/giy160

Labadessa, R., & Ancillotto, L. (2022) A tale of two crickets: global climate and local competition shape the distribution of European Oecanthus species (Orthoptera, Gryllidae). *Frontiers of Biogeography*, 14, e57538. https://doi.org/10.21425/F5FBG57538

Landis, J. R., & Koch, G. G. (1977) The measurement of observer agreement for categorical data. *Biometrics*, 33, 159. https://doi.org/10.2307/2529310

Langkilde, T., & Shine, R. (2004) Competing for crevices: interspecific conflict influences retreat-site selection in montane lizards. *Oecologia*, 140, 684–691. https://doi.org/10.1007/s00442-004-1640-1

Langkilde, T., & Shine, R. (2007) Interspecific conflict in lizards: social dominance depends upon an individual's species not its body size. *Austral Ecology*, 32, 869–877. https://doi.org/10.1111/j.1442-9993.2007.01771.x

Lelièvre, H., Blouin-Demers, G., Bonnet, X., & Lourdais, O. (2010) Thermal benefits of artificial shelters in snakes: a radiotelemetric study of two sympatric colubrids. *Journal of Thermal Biology*, 35, 324–331. https://doi.org/10.1016/j.jtherbio.2010.06.011

Lenoir, J., Bertrand, R., Comte, L., Bourgeaud, L., Hattab, T., Murienne, J., & Grenouillet, G. (2020) Species better track climate warming in the oceans than on land. *Nature Ecology and Evolution*, 4, 1044–1059. https://doi.org/10.1038/s41559-020-1198-2

Liu, C., White, M., & Newell, G. (2013) Selecting thresholds for the prediction of species occurrence with presence-only data. *Journal of Biogeography*, 40, Article 4. https://doi.org/10.1111/jbi.12058

Ljungström, G., Wapstra, E., & Olsson, M. (2015) Sand lizard (Lacerta agilis) phenology in a warming world. *BMC Evolutionary Biology*, 15, 206. https://doi.org/10.1186/s12862-015-0476-0

Loarie, S. R., Duffy, P. B., Hamilton, H., Asner, G. P., Field, C. B., & Ackerly, D. D. (2009) The velocity of climate change. *Nature*, 462, Article 7276. https://doi.org/10.1038/nature08649

Lourdais, O., Guillon, M., DeNardo, D., & Blouin-Demers, G. (2013) Cold climate specialization: adaptive covariation between metabolic rate and thermoregulation in pregnant vipers. *Physiology and Behavior*, 119, 149–155. https://doi.org/10.1016/j.physbeh.2013.05.041

Marzahn, E., Mayer, W., Joger, U., et al. (2016) Phylogeography of the *Lacerta viridis* complex: mitochondrial and nuclear markers provide taxonomic insights. *Journal of Zoological Systematics and Evolutionary Research*, 54, 85–105. https://doi.org/10.1111/jzs.12115

Mezzasalma, M., Maio, N., & Guarino, F. M. (2014) To move or not to move: cranial joints in European Gekkotans and Lacertids, an osteological and histological perspective: Gekkotan and Lacertid cranial kinesis. The *Anatomical Record*, 297, 463–472. https://doi.org/10.1002/ar.22827

Mole, S. R. C. (2010) Changes in relative abundance of the western green lizard *Lacerta bilineata* and the common wall lizard Podarcis muralis introduced onto Boscombe Cliffs, Dorset, UK. *The Herpetological Journal*, 15, 6.

Moreno-Rueda, G., Pleguezuelos, J. M., Pizarro, M., & Montori, A. (2012) Northward shifts of the distributions of Spanish reptiles in association with climate change. *Conservation Biology*, 26, 278–283. https://doi.org/10.1111/j.1523-1739.2011.01793.x

Muscarella, R., Galante, P. J., Soley-Guardia, M., Boria, R. A., Kass, J. M., Uriarte, M., & Anderson, R. P. (2014) ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. *Methods in Ecology and Evolution*, 5, Article 11. https://doi.org/10.1111/2041-210X.12261

Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K., & Toxopeus, A. G. (2014) Where is positional uncertainty a problem for species distribution
modelling? Ecography, 37, Article 2. https://doi.org/10.1111/j.1600-0587.2013.00205.x

Nemes, S., Vogrin, M., Hartel, T., & Öllerer, K. (2006) Habitat selection at the sand lizard (Lacerta agilis) ontogenetic shifts. North-Western Journal of Zoology, 2, 17-26.

Ngo, H. N., Nguyen, H. Q., Phan, T. Q., Nguyen, T. Q., Gewiss, L. R., Rödder, D., & Ziegler, T. (2022) Modeling the environmental refugia of the endangered Lichtenfelder's Tiger Gecko (Goniurosaurus lichtenfelderi) towards implementation of transboundary conservation. Frontiers of Biogeography, 14, e51167. https://doi.org/10.21425/F5FBG51167

O'Brien, R. M. (2007) A caution regarding rules of thumb for variance inflation factors. Quality and Quantity, 41, 673-690. https://doi.org/10.1007/s11135-006-9018-6

Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E., & Blair, M. E. (2017) Opening the black box: an open-source release of Maxent. Ecography, 40, Article 7. https://doi.org/10.1111/ecog.03049

Pikalik, M., Scali, S., Sacchi, R., Majlathova, V., & Majlath, I. (2020) I love my ecotone: habitat selection by the Western green lizard (Lacerta bilineata). Herpetology Notes, 13, 1027-1033.

Préau, C., Grandjean, F., Sellier, Y., Gailledrat, M., Bertrand, R., & Isselin-Nondedeu, F. (2020) Habitat patches for newts in the face of climate change: local scale assessment combining niche modelling and graph theory. Scientific Reports, 10, 3570. https://doi.org/10.1038/s41598-020-60479-4

Radosavljevic, A., & Anderson, R.P. (2014) Making better Maxent models of species distributions: Complexity, overfitting and evaluation. Journal of Biogeography, 41, 629-643. https://doi.org/10.1111/jbi.12227

Ramirez-Villegas, J., Cuesta, F., Devenish, C., Peralvo, M., Jarvis, A., & Arnillas, C. A. (2014) Using species distributions models for designing conservation strategies of Tropical Andean biodiversity under climate change. Journal for Nature Conservation, 22, 391-404. https://doi.org/10.1016/j.jnc.2014.03.007

Real, R., Báez, J. C., Fa, J. E., Olivero, J., & Acevedo, P. (2022) Making the competitive exclusion principle operational at the biogeographical

scale using fuzzy logic. Frontiers in Ecology and Evolution, 10, 991344. https://doi.org/10.3389/fevo.2022.991344

Reedy, A. M., Zaragoza, D., & Warner, D. A. (2013) Maternally chosen nest sites positively affect multiple components of offspring fitness in a lizard. Behavioral Ecology, 24, 39-46. https://doi.org/10.1093/beheco/ars133

Roitberg, E. S., Eplanova, G. V., Kotenko, T. I., Amat, F., Carretero, M. A., Kuranova, V. N., Bulakhova, N. A., Zinenko, O. I., & Yakovlev, V. A. (2015) Geographic variation of life-history traits in the sand lizard, *Lacerta agilis*: testing Darwin's fecundity-advantage hypothesis. Journal of Evolutionary Biology, 28, 613-629. https://doi.org/10.1111/jeb.12594

Rozen-Rechels, D., Dupoué, A., Lourdais, O., Chamaillé-Jammes, S., Meylan, S., Clobert, J., & Le Galliard, J. (2019) When water interacts with temperature: Ecological and evolutionary implications of thermo-hydroregulation in terrestrial ectotherms. Ecology and Evolution, 9, 10029-10043. https://doi.org/10.1002/ece3.5440

Ruiz-Labourdette, D., Nogués-Bravo, D., Ollero, H. S., Schmitz, M. F., & Pineda, F. D. (2012) Forest composition in Mediterranean mountains is projected to shift along the entire elevational gradient under climate change: Forest dynamics under climate change. Journal of Biogeography, 39, 162-176. https://doi.org/10.1111/j.1365-2699.2011.02592.x

Shine, R., & Langkilde, T. (2005) How do water skinks avoid shelters already occupied by other lizards? Behaviour, 142, 203-216. https://doi.org/10.1163/1568539053627659

Struck, T. H., Feder, J. L., Bendiksby, M., et al. (2018) Finding evolutionary processes hidden in cryptic species. Trends in Ecology and Evolution, 33, 153-163. https://doi.org/10.1016/j.tree.2017.11.007

Tan, W. C., Ginal, P., Rhodin, A. G. J., Iverson, J. B., & Rödder, D. (2022) A present and future assessment of the effectiveness of existing reserves in preserving three critically endangered freshwater turtles in Southeast Asia and South Asia. Frontiers of Biogeography, 14, e50928. https://doi.org/10.21425/F5FBG50928

Taylor, E. N., Diele-Viegas, L. M., Gangloff, E. J., et al. (2021) The thermal ecology and physiology of reptiles and amphibians: A user's guide.
Journal of Experimental Zoology A: Ecological and Integrative Physiology, 335, 13-44. https://doi.org/10.1002/jez.2396

Trochet, A., Dupoué, A., Souchet, J., et al. (2018) Variation of preferred body temperatures along an altitudinal gradient: a multi-species study. Journal of Thermal Biology, 77, 38-44. https://doi.org/10.1016/j.jtherbio.2018.08.002

Trouet, V., & Van Oldenborgh, G. J. (2013) KNMI climate explorer: a web-based research tool for high-resolution paleoclimatology. Tree-Ring Research, 69, 3-13. https://doi.org/10.3959/1536-1098-69.1.3

Tvrtkovi, N., Lazar, B., Tome, A., & Grbac, I. (1998) The western green lizard *Lacerta (viridis) bilineata* Daudin, 1804 (Sauria: Lacertidae) Natura Croatica, 7, 363-369.

Urban, M. C., Tewksbury, J. J., & Sheldon, K. S. (2012) On a collision course: competition and dispersal differences create no-analogue communities and cause extinctions during climate change. Proceedings of the Royal Society B, 279, 2072-2080. https://doi.org/10.1098/rspb.2011.2367

van Oldenborgh, G. J., Doblas Reyes, F. J., Drijfhout, S. S., & Hawkins, E. (2013) Reliability of regional climate model trends. Environmental Research Letters, 8, 014055. https://doi.org/10.1088/1748-9326/8/1/014055

Wang, Z.-L., Zhang, D.-Y., & Wang, G. (2005) Does spatial structure facilitate coexistence of identical competitors? Ecological Modelling, 181, 17-23. https://doi.org/10.1016/j.ecolmodel.2004.06.020

Žagar, A., Carretero, M. A., Osojnik, N., Sillero, N., & Vrezec, A. (2015) A place in the sun: interspecific interference affects thermoregulation in coexisting lizards. Behavioral Ecology and Sociobiology, 69, 1127-1137. https://doi.org/10.1007/s00265-015-1927-8

Submitted: 5 December 2022

First decision: 25 January 2023

Accepted: 20 June 2023

Edited by Janet Franklin