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

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


| INTRODUC TI ON
Biodiversity is facing a global loss of habitats resulting in altered species and biotic interactions worldwide (Ceballos et al., 2017;Plotnick et al., 2016;Wiens, 2016). Causes are multiple, often synergetic, and have various footprints from genes to ecosystems (Pecl et al., 2017).
As a consequence, many species' distributions are shifting (Lenoir et al., 2020;VanDerWal et al., 2013) and these range shifts have direct impacts on processes and patterns of population genetics (Carvalho et al., 2019;Hampe & Petit, 2005;Templeton et al., 2001). In particular, once populations at the rear edge become isolated, the loss of connectivity and immigration progressively leads to an increase in inbreeding through genetic drift, one of the main precursors of decline in small populations (Frankham, 2010;Hampe & Petit, 2005).
In contrast, at the leading edge, newly favourable habitats should be colonized by few individuals resulting in founder events, which implies lower genetic diversity associated with small groups of emigrants and recurrent extinction events (Hampe et al., 2013;Hampe & Petit, 2005;Nadeau & Urban, 2019;Vilà-Cabrera et al., 2019;Waters et al., 2013). Studies have challenged these general assumptions since genetic diversity does not necessarily correlate with biogeography (Eckert et al., 2008;Pironon et al., 2017) nor does it explain population decline (Tobler et al., 2013). This prompts for empirically testing patterns of genetic variation and phylogenetic history along colonization gradients to better understand how population genetic diversity or structure correlate with range shifts of wild populations.
Terrestrial ectotherms are well-suited organisms for this purpose since their ecology and demography tightly depend on climate conditions and landscape structure (Le Galliard et al., 2012). These organisms are characterized by a strong thermal dependence of their performance, which may explain their spatial distribution (Buckley et al., 2012). Under temperate climates, ectotherm distributions at a continental scale have been strongly 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 characterize 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., 2018Le Galliard et al., 2012;Rozen-Rechels et al., 2019;Sinervo et al., 2010;Stier et al., 2017). 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;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 south-west 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 and Petit (2005). This allows for the examination of recolonization dynamics and how multiple demographic trajectories shape genetic variation since the Last Glacial Maximum (LGM). In the rear edge within 100 m above sea level (ASL), manmade draining of natural habitats (peatbogs and marshes) occurred in the mid-19th century, which has isolated populations into a few suitable patches of humid areas [boggy islets (Berroneau, 2014)].
Surrounding environments became inhospitable because of human landscape modification and recent climate warming, thus hindering any emigration from these populations. Similar anthropogenic-driven isolations have occurred in the foothills of the French Pyrenees financially supported by the Australian Research Council (no. DP180100058 to A/ Prof Anne Peters).
Editor: Orly Razgour and genetic trends were better explained by inhospitable (warm and dry) climate conditions and forest cover.

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

K E Y W O R D S
abundance, climate, colonization, ectotherm, forest cover, gene flow, inbreeding, population decline, structuration Mountains between 100 and 500 m ASL, where populations are enclaved within forest-isolated peatbogs. This area is located between populations from the rear edges and those of the continuous range, now forming the Admixture zone [sensu (Hampe & Petit, 2005)], a lineage mixing area between previously allopatric subclades (Rius & Darling, 2014). On the other hand, in its continuous range from 500 to 1,300 m ASL, populations occupy various habitat types mostly composed by river ecosystems and forest borders. Finally, populations are now expanding above 1,300 m ASL and colonizing highland meadows, which represent the leading edge.
We used microsatellite polymorphisms to assess how 42 natural populations across this cline differed in genetic structure, diversity and connectivity. We anticipated that Pyrenean ascension occurred a few times after the LGM, similar to the latitudinal expansion of the viviparous form of this species (Horreo et al., 2018). Additionally, we examined how human-induced environmental threats may additively affect lizard populations by testing the specific influence of recent climate conditions and habitat characteristics on demographic and genetic markers. Because rear edge populations are exposed to increasingly hostile (warm and dry) thermal conditions, they should also experience high levels of differentiation and inbreeding reflecting range contraction (Hampe & Petit, 2005). Conversely, in the wave front of the expanding margin, we predicted higher segregation and lower genetic diversity if colonization 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 Lorenzon et al., 1999;Rozen-Rechels et al., in press;Rutschmann et al., 2016Rutschmann et al., , 2020.

| 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%) had no statistical influence on abundance estimate (p = .626) and a similar index provided a reliable proxy of absolute density in another study . Populations were divided into the four ecological units based on their altitude thresholds as described above (100, 500, and 1,300 m) and the associated landscapes (Table S1)

| Climatic conditions and habitat characteristics in sampled populations
For each population, we extracted two climate variables (temperature and precipitations) of recent climate conditions from CHELSA database, where climatology is averaged between 1981 and 2005 (Karger et al., 2017(Karger et al., , 2020 and is available at a high resolution of 30 Arcsec (~1 km 2 ). For each population, we considered the annual near-surface (2 m) air temperature ( Figure S1a) and annual cumulative precipitation flux corrected by orographic effects ( Figure S1b).
We also considered two descriptors of habitat quality relevant for this species: forest cover and wetland potential. We estimated forest cover at a precision of 30 m from the aboveground live woody biomass density for the year 2000 ( Figure S1c), derived from many studies (http://www.globa lfore stwat ch.org/). Wetland potential (the probability for the habitat to be a humid zone ranging from 0 to 3; 0: dry environments; 1: good probability of wetland area; 2: strong probability of wetland area; 3: very strong probability of wetland area) was obtained at a precision of 50 m ( Figure S1d), calculated from pedological and hydrogeological information (Berthier et al., 2014). This index provides a continuous measure of the hydric environment in each population, which fitted well with our direct observation in the field (A. Dupoué, pers. obs.).

| Genetic data
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 amplified sixteen microsatellite markers (Table S2) in 5 µl reaction volumes with 1 µl of PCR master Kit (MPBiomedicals) and 5-10 ng of template DNA. PCR products were run on an ABI 3730 DNA Analyser (Applied Biosystems) with Genescan-600 Liz size standard. Genotyping was performed with Genemapper 5.0 (Applied Biosystems). We performed an initial examination on the dataset to retain the loci with limited allele dropout (<18% of missing data). Ambiguous genotypes were amplified and sized a second time. Out of 16 initial loci, 5 loci were rejected from further analyses due to an excess of missing data (Table S2). Null alleles and linkage disequilibrium among remaining loci were tested using GENEPOP 4.7.0 (Rousset, 2008) following 10,000 iterations. The 11 loci did not present significant evidence for null alleles (Table S2), nor linkage disequilibrium (Table S3) and were used in the following analyses. We used COLONY software (Jones & Wang, 2010) to identify and remove relatives in each population (Table S4)

| Genetic metrics
We used R software (R Core Team, 2020) and the DiveRsity R package (Keenan et al., 2013) to estimate genetic diversity in each population including expected heterozygosity (He), observed heterozygosity (H o ), deviation test to Hardy-Weinberg equilibrium (HWE) and allelic richness (Ar). We used the same package to estimate heterozygote deficit within a population (F IS ) and genetic differentiation between populations (F ST ) after 1,000 bootstraps.
In all these analyses, we excluded one population (LAG , Table S1) given that the sample size (n = 6 lizards) was too small to obtain reliable estimates of population genetic diversity. We used the poppr R package (Kamvar et al., 2015) to analyse the molecular variance (AMOVA) and check how variation in genetic diversity differed (a) among ecological units, (b) among or within populations, and (c) within individuals. We used the program BOTTLENECK (Piry et al., 1999) to check whether any populations experienced a recent bottleneck. Shortly after a bottleneck event, heterozygosity increases more rapidly than allelic diversity (Piry et al., 1999). The programme tests the heterozygosity excess and compares H o and H e (Wilcoxon test) after 1,000 iterations. We considered a two phase model including a 95% of stepwise mutation model and a variance of 12 as recommended (Piry et al., 1999).

| Genetic structure
We used the Bayesian clustering analysis of the program STRUCTURE 2.3.3 (Pritchard et al., 2000) to determine whether populations could be regrouped without a priori information. We determined how many clusters (K) best fit with our data, based on ΔK method (Earl & VonHoldt, 2012). To do so, we ran 10 independent simulations for each K ranging from 1 to 42, with 500,000 iterations and a burn-in period of 100,000 steps. We adjusted the burn-in period to 200,000 steps in another simulation with K ranging from 1 to 10, to check that the initial burn-in period was sufficiently long to ensure model convergence. We compared two outputs (alpha and likelihood) from the programme, illustrating relatively few admixtures between population (alpha < 1) and good model convergence ( Figure S2). We estimated structuration of each population as the highest probability (maximal Q-value) for an individual to be affiliated to a cluster ( Figure 1).
We complemented this analysis by using another method for the examination of population structure. We considered a ordination method based on Discriminant Analysis of Principal Component (DAPC) proposed by Jombart et al., (2010) and implemented in the adegenet R package (Jombart, 2008). This multivariate method generates models for a range of K clusters and calculate Bayesian information criterion (BIC) after 10,000 iterations (Jombart et al., 2010). Optimal K is chosen following a stepping stones approach, for the minimum K when BIC decreases by a negligible amount ( Figure S3).

| Genetic connectivity
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 C ) 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.

| Effect of environmental conditions on demography and genetic diversity
We used AICc based model selection and averaging from the MuMIn 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 IS and F ST ) 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 covariates. We verified that residuals of all selected models (

| Pyrenean colonization scenarios
We used Approximate Bayesian Computation (ABC) to test different Pyrenean colonization scenarios using DIYABC 2.1.0 (Cornuet et al., 2010). First, we determined if variation in population genetics within and between ecological units coincided to an ancient or a recent Pyrenean colonization. We followed the same methodology as described in (Ursenbacher et al., 2015) and generated 3.10 6 simulations to compare the posterior probabilities (logistic regression) of three scenarios (Figure 2) including: a simultaneous split among the four ecological units (scenario 1), an altitudinal progressive colonization (scenario 2) and a simultaneous split in the Pyrenees regions differing from the rear edge (scenario 3). We restricted the number of tested scenarios by ordering ecological units based on their altitudinal class ( Figure 2). However, we did not test for scenarios with higher proximity between the rear and the leading edge compared to others units since it has no biological meaning. Once the most likely scenario was identified, we further ran 4.10 6 iterations to estimate the parameters (effective population size: N, and divergence time: t). We defined N1 < N2 < N3 > N4 and t3 > t2> F I G U R E 1 Genetic structuration of 42 natural populations of the common lizard (Z. v. louislantzi) gathered into four ecological units (typical habitats illustrated on the top pictures) associated with the recent shift in distribution of the Western clade. Individuals (n = 599) were clustered without a priori into K = 6 clusters (one colour per K) following a Bayesian approach (top chart). We used the maximal proportion of individuals affiliated to a cluster within a population (as illustrated by the 42 pie charts) to represent an index of population structuration. Although structuration depended on ecological units (mean ± SE on the right panel), the different clusters depended on geographic localization rather than ecological units (bottom panels) t1 as conditional priors (Figure 2), and we considered the mean number of alleles, the mean gene diversity, the mean size variance, and, between two samples, the F ST , the mean index of classification, the shared allele distance, and the (dµ) 2 distance as summary statistics. We checked the performance of each ABC analysis (scenario comparison and parameters estimation) by simulating 500 pseudo-observed datasets, comparing their posterior distribution and compiling relative error (Cornuet et al., 2010). We used a generation time of 3 years based on the equation T = α + [s/1(1-s)] (Lande et al., 2003) where α is the age at maturity and s is the an- and a recruitment in reproduction at 1 year old (low and intermediate altitudes) or 2 years old (high altitudes).

| Genetic diversity
Analyses of molecular variance showed that genetic diversity varied mostly within individuals (86.8%), but a substantial and significant part of variation was also explained within populations (4.1%) and between populations (7.4%) followed by a minor variation between ecological units (1.7%, see Table S5). Within-population diversity and population differentiation varied among ecological units (Ar: F 3,37 = 3.7, p = .020; F IS : F 3,37 = 4.6, p = .008; F ST : F 3,37 = 4.2, p = .012). Ar was the highest in the admixture zone, and F IS was highest in the rear edge and decreased along altitudinal clines (Figure 3a,b). F ST was higher in both the rear and leading edge but the lowest in the admixture zone (Figure 3c). We did not find any evidence for a recent bottleneck (p > .216), except in one population (BTM, p = .002).

F I G U R E 2
Three demographic scenarios of Pyrenean ascension in the common lizards. For each ecological unit, the total population size (N) was kept constant. Posterior estimates of 10 6 simulations for each scenario were compared using logistic regression, and scenario 2 was retained with the highest posterior probability (mean ± 95% CI: 0.74 ± 0.06), compared to scenario 1 (0.10 ± 0.04) and scenario 3 (0.16 ± 0.17

| Genetic structure
Both the Bayesian spatial modelling of genetic population structure and the ordination method gathered populations into K = 6 clusters (Table S6 & Figure S2). Clustering was relatively heterogeneous between and within ecological units, and rather dependent on population locality (Figure 1).
Yet, structuration (i.e., the probability for an individual to be affiliated to a cluster) differed significantly between ecological units (F 3,38 = 3.5, p = .024) since there was a lower structuration in the admixture zone compared to the 3 other ecological units (pairwise comparisons, all p < .032).

| Genetic isolation
Isolation by distance (IBD) was highly significant globally (p = .001, Figure 4a). Locally, however, IBD differed among ecological units and according to isolation level. IBD was non-significant at both the rear edge (p = .276, Figure 4b) and the admixture zone (p = .320, Figure 4c), where all pairwise geographic distances were below 60 km. Instead, significant IBD was found in both the continuous range (p = .006, Figure 4d) and the leading edge (p = .003, Figure 4e) where many pairwise geographic distances were above 60 km.

| Effect of past and present environmental conditions on demography and genetic diversity
Spatial variation in lizard abundance was correlated solely to air temperature (Table 1, Figure 5a), with the most scattered populations associated to warm environments. Ar was the lowest under warm and dry climates, and in open habitats (Table 1, Figure 5b).

F I G U R E 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 Variation in both F IS and F ST was better explained by additive constraints in temperature and forest cover (Table 1, Figure 5c,d).
That is, inbreeding and differentiation is higher in warm and open habitats (Table 1).

| Pyrenean colonization scenarios
The most likely scenario was the colonization scenario 2 with a posterior probability of 0.86 ± 0.04 (Figure 2). We compared logistic distribution in 500 simulated datasets to estimate a global posterior predictive error of 0.17, a type I error (probability of rejecting scenario 2 when it is the correct scenario) of 0.33 and a type II error (probability of choosing scenario 2 when it is an incorrect scenario) of 0.22. Effective population size considerably varied between ecological units ranging from 1,600 in the leading edge to 9,000 reproductive individuals in the continuous range ( Table 2). The split between the rear edge and the other populations occurred 750 generations ago, thereby placing the start of the Pyrenean ascension to around 2,000 years ago ( Table 2). The divergence in the other ecological units occurred very recently in the last millennium, with evidence of high elevation colonization between 100 and 700 years before present (Table 2).

| D ISCUSS I ON
Genetic diversity was relatively high among 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 F I G U R E 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 colour gradient highlights the values from low (blue) to high range (orange). (a) Lizard abundance was negatively correlated solely with air temperature, while (b) withinpopulation diversity, (c) inbreeding or (d) population differentiation were influenced by additive effects of air temperature and forest cover  (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 optimize mate choice for that purpose (Richard et al., 2009).
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, 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 among high altitude populations following secondary contacts associated with range shifts in warmer climates 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., 2018Dupoué et al., , 2020Massot 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 TA B L E 2 Posterior estimations of demographic parameters (means and 95% CI) following 4.10 6 simulations and splits designed under scenario 2 (see Figure 2) 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 colonization 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 1,850. 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 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 of extinction. Ground-dwelling ectotherms usually display relatively low dispersal rates (Stevens et al., 2014), yet our study challenged this generality and suggested fast colonization capabilities in the common lizard likely facilitated by past habitat connectivity. Together, our results therefore illustrate the importance of integrating genetic diversity, differentiation and colonization capabilities to accurately characterize and anticipate the potential footprints of global change on wild species and ecosystems.

ACK N OWLED G EM ENTS
We are particularly grateful to the associations "Nature en Occitanie" and "Cistude Nature" who provided the database that allowed us to localize the great majority of lizard populations, and we thank the "Parc National des Pyrenees" for authorizing us to sample lizards. We also thank Dr Sylvain Ursenbacher for advice on genetic analyses.

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

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13202.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data supporting this article are available at the Zenodo