# Genetic and demographic trends from rear to leading edge are explained by climate and forest cover in a cold adapted ectotherm

Andréaz Dupoué, Audrey Trochet, Murielle Richard, Mahaut Victoire Sorlin, Michael Guillon, Jules Teulieres-Quillet, Clément Vallé, Cyrielle Rault, Maud Berroneau, Matthieu Berroneau, et al.

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

## To cite this version:

Andréaz Dupoué, Audrey Trochet, Murielle Richard, Mahaut Victoire Sorlin, Michael Guillon, et al.. Genetic and demographic trends from rear to leading edge are explained by climate and forest cover in a cold adapted ectotherm. Diversity and Distributions, 2021, 27 (2), pp.267-281. 10.1111/ddi.13202. hal-03009458

HAL Id: hal-03009458

https://hal.science/hal-03009458

Submitted on 17 Nov 2020

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<footnote n="20">20 * Authors correspondence: andreaz.dupoue@gmail.com</footnote>

<footnote n="21">21 # Co-senior authors</footnote>

<footnote>22</footnote>

<footnote>23</footnote>
24 **Abstract**

25 **Aim:** Determining whether altitudinal shifts in species distributions leave molecular
26 footprints on wild populations along their range margins from rear to leading edge.

27 **Location:** South-west France.

28 **Methods:** We compared the demographic and genetic variation in 42 wild populations of the
29 Western oviparous subclade B2 of a cold adapted lizard (*Zootoca vivipara louislantzi*). These
30 populations can be divided into four ecological units across altitudinal clines in South-west
31 France (rear edge: <100 m, admixture zone: 100-500 m, continuous range: 500-1300 m and
32 leading edge: >1300 m above sea level).

33 **Results:** Within the rear edge were found the highest levels of inbreeding, genetic
34 differentiation and evidence of interrupted gene flow compared to central or colonising areas.
35 Within the leading edge, altitudinal range expansion occurred over the last centuries and
36 populations showed relatively low genetic diversity. These demographic and genetic trends
37 were better explained by inhospitable (warm and dry) climate conditions and forest cover.

38 **Main conclusions:** This empirical evidence illustrates that molecular footprints of climate
39 conditions and habitat quality on wild population trends can be perceived after recent events,
40 which should be of particular importance to accurately understand and anticipate human-
41 induced global change on wild species and ecosystems.

42 **Keywords.** Abundance, climate, colonization, ectotherm, forest cover, gene flow, inbreeding,
43 population decline, structuration

44
45 **Introduction**

46 Biodiversity is facing a global loss of habitats resulting in altered species and biotic
47 interactions worldwide (Ceballos et al., 2017; Plotnick et al., 2016; Wiens, 2016). Causes are
48 multiple, often synergetic, and have various footprints from genes to ecosystems (Pecl et al.,
49 2017). As a consequence, many species' distributions are shifting (Lenoir et al., 2020;
50 VanDerWal et al., 2013) and these range shifts have direct impacts on processes and patterns
51 of population genetics (Carvalho et al., 2019; Hampe & Petit, 2005; Templeton et al., 2001).
52 In particular, once populations at the rear edge become isolated, the loss of connectivity and
53 immigration progressively leads to an increase in inbreeding through genetic drift, one of the
54 main precursor of decline in small populations (Frankham, 2010; Hampe & Petit, 2005). In
55 contrast, at the leading edge, newly favourable habitats should be colonized by few
56 individuals resulting in founder events, which implies lower genetic diversity associated with
57 small groups of emigrants and recurrent extinction events (Hampe et al., 2013; Hampe &
58 Petit, 2005; Nadeau & Urban, 2019; Vilà-Cabrera et al., 2019; Waters et al., 2013). Studies
59 have challenged these general assumptions since genetic diversity does not necessarily
60 correlate with biogeography (Eckert et al., 2008; Pironon et al., 2017) nor does it explain
61 population decline (Tobler et al., 2013). This prompts for empirically testing patterns of
62 genetic variation and phylogenetic history along colonization gradients to better understand
63 how population genetic diversity or structure correlate with range shifts of wild populations.

64 Terrestrial ectotherms are well suited organisms for this purpose since their ecology
65 and demography tightly depend on climate conditions and landscape structure (Le Galliard et
66 al., 2012). These organisms are characterized by a strong thermal dependence on their
67 performance, which may explain their spatial distribution (Buckley et al., 2012). Under
68 temperate climates, ectotherm distributions at a continental scale have been strongly
69 influenced by past Quaternary climatic oscillations resulting in successive range expansions
and range contractions in glacial refugia (Hewitt, 2000). These historical changes must be considered to address current patterns of genetic variation at a regional scale, especially to characterise the recent impacts of climate change. Indeed, terrestrial ectotherms in their trailing range edge (low latitudes and low altitudes) are now increasingly challenged by recurrent temperatures nearing their upper thermal safety margin (Sunday et al., 2014).

Despite a relatively high potential of phenotypic plasticity, once temperatures rise above maximal critical temperatures, individuals at the rear edge of the distribution are exposed to greater risks of overheating and dehydration, which can impair the activity time and incur physiological costs (Dupoué et al., 2018, 2020; Le Galliard et al., 2012; Rozen-Rechels et al., 2019; Sinervo et al., 2010). These extreme climate conditions may constrain lifetime reproductive success and/or longevity, eventually leading to a loss of reproductive individuals and population extirpation (Bestion et al., 2015; Dupoué, Rutschmann, Le Galliard, Clobert, et al., 2017; Lourdais et al., 2017).

The present comparative study assessed relative abundance and genetic profiles in natural populations of the Western oviparous subclade B2 of the common lizard (Zootoca vivipara louislantzi) along an altitudinal cline. This boreal Squamate is a cold adapted specialist combining freezing tolerance to supercooling physiological adaptations (Voituron et al., 2002), as well as reproductive bimodality with distinct oviparous and viviparous populations across its distribution range (Surget-Groba et al., 2006). The viviparous form has colonized the Eurosiberian region from Western Europe to Hokkaido island in Japan being restricted by the polar circle (Surget-Groba et al., 2006). In the Southwest margin of its range, the oviparous form (Z. v. louislantzi) experiences a spatially restricted distribution that has been structured across altitude most likely due to temperatures and anthropogenic habitat fragmentation. Populations were affiliated to four ecological units defined a priori and following the classification established by Hampe & Petit (2005). This allows for the
95 examination of recolonization dynamics and how multiple demographic trajectories shape
96 genetic variation since the Last-Glacial-Maximum (LGM). In the rear edge within 100 m
97 above sea level (ASL), man-made draining of natural habitats (peatbogs and marshes)
98 occurred in the mid-19th century, which has isolated populations into a few suitable patches of
99 humid areas [boggy islets (Berroneau, 2014)]). Surrounding environments became
100 inhospitable because of human landscape modification and recent climate warming, thus
101 hindering any emigration from these populations. Similar anthropogenic-driven isolations
102 have occurred in the foothills of the French Pyrenees Mountains between 100 and 500 m
103 ASL, where populations are enclaved within forest-isolated peat bogs. This area is located
104 between populations from the rear edges and those of the continuous range, now forming the
105 Admixture zone [sensu (Hampe & Petit, 2005)], a lineage mixing area between previously
106 allopatric sub-clades (Rius & Darling, 2014). On the other hand, in its continuous range from
107 500 to 1,300 m ASL, populations occupy various habitat types mostly composed by river
108 ecosystems and forest borders. Finally, populations are now expanding above 1,300 m ASL
109 and colonizing highland meadows, which represent the leading edge.

110 We used microsatellite polymorphisms to assess how 42 natural populations across
111 this cline differed in genetic structure, diversity and connectivity. We anticipated that
112 Pyrenean ascension occurred a few times after the LGM, similar to the latitudinal expansion
113 of the viviparous form of this species (Horreo et al., 2018). Additionally, we examined how
114 human-induced environmental threats may additively affect lizard populations by testing the
115 specific influence of recent climate conditions and habitat characteristics on demographic and
116 genetic markers. Because rear edge populations are exposed to increasingly hostile (warm and
117 dry) thermal conditions, they should also experience high levels of differentiation and
118 inbreeding reflecting range contraction (Hampe & Petit, 2005). Conversely, in the wave front
119 of the expanding margin, we predicted higher segregation and lower genetic diversity if
colonisation occurs through founder events (Hampe et al., 2013; Nadeau & Urban, 2019; Waters et al., 2013). Finally, we expected demographic and genetic variation to be negatively impacted by constraining climate conditions and low habitat quality. We tested this prediction by checking the influence of two climate variables (i.e., air temperature and precipitation flux) and two habitats characteristics (i.e., forest cover and wetland potential) that are known to drive ecophysiological adjustments and local adaptations in this species (Dupoué, Rutschmann, Le Galliard, Clobert, et al., 2017; Dupoué, Rutschmann, Le Galliard, Miles, et al., 2017; Lorenzon et al., 1999; Rozen-Rechels et al., in press; Rutschmann et al., 2016, 2020).

129

## 130 Material and methods

Sampled populations among four ecological units

We caught 627 individuals from 42 populations between April-August 2017 and in May 2018. Populations were selected and considered different based on natural barriers (rivers, forests, altitude) and a minimal 2 km plane geographic distance (~300 m of altitude differential). Following capture, lizards were swabbed in the buccal cavity (Beebee, 2008), and samples were immediately stored in TE buffer and frozen at -20°C the same day to optimize DNA conservation until extraction. In each population, we estimated the relative abundance as the number of lizards captured divided by the time spent on the population and the number of people (range: 2 – 4 persons). As a heliothermic thermoregulator, common lizards rely on solar radiation to regulate body temperature, and capture probability may depend on weather conditions. We limited this bias by sampling all populations under relatively similar weather conditions (< 50% cloud cover) during the activity peak time (10:30 - 17:30) as long as basking opportunities remained available. In addition, cloud cover (0-50
144 %) had no statistical influence on abundance estimate (p = 0.626) and a similar index
145 provided a reliable proxy of absolute density in another study (Dupoué, Rutschmann, Le
146 Galliard, Clobert et al., 2017). Populations were divided into the four ecological units based
147 on their altitude thresholds as described above (100, 500, and 1,300 m) and the associated
148 landscapes (Table S1). This included the boggy islets in the rear edge (n = 55 lizards among 4
149 populations), the forest enclaves of the admixture zone (137 lizards among 9 populations), the
150 mixed habitats of the continuous range (242 lizards among 16 populations), and the highland
151 meadows in the leading edge (198 lizards among 13 populations).

152 *Climatic conditions and habitat characteristics in sampled populations*

153 For each population, we extracted two climate variables (temperature and precipitations) of
154 recent climate conditions from CHELSA database, where climatology is averaged between
155 1981 and 2005 (Karger et al., 2017, 2020) and is available at a high resolution of 30 Arcsec
156 (~1 km²). For each population, we considered the annual near-surface (2 m) air temperature
157 (Fig. S1A) and annual cumulative precipitation flux corrected by orographic effects (Fig.
158 S1B). We also considered two descriptors of habitat quality relevant for this species: forest
159 cover and wetland potential. We estimated forest cover at a precision of 30 m from the
160 aboveground live woody biomass density for the year 2000 (Fig. S1C), derived from many
161 studies (http://www.globalforestwatch.org/). Wetland potential (the probability for the habitat
162 to be a humid zone ranging from 0 to 3; 0: dry environments; 1: good probability of wetland
163 area; 2: strong probability of wetland area; 3: very strong probability of wetland area) was
164 obtained at a precision of 50 m (Fig. S1D), calculated from pedological and hydrogeological
165 information (Berthier et al., 2014). This index provides a continuous measure of the hydric
166 environment in each population, which fitted well with our direct observation in the field
167 (Dupoué, *pers. obs.*).

168 *Genetic data*
169 Genomic DNA was extracted from swabs using the DNeasy 96 Blood and Tissue Purification Kit (QiagenTM, Hilden, Germany) according to the manufacturer's instructions. We 170 amplified sixteen microsatellite markers (Table S2) in 5µl reaction volumes with 1µl of PCR 171 master Kit (MPBiomedicals) and 5-10ng of template DNA. PCR products were run on an ABI 172 3730 DNA Analyser (Applied Biosystems) with Genescan-600Liz size standard. Genotyping 173 was performed with Genemapper 5.0 (Applied Biosystems). We performed an initial 174 examination on the dataset to retain the loci with limited allele dropout (< 18% of missing 175 data). Ambiguous genotypes were amplified and sized a second time. Out of 16 initial loci, 5 176 loci were rejected from further analyses due to an excess of missing data (Table S2). Null 177 alleles and linkage disequilibrium among remaining loci was tested using GENEPOP 4.7.0 178 (Rousset, 2008) following 10,000 iterations. The 11 loci did not present significant evidence 179 for null alleles (Table S2), nor linkage disequilibrium (Table S3) and were used in the 180 following analyses. We used COLONY software (Jones & Wang, 2010) to identify and 181 remove relatives in each population (Table S4) given that they may alter subsequent analyses 182 (Rodríguez-Ramilo & Wang, 2012). The final dataset included n = 599 individuals. 183
184 *Statistical analyses*

185 *Genetic metrics*

186 We used R software (R Development Core Team 2015) and the *DiveRsity* R-package (Keenan 187 et al., 2013) to estimate genetic diversity in each population including expected 188 heterozygosity (He), observed heterozygosity (Ho), deviation test to Hardy-Weinberg 189 equilibrium (HWE), and allelic richness (Ar). We used the same package to estimate 190 heterozygote deficit within a population ($F_{IS}$) and genetic differentiation between populations 191 ($F_{ST}$) after 1,000 bootstraps. In all these analyses, we excluded one population (Lagunas, 192 Table S1) given that the sample size (n = 6 lizards) was too small to obtain reliable estimates 193 of population genetic diversity. We used the *poppr* R-package (Kamvar et al., 2015) to
194 analyse the molecular variance (AMOVA) and check how variation in genetic diversity
195 differed i) among ecological units, ii) among or within populations, and iii) within
196 individuals. We used the program BOTTLENECK (Piry et al., 1999) to check if any
197 populations experienced a recent bottleneck. Shortly after a bottleneck event, heterozygosity
198 increases more rapidly than allelic diversity (Piry et al., 1999). The program tests the
199 heterozygosity excess and compares Ho and He (Wilcoxon test) after 1,000 iterations. We
200 considered a two phase model including a 95% of stepwise mutation model and a variance of
201 12 as recommended (Piry et al., 1999).

202 Genetic structure

203 We used the Bayesian clustering analysis of the program STRUCTURE 2.3.3 (Pritchard et al.,
204 2000) to determine whether populations could be regrouped without a priori information. We
205 determined how many clusters (K) best fit with our data, based on ΔK method (Earl &
206 VonHoldt, 2012). To do so, we ran 10 independent simulations for each K ranging from 1 to
207 42, with 500,000 iterations and a burn-in period of 100,000 steps. We adjusted the burn-in
208 period to 200,000 steps in another simulation with K ranging from 1 to 10, to check that the
209 initial burn-in period was sufficiently long to ensure model convergence. We compared two
210 outputs (alpha and likelihood) from the program, illustrating relatively few admixtures
211 between population (alpha < 1) and good model convergence (Fig. S2). We estimated
212 structuration of each population as the highest probability (maximal Q-value) for an
213 individual to be affiliated to a cluster (Fig. 1).

214 We complemented this analysis by using another method for the examination of
215 population structure. We considered a method based on Discriminant Analysis of Principal
216 Component (DAPC) proposed by Jombart et al., (2010) and implemented in the *adegenet* R
217 package (Jombart, 2008). This multivariate method generates models for a range of *K* clusters
218 and calculate Bayesian Information Criterion (BIC) after 10,000 iterations (Jombart et al.,
219

2010). Optimal *K* is chosen following a stepping stones approach, for the minimum *K* when

220

BIC decreases by a negligible amount (Fig. S3).

221

## Genetic connectivity

222

We examined the degree of isolation-by-distance (IBD) using Mantel tests implemented in the *adegenet* R package. We compared the relationships between pairwise genetic distances and pairwise Euclidian geographic distances either globally (i.e., analysis on the entire dataset) or locally (i.e., analyses within each ecological unit). We considered chord genetic distances (*D*<sub>C</sub>) between populations defined by (Cavalli-Sforza & Edwards, 1967), since these are more likely to retrieve correct relations among populations (Takezaki & Nei, 1996). P-values were generated by comparing observed and predicted distributions after 1,000 iterations.

223

224

225

226

227

228

229

## Effect of environmental conditions on demography and genetic diversity

230

We used AICc based model selection and averaging from the *MuMln* R-package (Barton, 2019) to examine the additive impacts of climate conditions and habitat characteristics on lizard abundance (calculated from field data) and genetic diversity [calculated with *DiveRsity* R package (Keenan et al., 2013)]. We ran four independent sets of 15 linear mixed model [R package *nlme* (Pinheiro et al., 2016)] comparisons, using either demographic (lizard abundance) or genetic markers (Ar, *F*<sub>IS</sub> and *F*<sub>ST</sub>) as response variables, and environmental covariates alone or in addition (see models detailed in Table 1). Level of correlation between environmental variables was below 0.42, thus limiting the risk of multicollinearity issue (Dormann et al., 2013). In all models, ecological unit was set as a random term to account for spatial interdependence of populations. For each demographic or genetic parameter, we inferred the best models from the AICc score [ΔAICc < 2 (Burnham & Anderson, 2002)], and calculated the conditional model average to test for effect size and significance of fixed co-

231

232

233

234

235

236

237

238

239

240

241
242

variates. We verified that residuals of all selected models (Table 1) respected normality

243

conditions (Shapiro test: all p > 0.104).

244

Pyrenean colonisation scenarios

245

We used Approximate Bayesian Computation (ABC) to test different Pyrenean colonisation

246

scenarios using DIYABC 2.1.0 (Cornuet et al., 2010). First, we determined if variation in

247

population genetics within and between ecological units coincided to an ancient or a recent

248

Pyrenean colonisation. We followed the same methodology as described in (Ursenbacher et

249

al., 2015) and generated 3·10⁶ simulations to compare the posterior probabilities (logistic

250

regression) of three scenarios (Fig. 2) including: a simultaneous split among the four

251

ecological units (scenario 1), an altitudinal progressive colonisation (scenario 2) and a

252

simultaneous split in the Pyrenees regions differing from the rear edge (scenario 3). We

253

restricted the number of tested scenarios by ordering ecological units based on their altitudinal

254

class (Fig. 2). However, we did not test for scenarios with higher proximity between the rear

255

and the leading edge compared to others units since it has no biological meaning. Once the

256

most likely scenario was identified, we further ran 4·10⁶ iterations to estimate the parameters

257

(effective population size: N, and divergence time: t). We defined N1 < N2 < N3 > N4 and t3

258

> t2 > t1 as conditional priors (Fig. 2) and we considered the mean number of alleles, the

259

mean gene diversity, the mean size variance, and, between two samples, the FST, the mean

260

index of classification, the shared allele distance, and the (dµ)² distance as summary statistics.

261

We checked the performance of each ABC analysis (scenario comparison and parameters

262

estimation) by simulating 500 pseudo-observed datasets, comparing their posterior

263

distribution, and compiling relative error (Cornuet et al., 2010). We used a generation time of

264

3 years based on the equation T = α + [s/1(1-s)] (Lande et al., 2003) where α is the age at

265

maturity and s is the annual survival rate. Survival and age at maturity are likely to vary along

266

altitude (Heulin et al., 1997), so we considered adult survival rate of 0.51 (low and
267 intermediate altitudes) or 0.32 (high altitudes) and a recruitment in reproduction at one year
268 old (low and intermediate altitudes) or 2 years old (high altitudes).

269 **Results**

270 <u>Genetic diversity</u>

271 Analyses of molecular variance showed that genetic diversity varied mostly within individuals
272 (86.8%), but a substantial and significant part of variation was also explained within
273 populations (4.1%) and between populations (7.4%) followed by a minor variation between
274 ecological units (1.7%, see Table S5). Within population diversity and population
differentiation varied among ecological units (Ar: F<sub>3,37</sub> = 3.7, p = 0.020; F<sub>IS</sub>: F<sub>3,37</sub> = 4.6, p =
275 0.008; F<sub>ST</sub>: F<sub>3,37</sub> = 4.2, p = 0.012). Ar was the highest in the admixture zone and F<sub>IS</sub> were
276 highest in the rear edge and decreased along altitudinal clines (Fig. 3a & 3b). F<sub>ST</sub> was higher
277 in both the rear and leading edge but the lowest in the admixture zone (Fig. 3c). We did not
278 find any evidence for a recent bottleneck (p > 0.216), except in one population (BTM, p =
279 0.002).
280

281 <u>Genetic structure</u>

282 Both the Bayesian spatial modelling of genetic population structure and the ordination method
283 gathered populations into *K* = 6 clusters (Table S6 & Fig. S2). Clustering was relatively
284 heterogeneous between and within ecological units, and rather dependent on population
285 locality (Fig. 1). Yet, structuration (i.e., the probability for an individual to be affiliated to a
286 cluster) differed significantly between ecological units (F<sub>3,38</sub> = 3.5, p = 0.024) since there was
287 a lower structuration in the admixture zone compared to the 3 other ecological units (pairwise
288 comparisons, all p < 0.032).

289 <u>Genetic isolation</u>
290 Isolation by distance (IBD) was highly significant globally (p = 0.001, Fig. 4a). Locally
291 however, IBD differed among ecological units and according to isolation level. IBD was non-
292 significant at both the rear edge (p = 0.276, Fig. 4b) and the admixture zone (p = 0.320, Fig.
293 4c), where all pairwise geographic distances were below 60 km. Instead, significant IBD was
294 found in both the continuous range (p = 0.006, Fig. 4d) and the leading edge (p = 0.003, Fig.
295 4e) where many pairwise geographic distances were above 60 km.

296 <u>Effect of past and present environmental conditions on demography and genetic diversity</u>

297 Spatial variation in lizard abundance was correlated solely to air temperature (Table 1, Fig.
298 5a), with the most scattered populations associated to warm environments. Ar was the lowest
299 under warm and dry climates, and in open habitats (Table 1, Fig. 5b). Variation in both FIS
300 and FST were better explained by additive constraints in temperature and forest cover (Table
301 1, Fig. 5c & 5d). That is, inbreeding and differentiation is higher in warm and open habitats
302 (Table 1).

303 <u>Pyrenean colonisation scenarios</u>

304 The most likely scenario was the colonization scenario 2 with a posterior probability of 0.86 ±
305 0.04 (Fig. 2). We compared logistic distribution in 500 simulated datasets to estimate a global
306 posterior predictive error of 0.17, a type I error (probability of rejecting scenario 2 when it is
307 the correct scenario) of 0.33 and a type II error (probability of choosing scenario 2 when it is
308 an incorrect scenario) of 0.22. Effective population size considerably varied between
309 ecological units ranging from 1,500 in the leading edge to 9,000 reproductive individuals in
310 the continuous range (Table 2). The split between the rear edge and the other populations
311 occurred 750 generations ago, thereby placing the start of the Pyrenean ascension to around
312 2,000 years ago (Table 2). The divergence in the other ecological units occurred very recently
313

in the last millennium, with evidence of high elevation colonization between 100 and 700

314

years before present (Table 2).

315

316

Discussion

317

Genetic diversity was relatively high amongst lizards or across populations in the
geographic area and reasonably follows the theoretical predictions of the rear-leading edges
model of population genetic diversity and structure (Hampe et al., 2013; Hampe & Petit,
2005). That is, populations at the rear edge experience strong levels of inbreeding (loss of
heterozygosity). Associated with this, pairwise genetic distances were not correlated to
geographic distances neither within the rear edge nor in the admixture zone, probably due to
recently interrupted gene flow. From the continuous range to the leading edge, populations are
mostly isolated by distance, which suggests continuous distance-limited gene flow over time.
However, lower allelic richness and high differentiation in the leading edge populations imply
that they likely originated via small groups of immigrants through founder events (Hampe et
al., 2013; Nadeau & Urban, 2019; Waters et al., 2013). Lower inbreeding levels in these
colonized areas might for instance reflect a purge of inbred individuals during recruitment
(Hampe et al., 2013), since this species can optimise mate choice for that purpose (Richard et
al., 2009).

318

319

320

321

322

323

324

325

326

327

328

329

330

331

Our results provide new insights on genetic diversity and structuring of the subclade
B2 of *Zootoca vivipara*. Concomitantly with the geographic distribution of six genetic units
highlighted here, previous studies have shown that the strong genetic structuration of the
clade B is strongly shaped by mountain barriers and phylogeography (Horreo, Breedveld, et
al., 2019; Horreo, Peláez, et al., 2019; Milá et al., 2013). Genetic introgression between
southern Pyrenees (in Spain) and northern Pyrenees (in France) subclades is also common

332

333

334

335

336
among high altitude populations following secondary contacts associated with range shifts in warmer climates (Horreo, Breedveld, et al., 2019; Milá et al., 2013). This pattern of genetic introgression in the high elevation mountain range may explain why a population such as ARA sampled at high altitude and near a contact zone (Table S1) may have its own genetic structure. In our study, we further found that the highest allelic richness and lowest differentiation or structuration was found in the admixture zone, as this area may represent a mixing zone between genetic lineages of the rear edge populations and those from the continuous range (Hampe & Petit, 2005). Yet, mixing between populations arises despite obvious isolation and interrupted gene flow. Therefore, we hypothesize this specific pattern in the admixture zone to reflect ancient (structuration) versus recent (isolation) events caused by climate conditions and habitat fragmentation (Päckert et al., 2019). Surprisingly however, we found little evidence for bottlenecks (only one population likely due to a local isolation by a road), contrary to Horreo, Peláez, et al., (2019) who documented recent bottlenecks and significant heterozygote deficiencies in most clade B populations. This difference can be caused by methodological bias [n = 11 loci in the present study, n = 28 loci in (Horreo, Peláez, et al., 2019)] and a resulting lack of statistical power to detect recent bottlenecks. Further investigations including more microsatellites will help to better track recent genetic changes in the studied area.

A complementary study of the clade B populations further suggested that recent global warming may explain low migration rates and recent bottlenecks (Horreo, Peláez, et al., 2019). In support of a scenario of range contraction and expansion together with genetic structuring shaped by global changes, we showed here that climate conditions and forest cover explained additively variation in lizard demography and population genetics. Specifically, warmer temperatures in lowland populations were associated with lower lizard abundance and higher levels of inbreeding and genetic differentiation. Therefore, in addition
to relatively low effective population size, populations at the hot margin of the distribution
may be at greater risk of extinction (Cornetti et al., 2015; Templeton et al., 2001). Indeed,
comparative studies in contemporary populations of the viviparous form of this species have
shown that heat and drought represent important causal determinants of physiological stress,
reduction of behavioural activity and eventually population collapse (Bestion et al., 2015;
Dupoué et al., 2018, 2020; Dupoué, Rutschmann, Le Galliard, Clobert, et al., 2017; Massot et
al., 2008). Furthermore, forest cover systematically exhibited an additive influence on within-
population diversity, inbreeding, and genetic differentiation. That is, contemporary
populations located in the most open habitats exhibit higher inbreeding under warm climates
(low altitude) and experience poorer allelic richness and higher differentiation under cold
climates (high altitude). Although this species typically occupies open habitats (meadows,
peat-bogs, heathlands) to perform daily activities, mixed habitats with proximity to dense
forests might represent more intact habitat landscapes. In addition, lizards in populations next
to or mixed with forest patches are better capable of adjusting their breeding phenology
(Rutschmann et al., 2016) and their thermoregulation behaviour (Rozen-Rechels et al., in
press). Mixed forest habitats might therefore allow lizards to buffer negative consequences of
extreme weather events. Our population genetics study thus suggests that this cold-adapted
ectotherm may maintain small populations, especially at low altitudes, depending on
proximity to forest cover, possibly serving as 'thermal refuges' during resting periods. Hence,
the preservation of these habitats might represent a key conservation action to protect this
lizard and the associated cortege of species (Jofré et al., 2016). There is now a critical need
for developing genetically informed ecological niche models to better predict future species
distributions under cumulative constraints of climate change and landscape structuration
(Carvalho et al., 2019; Ikeda et al., 2017; Razgour et al., 2019).
Our results also revealed that the altitudinal colonisation of the Pyrenees started over the past two millennia (< 2 kya), much more recently than we expected (~10 kya following the LGM). Splits were particularly recent between the foothill and upward populations (~400 years), and between highland populations (~100 years), thus confirming a rapid and recent range expansion. Errors in scenario selection and parameters estimations were relatively low and in the range of acceptable error to be confident in observed results (Cornuet et al., 2010).

The most likely and parsimonious scenario for such recent ascension relates to the combined effects of habitat restructuration and recent climate change. Following the LGM, the present continuous range and leading edges were likely impassable for such a heliothermic lizard given that above 1,000 m, the Pyrenees were at the time entirely covered by dense pine forest (Davasse et al., 1996). As unravelled by palynology, this area was mostly deforested during the Middle Ages to develop pastoralism (Galop et al., 2003; Zanon et al., 2018), which probably opened new habitats and a new climbing path for this species. Additionally, climate conditions in the Pyrenees likely became suitable only after the Little Ice Age when the mountain glaciers retreated starting since approximately 1850. Our findings are consistent with similar insights related to the recent genetic differentiation (~200 years) in the same geographic area of the Pyrenean Desman (*Galemys pyrenaicus*), an endemic semiaquatic mammal affiliated to rivers and streams (Gillet et al., 2017).

Retracing population history of wild species is important if we are going to correctly interpret current distributions and foresee future impacts of global warming on demographic trends. In short, this study provides a rare example supporting general theoretical predictions stemming from the rear-leading edge conceptual framework (Hampe et al., 2013; Hampe & Petit, 2005; Nadeau & Urban, 2019; Razgour et al., 2013; Vilà-Cabrera et al., 2019; Waters et al., 2013). Rear edge populations showed the alarming lowest genetic diversity and highest inbreeding, interrupted gene flows and lower population size, which put them at greater risk
411 of extinction. Ground-dwelling ectotherms usually display relatively low dispersal rates
412 (Stevens et al., 2014), yet our study challenged this generality and suggested fast colonisation
413 capabilities in the common lizard likely facilitated by past habitat connectivity. Together, our
414 results therefore illustrate the importance of integrating genetic diversity, differentiation and
415 colonisation capabilities to accurately characterise and anticipate the potential footprints of
416 global change on wild species and ecosystems.

417 References

418 Barton, K. (2019). *MuMIn: Multi-Model Inference*. R Package Version 1.43.6. https://cran.r-project.org/package=MuMIn

419

420 Beebee, T. J. C. (2008). Buccal swabbing as a source of DNA from squamate reptiles.
*Conservation Genetics*, 9(4), 1087–1088. https://doi.org/10.1007/s10592-007-9464-2

421

422 Berroneau, M. (2014). Lézard vivipare *Zootoca vivipara* (Lichtenstein, 1823). In C. Nature
423 (Ed.), *Atlas des amphibiens et reptiles d'Aquitaine* (pp. 106–109). Association Cistude
424 Nature.

425

426 Berthier, L., Bardy, M., Chenu, J.-P., Guzmova, L., Laroche, B., Lehmann, S., Lemercier, B.,
427 Martin, M., Mérot, P., Squividant, H., Thiry, E., & Walter, C. (2014). *Enveloppes des
428 milieux potentiellement humides de la France métropolitaine - notice
429 d'accompagnement. Programme de modélisation des milieux potentiellement humides de
France*. Ministère d'Ecologie, du Développement Durable et de l'Energie.

425

426

427

428

429

430

431

432

433

Bestion, E., Teyssier, A., Richard, M., Clobert, J., & Cote, J. (2015). Live fast, die young: Experimental evidence of population extinction risk due to climate change. *PLoS Biology*, 13(10), e1002281. https://doi.org/10.1371/journal.pbio.1002281

Buckley, L. B., Hurlbert, A. H., & Jetz, W. (2012). Broad-scale ecological implications of
434

ectothermy and endothermy in changing environments. *Global Ecology and Biogeography*, 21(9), 873–885. https://doi.org/10.1111/j.1466-8238.2011.00737.x

435

436

Burnham, K. P., & Anderson, D. R. (2002). *Model selection and multimodel inference: a practical information-theoretic approach*, 2nd edn. Springer.

437

438

Carvalho, S. B., Torres, J., Tarroso, P., & Velo-Antón, G. (2019). Genes on the edge: A framework to detect genetic diversity imperiled by climate change. *Global Change Biology*, 25(12), 4034–4047. https://doi.org/10.1111/gcb.14740

439

440

441

Cavalli-Sforza, L. L., & Edwards, A. W. F. (1967). Phylogenetic analysis models and estimation procedures. *American Journal of Human Genetics*, 19(012), 233–257.

442

443

Ceballos, G., Ehrlich, P. R., & Dirzo, R. (2017). Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. *Proceedings of the National Academy of Sciences of the United States of America*, 114(30), E6089–E6096. https://doi.org/10.1073/pnas.1704949114

444

445

446

447

Cornetti, L., Ficetola, G. F., Hoban, S., & Vernesi, C. (2015). Genetic and ecological data reveal species boundaries between viviparous and oviparous lizard lineages. *Heredity*, 115(6), 517–526. https://doi.org/10.1038/hdy.2015.54

448

449

450

Cornuet, J. M., Ravigné, V., & Estoup, A. (2010). Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0). *BMC Bioinformatics*, 11, 401. https://doi.org/10.1186/1471-2105-11-401

451

452

453

Davasse, B., Galop, D., & Rendu, C. (1996). Paysages du Néolithique à nos jours dans les Pyrénées de l’est d’après l’écologie historique et l’archéologie pastorale. *XVIIe Rencontres Interna- Tionales d’archéologie et d’histoire d’Antibes*, 577–599.

454

455

456

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G.,
457 Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E.,
458 Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013).
459 Collinearity: A review of methods to deal with it and a simulation study evaluating their
460 performance. *Ecography*, 36(1), 027–046. https://doi.org/10.1111/j.1600-461 0587.2012.07348.x

462 Dupoué, A., Blaimont, P., Rozen-Rechels, D., Richard, M., Meylan, S., Clobert, J., Miles, D.
463 B., Martin, R., Decencière, B., Agostini, S., & Le Galliard, J.-F. (2020). Water
464 availability and temperature induce changes in oxidative status during pregnancy in a
465 viviparous lizard. *Functional Ecology*, 34(2), 475–485. https://doi.org/10.1111/1365-466 2435.13481

467 Dupoué, A., Rutschmann, A., Le Galliard, J.-F., Clobert, J., Angelier, F., Marciau, C., Ruault,
468 S., Miles, D. B., & Meylan, S. (2017). Shorter telomeres precede population extinction in
469 wild lizards. *Scientific Reports*, 7, 16976. https://doi.org/10.1038/s41598-017-17323-z

470 Dupoué, A., Rutschmann, A., Le Galliard, J.-F., Clobert, J., Blaimont, P., Sinervo, B., Miles,
471 D. B., Haussy, C., & Meylan, S. (2018). Reduction of baseline corticosterone secretion
472 correlates with climate warming and drying across wild lizard populations. *Journal of
473 Animal Ecology*, 87(5), 1331–1341. https://doi.org/10.1111/1365-2656.12843

474 Dupoué, A., Rutschmann, A., Le Galliard, J.-F., Miles, D. B., Clobert, J., DeNardo, D. F.,
475 Brusch, G. A., & Meylan, S. (2017). Water availability and environmental temperature
476 correlate with geographic variation in water balance in common lizards. *Oecologia*, 185,
477 561–571. https://doi.org/10.1007/s00442-017-3973-6

478 Earl, D. A., & VonHoldt, B. M. (2012). STRUCTURE HARVESTER: A website and
479 program for visualizing STRUCTURE output and implementing the Evanno method.
480 *Conservation Genetics Resources*, 4(2), 359–361. https://doi.org/10.1007/s12686-011-
481 9548-7

482 Eckert, C. G., Samis, K. E., & Lougheed, S. C. (2008). Genetic variation across species' geographical ranges: The central-marginal hypothesis and beyond. *Molecular Ecology*, 17(5), 1170–1188. https://doi.org/10.1111/j.1365-294X.2007.03659.x

485 Frankham, R. (2010). Inbreeding in the wild really does matter. *Heredity*, 104(2), 124–124. https://doi.org/10.1038/hdy.2009.155

487 Galop, D., Mazier, F., López Sáez, J. A., & Vannière, B. (2003). Palynologie et histoire des activités humaines en milieu montagnard. Bilan provisoire des recherches et nouvelles orientations méthodologiques sur le versant nord des Pyrénées. *Archéologie Du Midi Médiéval*, 21(1), 159–170. https://doi.org/10.3406/amime.2003.1407

491 Gillet, F., Cabria Garrido, M. T., Blanc, F., Fournier-Chambrillon, C., Némoz, M., Sourp, E., Vial-Novella, C., Zardoya, R., Aulagnier, S., & Michaux, J. R. (2017). Evidence of fine-scale genetic structure for the endangered Pyrenean desman (*Galemys pyrenaicus*) in the French Pyrenees. *Journal of Mammalogy*, 98(2), 523–532. https://doi.org/10.1093/jmammal/gyx002

496 Hampe, A., Pemonge, M. H., & Petit, R. J. (2013). Efficient mitigation of founder effects during the establishment of a leading-edge oak population. *Proceedings of the Royal Society B: Biological Sciences*, 280(1764), 20131070. https://doi.org/10.1098/rspb.2013.1070

500 Hampe, A., & Petit, R. J. (2005). Conserving biodiversity under climate change: the rear edge matters. *Ecology Letters*, 8(5), 461–467. https://doi.org/10.1111/j.1461-0248.2005.00739.x

503 Heulin, B., Osenegg-Leconte, K., & Michel, D. (1997). Demography of a bimodal
504

reproductive species of lizard (*Lacerta vivipara*): survival and density characteristics of
oviparous populations. *Herpetologica*, 53(4), 432–444.

505

506

Hewitt, G. (2000). The genetic legacy of the quaternary ice ages. *Nature*, 405(6789), 907–913. https://doi.org/10.1038/35016000

507

508

Horreo, J. L., Breedveld, M. C., Lindtke, D., Heulin, B., Surget-Groba, Y., & Fitze, P. S. (2019). Genetic introgression among differentiated clades is lower among clades exhibiting different parity modes. *Heredity*, 123(2), 264–272.
https://doi.org/10.1038/s41437-019-0201-7

509

510

511

512

Horreo, J. L., Peláez, M. L., Breedveld, M. C., Suárez, T., Urieta, M., & Fitze, P. S. (2019). Population structure of the oviparous South-West European common lizard. *European Journal of Wildlife Research*, 65(1), 11. https://doi.org/10.1007/s10344-018-1242-6

513

514

515

Horreo, J. L., Pelaez, M. L., Suárez, T., Breedveld, M. C., Heulin, B., Surget-Groba, Y., Oksanen, T. A., & Fitze, P. S. (2018). Phylogeography, evolutionary history and effects of glaciations in a species (*Zootoca vivipara*) inhabiting multiple biogeographic regions. *Journal of Biogeography*, 45(7), 1616–1627. https://doi.org/10.1111/jbi.13349

516

517

518

519

Ikeda, D. H., Max, T. L., Allan, G. J., Lau, M. K., Shuster, S. M., & Whitham, T. G. (2017). Genetically informed ecological niche models improve climate change predictions. *Global Change Biology*, 23(1), 164–176. https://doi.org/10.1111/gcb.13470

520

521

522

Jofré, G. M., Warn, M. R., & Reading, C. J. (2016). The role of managed coniferous forest in the conservation of reptiles. *Forest Ecology and Management*, 362, 69–78.
https://doi.org/10.1016/j.foreco.2015.11.044

523

524

525

Jombart, T. (2008). adegenet: A R package for the multivariate analysis of genetic markers.
*Bioinformatics*, 24(11), 1403–1405. https://doi.org/10.1093/bioinformatics/btn129

526
527

Jombart, T., Devillard, S., & Balloux, F. (2010). Discriminant analysis of principal

528

components: a new method for the analysis of genetically structured populations. BMC

529

Genetics, 11, 94. https://doi.org/10.1186/1471-2156-11-94

530

Jones, O. R., & Wang, J. (2010). COLONY: A program for parentage and sibship inference

531

from multilocus genotype data. Molecular Ecology Resources, 10(3), 551–555.

532

https://doi.org/10.1111/j.1755-0998.2009.02787.x

533

Kamvar, Z. N., Brooks, J. C., & Grünwald, N. J. (2015). Novel R tools for analysis of

534

genome-wide population genetic data with emphasis on clonality. Frontiers in Genetics,

535

6(JUN), 208. https://doi.org/10.3389/fgene.2015.00208

536

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W.,

537

Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high

538

resolution for the earth’s land surface areas. Scientific Data, 4, 170122.

539

https://doi.org/10.1038/sdata.2017.122

540

Karger, D. N., Dabaghchian, B., Lange, S., Thuiller, W., Zimmermann, N. E., & Graham, C.

541

H. (2020). High resolution climate data for Europe. In EnviDat (Version 1.0).

542

https://doi.org/doi:10.16904/envidat.150

543

Keenan, K., Mcginnity, P., Cross, T. F., Crozier, W. W., & Prodöhl, P. A. (2013). DiveRsity:

544

An R package for the estimation and exploration of population genetics parameters and

545

their associated errors. Methods in Ecology and Evolution, 4(8), 782–788.

546

https://doi.org/10.1111/2041-210X.12067

547

Lande, R., Engen, S., & Saether, B. E. (2003). Stochastic population dynamics in ecology and

548

conservation. Oxford University Press.

549

Le Galliard, J.-F., Massot, M., Baron, J. P., & Clobert, J. (2012). Ecological effects of climate
change on European reptiles. In J. Brodie, E. Post, & D. Doak (Eds.), *Wildlife conservation in a changing climate* (pp. 179–203). University of Chicago Press.

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

Lorenzon, P., Clobert, J., Oppliger, A., & John-Alder, H. (1999). Effect of water constraint on growth rate, activity and body temperature of yearling common lizard (*Lacerta vivipara*). *Oecologia*, 118(4), 423–430. https://doi.org/10.1007/s004420050744

Lourdais, O., Dupoué, A., Guillon, M., Guiller, G., Michaud, B., & DeNardo, D. F. (2017). Hydric “costs” of reproduction: pregnancy increases evaporative water loss in the snake *Vipera aspis*. *Physiological and Biochemical Zoology*, 90(6), 663–672.
https://doi.org/10.1086/694848

Massot, M., Clobert, J., & Ferrière, R. (2008). Climate warming, dispersal inhibition and extinction risk. *Global Change Biology*, 14(3), 461–469. https://doi.org/10.1111/j.1365-2486.2007.01514.x

Milá, B., Surget-Groba, Y., Heulin, B., Gosá, A., & Fitze, P. S. (2013). Multilocus phylogeography of the common lizard *Zootoca vivipara* at the Ibero-Pyrenean suture zone reveals lowland barriers and high-elevation introgression. *BMC Evolutionary Biology*, 13(1), 192. https://doi.org/10.1186/1471-2148-13-192

Nadeau, C. P., & Urban, M. C. (2019). Eco-evolution on the edge during climate change. *Ecography*, 42(7), 1280–1297. https://doi.org/10.1111/ecog.04404

Päckert, M., Ait Belkacem, A., Wolfgramm, H., Gast, O., Canal, D., Giacalone, G., Lo Valvo, M., Vamberger, M., Wink, M., Martens, J., & Stuckas, H. (2019). Genetic admixture
despite ecological segregation in a North African sparrow hybrid zone (Aves, 573 *Passeriformes, Passer domesticus × Passer hispaniolensis*). *Ecology and Evolution*, 574 9(22), 12710–12726. https://doi.org/10.1002/ece3.5744

Pecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen, I.-C., Clark, T. 576 D., Colwell, R. K., Danielsen, F., Evengård, B., Falconi, L., Ferrier, S., Frusher, S., 577 Garcia, R. A., Griffis, R. B., Hobday, A. J., Janion-Scheepers, C., Jarzyna, M. A., 578 Jennings, S., ... Williams, S. E. (2017). Biodiversity redistribution under climate change: 579 Impacts on ecosystems and human well-being. *Science*, 355(6332), eaai9214. 580 https://doi.org/10.1126/science.aai9214

Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., & R, C. T. (2016). *nlme: Linear and 582 Nonlinear Mixed Effects Models* (R package version 3.1-128). http://cran.r-project.org/package=nlme 583

Pironon, S., Papuga, G., Villellas, J., Angert, A. L., García, M. B., & Thompson, J. D. (2017). 585 Geographic variation in genetic and demographic performance: new insights from an old 586 biogeographical paradigm. *Biological Reviews*, 92(4), 1877–1909. 587 https://doi.org/10.1111/brv.12313 588

Piry, S., Luikart, G., & Cornuet, J.-M. (1999). BOTTLENECK: A computer program for 589 detecting recent reductions in the effective population size using allele frequency data. 590 *Journal of Heredity*, 90(4), 502–503. https://doi.org/10.1093/jhered/90.4.502

Plotnick, R. E., Smith, F. A., & Lyons, S. K. (2016). The fossil record of the sixth extinction. 592 *Ecology Letters*, 19(5), 546–553. https://doi.org/10.1111/ele.12589 593

Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using 594 multilocus genotype data. *Genetics*, 155(2), 945–959. 595
596 Razgour, O., Forester, B., Taggart, J. B., Bekaert, M., Juste, J., Ibáñez, C., Puechmaille, S. J.,
597 Novella-Fernandez, R., Alberdi, A., & Manel, S. (2019). Considering adaptive genetic
598 variation in climate change vulnerability assessment reduces species range loss
599 projections. *Proceedings of the National Academy of Sciences of the United States of
600 America*, **116**(21), 10418–10423. https://doi.org/10.1073/pnas.1820663116

601 Razgour, O., Juste, J., Ibáñez, C., Kiefer, A., Rebelo, H., Puechmaille, S. J., Arlettaz, R.,
602 Burke, T., Dawson, D. A., Beaumont, M., & Jones, G. (2013). The shaping of genetic
603 variation in edge-of-range populations under past and future climate change. *Ecology
604 Letters*, **16**(10), 1258–1266. https://doi.org/10.1111/ele.12158

605 Richard, M., Losdat, S., Lecomte, J., De Fraipont, M., & Clobert, J. (2009). Optimal level of
606 inbreeding in the common lizard. *Proceedings of the Royal Society B: Biological
607 Sciences*, **276**(1668), 2779–2786. https://doi.org/10.1098/rspb.2009.0319

608 Rius, M., & Darling, J. A. (2014). How important is intraspecific genetic admixture to the
609 success of colonising populations? *Trends in Ecology & Evolution*, **29**(4), 233–242.
610 https://doi.org/10.1016/j.tree.2014.02.003

611 Rodríguez-Ramilo, S. T., & Wang, J. (2012). The effect of close relatives on unsupervised
612 Bayesian clustering algorithms in population genetic structure analysis. *Molecular
613 Ecology Resources*, **12**(5), 873–884. https://doi.org/10.1111/j.1755-0998.2012.03156.x

614 Rousset, F. (2008). GENEPOP'007: A complete re-implementation of the GENEPOP
615 software for Windows and Linux. *Molecular Ecology Resources*, **8**(1), 103–106.
616 https://doi.org/10.1111/j.1471-8286.2007.01931.x

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

621 Rozen-Rechels, D., Rutschmann, A., Dupoué, A., Blaimont, P., Chauveau, V., Miles, D. B.,

622 Guillon, M., Richard, M., Badiane, A., Meylan, S., Clobert, J., & Le Galliard, J. (n.d.).

623 Interaction of hydric and thermal conditions drive geographic variation in thermoregulation in a widespread lizard. *Ecological Monographs*.

625 Rutschmann, A., Miles, D. B., Le Galliard, J.-F., Richard, M., Moulherat, S., Sinervo, B.,

626 Clobert, J. (2016). Climate and habitat interact to shape the thermal reaction norms of breeding phenology across lizard populations. *Journal of Animal Ecology*, 85(2), 457–466. https://doi.org/10.1111/1365-2656.12473

629 Rutschmann, A., Rozen-Rechels, D., Dupoué, A., Blaimont, P., de Villemereuil, P., Miles, D.

630 B., Richard, M., & Clobert, J. (2020). Climate dependent heating efficiency in the common lizard. *Ecology and Evolution*, 10(15), 8007–8017.

632 https://doi.org/10.1002/ece3.6241

633 Sinervo, B., Mendez-de-la-Cruz, F., Miles, D. B., Heulin, B., Bastiaans, E., Villagrán-Santa

634 Cruz, M., Lara-Resendiz, R., Martinez-Mendez, N., Calderón-Espinosa, M. L., Meza-

635 Lázaro, R. N., Gadsden, H., Avila, L. J., Morando, M., De la Riva, I. J., Victoriano

636 Sepulveda, P., Rocha, C. F. D., Ibargüengoytía, N. R., Aguilar Puntriano, C., Massot, M., ... Sites, J. W. (2010). Erosion of lizard diversity by climate change and altered thermal niches. *Science*, 328, 894–899. https://doi.org/10.1126/science.1184695

639 Stevens, V. M., Whitmee, S., Le Galliard, J.-F., Clobert, J., Böhning-Gaese, K., Bonte, D.,

640 Brändle, M., Matthias Dehling, D., Hof, C., Trochet, A., & Baguette, M. (2014). A

641 comparative analysis of dispersal syndromes in terrestrial and semi-terrestrial animals.

642 *Ecology Letters*, 17(8), 1039–1052. https://doi.org/10.1111/ele.12303

643 Sunday, J. M., Bates, A. E., Kearney, M. R., Colwell, R. K., Dulvy, N. K., Longino, J. T., &
644 Huey, R. B. (2014). Thermal-safety margins and the necessity of thermoregulatory
645 behavior across latitude and elevation. *Proceedings of the National Academy of Sciences
646 of the United States of America, 111*(15), 5610–5615.
647 https://doi.org/10.1073/pnas.1316145111

648 Surget-Groba, Y., Heulin, B., Guillaume, C.-P., Puky, M., Semenov, D., Orlova, V.,
649 Kupriyanova, L., Ghira, I., & Smajda, B. (2006). Multiple origins of viviparity, or
650 reversal from viviparity to oviparity? The European common lizard (Zootoca vivipara,
651 Lacertidae) and the evolution of parity. *Biological Journal of the Linnean Society,
652 87(1996), 1–11. https://doi.org/10.1111/j.1095-8312.2006.00552.x

653 Takezaki, N., & Nei, M. (1996). Genetic distances and reconstruction of phylogenetic trees
654 from microsatellite DNA. *Genetics*, **144**(October), 389–399.

655 Templeton, A. R., Robertson, R. J., Brisson, J., & Strasburg, J. (2001). Disrupting
656 evolutionary processes: The effect of habitat fragmentation on collared lizards in the
657 Missouri Ozarks. *Proceedings of the National Academy of Sciences of the United States
658 of America*, 98(10), 5426–5432. https://doi.org/10.1073/pnas.091093098

659 Tobler, U., Garner, T. W. J., & Schmidt, B. R. (2013). Genetic attributes of midwife toad
660 (*Alytes obstetricans*) populations do not correlate with degree of species decline. *Ecology
661 and Evolution*, 3(9), 2806–2819. https://doi.org/10.1002/ece3.677

662 Ursenbacher, S., Guillon, M., Cubizolle, H., Dupoué, A., Blouin-Demers, G., & Lourdais, O.
663 (2015). Postglacial recolonization in a cold climate specialist in western Europe: patterns
664 of genetic diversity in the adder (*Vipera berus*) support the central–marginal hypothesis.
665 *Molecular Ecology*, 24(14), 3639–3651. https://doi.org/10.1111/mec.13259

666 VanDerWal, J., Murphy, H. T., Kutt, A. S., Perkins, G. C., Bateman, B. L., Perry, J. J., &
667 Reside, A. E. (2013). Focus on poleward shifts in species’ distribution underestimates
668 the fingerprint of climate change. *Nature Climate Change*, 3(3), 239–243.

669 https://doi.org/10.1038/nclimate1688

670 Vilà-Cabrera, A., Premoli, A. C., & Jump, A. S. (2019). Refining predictions of population

671 decline at species’ rear edges. *Global Change Biology*, 25(5), 1549–1560.

672 https://doi.org/10.1111/gcb.14597

673 Voituron, Y., Mouquet, N., De Mazancourt, C., & Clobert, J. (2002). To freeze or not to

674 freeze? An evolutionary prespective on the cold hardiness strategies of overwinrering

675 ectotherms. *The American Naturalist*, 160(2), 255–270. https://doi.org/10.1086/341021

676 Waters, J. M., Fraser, C. I., & Hewitt, G. M. (2013). Founder takes all: Density-dependent

677 processes structure biodiversity. *Trends in Ecology & Evolution*, 28(2), 78–85.

678 https://doi.org/10.1016/j.tree.2012.08.024

679 Wiens, J. J. (2016). Climate-related local extinctions are already widespread among plant and

680 animal species. *PLoS Biology*, 14(12), e2001104.

681 https://doi.org/10.1371/journal.pbio.2001104

682 Zanon, M., Davis, B. A. S., Marquer, L., Brewer, S., & Kaplan, J. O. (2018). European forest

683 cover during the past 12,000 years: A palynological reconstruction based on modern

684 analogs and remote sensing. *Frontiers in Plant Science*, 9(March), 253.

685 https://doi.org/10.3389/fpls.2018.00253

686

687
688 Acknowledgements

689 We are particularly grateful to the associations ‘Nature en Occitanie’ and ‘Cistude Nature’
690 who provided the database that allowed us to localize the great majority of lizard populations,
691 and we thank the ‘Parc National des Pyrenees’ for authorizing us to sample lizards. We also
692 thank Dr Sylvain Ursenbacher for advice on genetic analyses.

693 **Funding**

694 This study received the financial and technical support from the Centre National de la
695 Recherche Scientifique (CNRS), and was funded by the Interreg Poctefa under the project
696 ‘Ectopyr’ (no. EFA031/15 to FA), by the Agence Nationale de la Recherche under the
697 ‘Aquatherm’ project (ANR-17-CE02-0013 to JFLG), by the European Union (Feder), the
698 regional council of Nouvelle Aquitaine, and the departmental council of Gironde & Pyrénées-
699 Atlantiques under the ‘Climate Sentinels’ project (coordinators: Cistude Nature, Fanny
700 Mallard), and under the ‘Aquastress’ project (2018-1R20214 to OL).

701 At the time of writing, AD was financially supported by the Australian Research Council (no.
702 DP180100058 to A/Prof Anne Peters).

703 Ethics

704 The collection of samples was approved by the "Conseil Scientifique Régional du Patrimoine
705 Naturel (CSRPN)" of the region Occitanie. Procedures complied with current French
706 regulations permit (#2017-s-02) relating to an authorization of capture, marking, transport,
707 detention, use and release of protected lizard species.

708 Author contributions

709 AD and AT conceived the study together with JC and FA. AD, MR, MS, JTQ, CV, CR,
710 MauB, MatB, OL and PB captured lizards. MR and MS performed molecular assays. AD, AT
and MG performed the data analyses. RB and MG participated in environmental data
extraction. MauB, MatB and GP initially identified lizard populations sampled in this study.
AT, RB, OC, OG, HLC, JS, JFLG, JC, and FA were involved in funding and planning this
work. AD led the sampling design and draft the manuscript. All authors commented the
article and gave final approval for publication.

716 **Biosketch**

Andréaz Dupoué is interested in depicting which are the proximate mechanisms driving
population extirpation. His main researches at the interface between molecular ecology and
ecological physiology emphasize the importance of water based trade-offs in the evolution of
wild squamate reptiles.

The research team lead multi-specific studies on Pyrenean and European ectotherms.

723 **Data accessibility**

The data supporting this article are available at the Zenodo Repository at
http://doi.org/10.5281/zenodo.4266021
<table><caption>Table 1. Model selection comparing the impacts of climate (temperature and precipitations) and habitat quality (forest cover and wetland potential) on lizard demography (abundance) and genetic profile [allelic richness (Ar), inbreeding (FIS) and differentiation (FST)]. For each population parameter, the table reports the ΔAICc based linear mixed model comparison and full averages of fixed covariates. Ecological unit was set as a random term in all models to account for spatial interdependence.</caption><thead><tr><th>Abundance</th><th>model</th><th>df</th><th>logLik</th><th>AICc</th><th>ΔAICc</th><th>w<sub>i</sub></th><th>term</th><th>estimate</th><th>SE</th><th>z value</th><th>Pr(&gt;|z|)</th></tr></thead><tbody><tr><td></td><td><strong>temperature</strong></td><td><strong>4</strong></td><td><strong>-73.9</strong></td><td><strong>156.87</strong></td><td><strong>0</strong></td><td><strong>0.45</strong></td><td>(Intercept)</td><td>-0.005</td><td>0.254</td><td>0.02</td><td>0.986</td></tr><tr><td></td><td>forest + temperature</td><td>5</td><td>-73.73</td><td>159.12</td><td>2.25</td><td>0.15</td><td><strong>temperature</strong></td><td><strong>-0.476</strong></td><td><strong>0.096</strong></td><td><strong>4.80</strong></td><td><strong>0.000</strong></td><td>***</td></tr><tr><td></td><td>temperature + wetland</td><td>5</td><td>-73.79</td><td>159.25</td><td>2.38</td><td>0.14</td><td>forest</td><td>-0.002</td><td>0.004</td><td>0.54</td><td>0.593</td></tr><tr><td></td><td>precipitation + temperature</td><td>5</td><td>-73.89</td><td>159.44</td><td>2.57</td><td>0.12</td><td>wetland</td><td>0.134</td><td>0.290</td><td>0.45</td><td>0.655</td></tr><tr><td></td><td>forest + temperature + wetland</td><td>6</td><td>-73.65</td><td>161.71</td><td>4.84</td><td>0.04</td><td>precipitations</td><td>0.056</td><td>0.309</td><td>0.18</td><td>0.860</td></tr><tr><td></td><td>forest + precipitation + temperature</td><td>6</td><td>-73.72</td><td>161.83</td><td>4.96</td><td>0.04</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>precipitation + temperature + wetland</td><td>6</td><td>-73.75</td><td>161.9</td><td>5.02</td><td>0.04</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>forest + precipitation + temperature + wetland</td><td>7</td><td>-73.62</td><td>164.53</td><td>7.66</td><td>0.01</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>(Null)</td><td>3</td><td>-79.21</td><td>165.06</td><td>8.19</td><td>0.01</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>precipitation</td><td>4</td><td>-78.75</td><td>166.59</td><td>9.72</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>wetland</td><td>4</td><td>-78.96</td><td>167</td><td>10.13</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>forest</td><td>4</td><td>-79.03</td><td>167.14</td><td>10.27</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>precipitation + wetland</td><td>5</td><td>-78.35</td><td>168.36</td><td>11.49</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>forest + precipitation</td><td>5</td><td>-78.71</td><td>169.09</td><td>12.22</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>forest + wetland</td><td>5</td><td>-78.74</td><td>169.16</td><td>12.28</td><td>0</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
<table><tbody><tr><td></td><td>forest + precipitation + wetland</td><td>6</td><td>-78.32</td><td>171.03</td><td>14.16</td><td>0</td></tr><tr><td>Ar</td><td>model</td><td>df</td><td>logLik</td><td>AICc</td><td>ΔAICc</td><td>w<sub>i</sub></td></tr><tr><td></td><td><strong>forest + precipitation + temperature + wetland</strong></td><td><strong>7</strong></td><td><strong>-12.03</strong></td><td><strong>41.45</strong></td><td><strong>0</strong></td><td><strong>0.38</strong></td></tr><tr><td></td><td><strong>forest + precipitation + temperature</strong></td><td><strong>6</strong></td><td><strong>-13.54</strong></td><td><strong>41.55</strong></td><td><strong>0.1</strong></td><td><strong>0.36</strong></td></tr><tr><td></td><td>precipitation + temperature + wetland</td><td>6</td><td>-14.65</td><td>43.76</td><td>2.31</td><td>0.12</td></tr><tr><td></td><td>precipitation + temperature</td><td>5</td><td>-16.11</td><td>43.93</td><td>2.48</td><td>0.11</td></tr><tr><td></td><td>forest + temperature + wetland</td><td>6</td><td>-16.5</td><td>47.48</td><td>6.03</td><td>0.02</td></tr><tr><td></td><td>temperature + wetland</td><td>5</td><td>-18.63</td><td>48.97</td><td>7.52</td><td>0.01</td></tr><tr><td></td><td>forest + temperature</td><td>5</td><td>-19.2</td><td>50.11</td><td>8.66</td><td>0</td></tr><tr><td></td><td>temperature</td><td>4</td><td>-21.24</td><td>51.6</td><td>10.15</td><td>0</td></tr><tr><td></td><td>forest + precipitation</td><td>5</td><td>-20.84</td><td>53.39</td><td>11.94</td><td>0</td></tr><tr><td></td><td>forest + precipitation + wetland</td><td>6</td><td>-19.77</td><td>54.01</td><td>12.56</td><td>0</td></tr><tr><td></td><td>wetland</td><td>4</td><td>-22.59</td><td>54.29</td><td>12.83</td><td>0</td></tr><tr><td></td><td>forest + wetland</td><td>5</td><td>-21.39</td><td>54.49</td><td>13.04</td><td>0</td></tr><tr><td></td><td>precipitation</td><td>4</td><td>-22.85</td><td>54.82</td><td>13.37</td><td>0</td></tr><tr><td></td><td>forest</td><td>4</td><td>-22.93</td><td>54.97</td><td>13.52</td><td>0</td></tr><tr><td></td><td>(Null)</td><td>3</td><td>-24.23</td><td>55.12</td><td>13.67</td><td>0</td></tr><tr><td></td><td>precipitation + wetland</td><td>5</td><td>-21.79</td><td>55.29</td><td>13.84</td><td>0</td></tr><tr><td>F<sub>IS</sub></td><td>model</td><td>df</td><td>logLik</td><td>AICc</td><td>ΔAICc</td><td>w<sub>i</sub></td></tr><tr><td></td><td><strong>forest + temperature</strong></td><td><strong>5</strong></td><td><strong>44.26</strong></td><td><strong>-76.8</strong></td><td><strong>0</strong></td><td><strong>0.44</strong></td></tr></tbody></table>

<table><thead><tr><th>term</th><th>estimate</th><th>SE</th><th>z value</th><th>Pr(&gt;|z|)</th></tr></thead><tbody><tr><td>(Intercept)</td><td>-0.009</td><td>0.057</td><td>0.16</td><td>0.875</td></tr><tr><td><strong>temperature</strong></td><td><strong>0.117</strong></td><td><strong>0.026</strong></td><td><strong>4.27</strong></td><td><strong>0.000</strong> ***</td></tr><tr><td><strong>precipitations</strong></td><td><strong>0.237</strong></td><td><strong>0.077</strong></td><td><strong>2.97</strong></td><td><strong>0.003</strong> **</td></tr><tr><td><strong>forest</strong></td><td><strong>0.002</strong></td><td><strong>0.001</strong></td><td><strong>2.13</strong></td><td><strong>0.033</strong> *</td></tr><tr><td>wetland</td><td>-0.122</td><td>0.073</td><td>1.61</td><td>0.108</td></tr></tbody></table>

<table><thead><tr><th>term</th><th>estimate</th><th>SE</th><th>z value</th><th>Pr(&gt;|z|)</th></tr></thead><tbody><tr><td>(Intercept)</td><td>0.008</td><td>0.014</td><td>0.52</td><td>0.606</td></tr></tbody></table>
<table><tr><td>forest + temperature + wetland</td><td>6</td><td>44.86</td><td>-75.24</td><td>1.55</td><td>0.2</td></tr><tr><td>forest + precipitation + temperature</td><td>6</td><td>44.76</td><td>-75.05</td><td>1.75</td><td>0.18</td></tr><tr><td>forest + precipitation + temperature + wetland</td><td>7</td><td>45.14</td><td>-72.89</td><td>3.91</td><td>0.06</td></tr><tr><td>forest + precipitation</td><td>5</td><td>41.59</td><td>-71.47</td><td>5.32</td><td>0.03</td></tr><tr><td>forest</td><td>4</td><td>39.92</td><td>-70.72</td><td>6.08</td><td>0.02</td></tr><tr><td>forest + wetland</td><td>5</td><td>41.22</td><td>-70.72</td><td>6.08</td><td>0.02</td></tr><tr><td>forest + precipitation + wetland</td><td>6</td><td>42.28</td><td>-70.09</td><td>6.71</td><td>0.02</td></tr><tr><td>temperature</td><td>4</td><td>39.04</td><td>-68.98</td><td>7.82</td><td>0.01</td></tr><tr><td>temperature + wetland</td><td>5</td><td>39.48</td><td>-67.25</td><td>9.55</td><td>0</td></tr><tr><td>precipitation + temperature</td><td>5</td><td>39.32</td><td>-66.93</td><td>9.87</td><td>0</td></tr><tr><td>(Null)</td><td>3</td><td>36.42</td><td>-66.2</td><td>10.6</td><td>0</td></tr><tr><td>precipitation</td><td>4</td><td>37.31</td><td>-65.51</td><td>11.29</td><td>0</td></tr><tr><td>wetland</td><td>4</td><td>37.26</td><td>-65.41</td><td>11.39</td><td>0</td></tr><tr><td>precipitation + temperature + wetland</td><td>6</td><td>39.7</td><td>-64.94</td><td>11.86</td><td>0</td></tr><tr><td>precipitation + wetland</td><td>5</td><td>38.03</td><td>-64.34</td><td>12.46</td><td>0</td></tr></table>

<table><thead><tr><th>F<sub>ST</sub></th><th>model</th><th>df</th><th>logLik</th><th>AICc</th><th>ΔAICc</th><th>w<sub>i</sub></th></tr></thead><tbody><tr><td></td><td><strong>forest + precipitation + temperature</strong></td><td><strong>6</strong></td><td><strong>98.41</strong></td><td><strong>-182.34</strong></td><td><strong>0</strong></td><td><strong>0.38</strong></td></tr><tr><td></td><td><strong>forest + precipitation + temperature + wetland</strong></td><td><strong>7</strong></td><td><strong>98.89</strong></td><td><strong>-180.38</strong></td><td><strong>1.97</strong></td><td><strong>0.14</strong></td></tr><tr><td></td><td>forest + temperature</td><td>5</td><td>95.95</td><td>-180.19</td><td>2.16</td><td>0.13</td></tr><tr><td></td><td>forest + temperature + wetland</td><td>6</td><td>97.01</td><td>-179.56</td><td>2.79</td><td>0.09</td></tr></tbody></table>

<table><tr><td>temperature</td><td>0.017</td><td>0.006</td><td>2.67</td><td>0.008</td><td>**</td></tr><tr><td>forest</td><td>-0.001</td><td>0.000</td><td>3.19</td><td>0.001</td><td>**</td></tr><tr><td>wetland</td><td>0.018</td><td>0.018</td><td>1.00</td><td>0.317</td><td></td></tr><tr><td>precipitations</td><td>-0.019</td><td>0.019</td><td>0.94</td><td>0.346</td><td></td></tr></table>

<table><thead><tr><th>term</th><th>estimate</th><th>SE</th><th>z value</th><th>Pr(&gt;|z|)</th><th></th></tr></thead><tbody><tr><td>(Intercept)</td><td>-0.0003</td><td>0.0043</td><td>0.07</td><td>0.947</td><td></td></tr><tr><td><strong>temperature</strong></td><td><strong>-0.0046</strong></td><td><strong>0.0018</strong></td><td><strong>2.53</strong></td><td><strong>0.012</strong></td><td><strong>*</strong></td></tr><tr><td>precipitations</td><td>-0.0097</td><td>0.0051</td><td>1.83</td><td>0.068</td><td></td></tr><tr><td><strong>forest</strong></td><td><strong>-0.0002</strong></td><td><strong>0.0001</strong></td><td><strong>2.54</strong></td><td><strong>0.011</strong></td><td><strong>*</strong></td></tr></tbody></table>
<table><tr><td>forest</td><td>4</td><td>93.95</td><td>-178.79</td><td>3.55</td><td>0.06</td></tr><tr><td>forest + precipitation</td><td>5</td><td>94.67</td><td>-177.63</td><td>4.71</td><td>0.04</td></tr><tr><td>precipitation + temperature</td><td>5</td><td>94.65</td><td>-177.58</td><td>4.76</td><td>0.04</td></tr><tr><td>forest + wetland</td><td>5</td><td>94.31</td><td>-176.91</td><td>5.44</td><td>0.03</td></tr><tr><td>(Null)</td><td>3</td><td>91.43</td><td>-176.22</td><td>6.13</td><td>0.02</td></tr><tr><td>temperature</td><td>4</td><td>92.62</td><td>-176.13</td><td>6.22</td><td>0.02</td></tr><tr><td>precipitation + temperature + wetland</td><td>6</td><td>95.13</td><td>-175.79</td><td>6.55</td><td>0.01</td></tr><tr><td>temperature + wetland</td><td>5</td><td>93.55</td><td>-175.39</td><td>6.96</td><td>0.01</td></tr><tr><td>forest + precipitation + wetland</td><td>6</td><td>94.79</td><td>-175.11</td><td>7.23</td><td>0.01</td></tr><tr><td>wetland</td><td>4</td><td>91.83</td><td>-174.55</td><td>7.79</td><td>0.01</td></tr><tr><td>precipitation</td><td>4</td><td>91.69</td><td>-174.27</td><td>8.07</td><td>0.01</td></tr><tr><td>precipitation + wetland</td><td>5</td><td>91.94</td><td>-172.16</td><td>10.18</td><td>0</td></tr></table>

<table><tr><td>wetland</td><td>0.0051</td><td>0.0049</td><td>0.99</td><td>0.321</td></tr></table>

730

731 df: degree of freedom, logLik: log likelihood, ΔAICc: difference in corrected Akaike Information Criterion relative to model with lowest AICs.

732 Highlighted models (grey area) are based on ΔAICc < 2. Highlighted variables (bolded characters) are based on α < 0.05.

733
734

<caption>Table 2. Posterior estimations of demographic parameters (means and 95% CI)</caption>

735

following 4·10⁶ simulations and splits designed under scenario 2 (see Fig. 2). Table

736

reports the effective population size (Ne) of each ecological unit, and the generation time

737

before sampling associated with each out three divergence times (split of the rear edge: t1,

738

split of the admixture zone: t2, and separation of the leading edge from continuous range: t3;

739

see Fig. 2). For each parameter, we estimated estimation accuracy by generating 500 datasets

740

and compile the Relative Median Absolute Error (RMAE) based on their posterior

741

distributions. The corresponding number of years were assumed considering a generation time

742

of three years. See text for details.

743

<table><thead><tr><th>Ecological units</th><th>Ne</th><th>95% CI</th><th>RMAE</th></tr></thead><tbody><tr><td>Rear edge (&lt; 100 m)</td><td>6,100</td><td>2,170 – 10,030</td><td>0.169</td></tr><tr><td>Admixture zone (100-500 m)</td><td>6,390</td><td>2,580 – 10,200</td><td>0.146</td></tr><tr><td>Continuous range (500-1300 m)</td><td>8,730</td><td>2,710 – 14,750</td><td>0.104</td></tr><tr><td>Leading Edge (&gt; 1300 m)</td><td>1,600</td><td>1,014 – 2,186</td><td>0.171</td></tr></tbody></table><table><thead><tr><th></th><th>Generations</th><th>95%CI</th><th>RMAE</th><th>Generation time</th><th>Years</th><th>95%CI</th></tr></thead><tbody><tr><td>t1</td><td>37</td><td>21 – 53</td><td>0.262</td><td>2.5</td><td>92.5</td><td>52.5 – 132.5</td></tr><tr><td>t2</td><td>191</td><td>108 – 274</td><td>0.211</td><td>2.0</td><td>382</td><td>216 – 548</td></tr><tr><td>t3</td><td>748</td><td>536 – 960</td><td>0.120</td><td>2.0</td><td>1,496</td><td>1,072 – 1,920</td></tr></tbody></table>

744
745 **Figure captions**

746 **Fig. 1.** Genetic structuration of 42 natural populations of the common lizard (Z. v.
747 *louislantzi*) gathered into 4 ecological units associated with the recent shift in
748 **distribution of the Western clade.** Individuals (n = 599) were clustered without *a priori* into
749 K = 6 clusters (one colour per K) following a Bayesian approach (top panel). We used the
750 maximal proportion of individuals affiliated to a cluster within a population (as illustrated by
751 the 42 pie charts) to represent an index of population structuration. Although structuration
752 depended on ecological units (mean ± SE on the right panel), the different clusters depended
753 on geographic localisation rather than ecological units (bottom panels).

754

755 **Fig. 2.** **Three demographic scenarios of Pyrenean ascension in the common lizards.** For
756 each ecological unit, the total population size (N) was kept constant. Posterior estimates of
757 $10^6$ simulations for each scenario were compared using logistic regression and scenario 2 was
758 retained with the highest posterior probability (mean ± 95% CI: 0.74 ± 0.06), compared to
759 scenario 1 (0.10 ± 0.04) and scenario 3 (0.16 ± 0.17). The selected scenario associates three
760 divergence times between ecological units across Pyrenean ascension (split of the rear edge:
761 t1, split of the admixture zone: t2, and separation of the leading edge from continuous range:
762 t3).

763

764 **Fig. 3.** **Comparison of population genetics (means ± SE) among ecological units and associated predictions (grey polygons).** (a) Within-populations genetic diversity decreases across altitude much like (b) inbreeding, while (c) genetic differentiation follow a non-linear trend following our predictions.

768
769

770 **Fig. 4.** Population connectivity tested by isolation by distance (IBD) both globally and locally. Significant relationships between genetic distance and geographic distance are symbolized with line types (solid: significant, dashed: non-significant) together with a 2- dimensional kernel density estimation to illustrate population continuity. (a) Globally, genetic distances are associated with geographic distance and continuous, (b,c) IBD was non- significant in isolated populations, and (d,e) IBD was significant but discontinuous in connected populations.

777

778 **Fig. 5.** Variation of demographic and genetic markers with climate conditions and habitat characteristics. The line and surfaces are drawn from predictions of the best linear mixed models selected by the best conditional average (Table 1). The color gradient highlights the values from low (blue) to high range (orange). a) Lizard abundance was negatively correlated solely to air temperature, while b) within-population diversity, c) inbreeding or d) population differentiation were influenced by additive effects of air temperature and forest cover.

784
<figure><img src="image_3.png" /></figure>

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

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

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

Continuous
range

500-1300m

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

Admixture
zone

100-500m

<100m

Scenario 2

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

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

Leading
edge

>1300m

Continuous
range

500-1300m

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

Admixture
zone

100-500m

<100m

Scenario 3

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

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

Leading
edge

>1300m

Continuous
range

500-1300m

Admixture
zone

100-500m

<100m
Within-population diversity (Ar)

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

Inbreeding (F<sub>IS</sub>)

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

Differentiation (F<sub>ST</sub>)

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

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

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

<figure><img src="image_18.png" /></figure>
<figure><img src="image_19.png" /><figcaption>Geographic distance (km)</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>c)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>d)</figcaption></figure>

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