Andréaz Dupoué
email: andreaz.dupoue@gmail.com

Audrey Trochet

Mahaut Murielle Richard

Victoire Sorlin

Michael Guillon

Jules Teulieres-Quillet

Clément Vallé

Cyrielle Rault

Maud Berroneau

Matthieu Berroneau

Mahaut Sorlin

Michaël Guillon

Jules Teulière Quillet

Maud Berroneau

Matthieu Berroneau

Olivier Lourdais

Pauline Blaimont

Romain Bertrand

Gilles Pottier

Olivier Calvez

Olivier Guillaume

Hugo Le Chevalier

Jérémie Souchet

Jean François Le Galliard

Jean Clobert

Fabien Aubret

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

Keywords: South-west France Abundance, climate, colonization, ectotherm, forest cover, gene flow, inbreeding, population decline, structuration

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.

Introduction

Biodiversity is facing a global loss of habitats resulting in altered species and biotic interactions worldwide [START_REF] Ceballos | Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines[END_REF][START_REF] Plotnick | The fossil record of the sixth extinction[END_REF][START_REF] Wiens | Climate-related local extinctions are already widespread among plant and animal species[END_REF]. Causes are multiple, often synergetic, and have various footprints from genes to ecosystems [START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being[END_REF]. As a consequence, many species' distributions are shifting [START_REF] Lenoir | Species better track the shifting isotherms in the oceans than on lands[END_REF][START_REF] Vanderwal | Focus on poleward shifts in species' distribution underestimates the fingerprint of climate change[END_REF] and these range shifts have direct impacts on processes and patterns of population genetics [START_REF] Carvalho | Genes on the edge: A framework to detect genetic diversity imperiled by climate change[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF][START_REF] Templeton | Disrupting evolutionary processes: The effect of habitat fragmentation on collared lizards in the Missouri Ozarks[END_REF].

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 precursor of decline in small populations [START_REF] Frankham | Inbreeding in the wild really does matter[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. 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 [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF][START_REF] Nadeau | Eco-evolution on the edge during climate change[END_REF][START_REF] Vilà-Cabrera | Refining predictions of population decline at species' rear edges[END_REF][START_REF] Waters | Founder takes all: Density-dependent processes structure biodiversity[END_REF]. Studies have challenged these general assumptions since genetic diversity does not necessarily correlate with biogeography [START_REF] Eckert | Genetic variation across species' geographical ranges: The central-marginal hypothesis and beyond[END_REF][START_REF] Pironon | Geographic variation in genetic and demographic performance: new insights from an old biogeographical paradigm[END_REF] nor does it explain population decline [START_REF] Tobler | Genetic attributes of midwife toad (Alytes obstetricans) populations do not correlate with degree of species decline[END_REF]. 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 [START_REF] Galliard | Ecological effects of climate change on European reptiles[END_REF]. These organisms are characterized by a strong thermal dependence on their performance, which may explain their spatial distribution [START_REF] Buckley | Broad-scale ecological implications of ectothermy and endothermy in changing environments[END_REF]. 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 [START_REF] Hewitt | The genetic legacy of the quaternary ice ages[END_REF]. 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 [START_REF] Sunday | Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation[END_REF].

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 [START_REF] Dupoué | Reduction of baseline corticosterone secretion correlates with climate warming and drying across wild lizard populations[END_REF][START_REF] Dupoué | Water availability and temperature induce changes in oxidative status during pregnancy in a viviparous lizard[END_REF][START_REF] Galliard | Ecological effects of climate change on European reptiles[END_REF][START_REF] Rozen-Rechels | When water interacts with temperature: Ecological and evolutionary implications of thermo-hydroregulation in terrestrial ectotherms[END_REF][START_REF] Sinervo | Erosion of lizard diversity by climate change and altered thermal niches[END_REF]. 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;[START_REF] Lourdais | Hydric "costs" of reproduction: pregnancy increases evaporative water loss in the snake Vipera aspis[END_REF].

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 [START_REF] Voituron | To freeze or not to freeze ? An evolutionary prespective on the cold hardiness strategies of overwinrering ectotherms[END_REF], as well as reproductive bimodality with distinct oviparous and viviparous populations across its distribution range [START_REF] Surget-Groba | Multiple origins of viviparity, or reversal from viviparity to oviparity? The European common lizard (Zootoca vivipara, Lacertidae) and the evolution of parity[END_REF]. The viviparous form has colonized the Eurosiberian region from Western Europe to Hokkaido island in Japan being restricted by the polar circle [START_REF] Surget-Groba | Multiple origins of viviparity, or reversal from viviparity to oviparity? The European common lizard (Zootoca vivipara, Lacertidae) and the evolution of parity[END_REF]. 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 [START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. 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), man-made draining of natural habitats (peatbogs and marshes) occurred in the mid-19 th century, which has isolated populations into a few suitable patches of humid areas [boggy islets [START_REF] Berroneau | Lézard vivipare Zootoca vivipara (Lichtenstein, 1823)[END_REF]]). 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 Mountains between 100 and 500 m ASL, where populations are enclaved within forest-isolated peat bogs. This area is located between populations from the rear edges and those of the continuous range, now forming the Admixture zone [sensu [START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]], a lineage mixing area between previously allopatric sub-clades [START_REF] Rius | How important is intraspecific genetic admixture to the success of colonising populations?[END_REF]. 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 [START_REF] Horreo | Phylogeography, evolutionary history and effects of glaciations in a species (Zootoca vivipara) inhabiting multiple biogeographic regions[END_REF]. 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 [START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. Conversely, in the wave front of the expanding margin, we predicted higher segregation and lower genetic diversity if colonisation occurs through founder events [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF][START_REF] Nadeau | Eco-evolution on the edge during climate change[END_REF][START_REF] Waters | Founder takes all: Density-dependent processes structure biodiversity[END_REF]. 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;[START_REF] Lorenzon | Effect of water constraint on growth rate , activity and body temperature of yearling common lizard (Lacerta vivipara)[END_REF]Rozen-Rechels et al., in press;[START_REF] Rutschmann | Climate and habitat interact to shape the thermal reaction norms of breeding phenology across lizard populations[END_REF][START_REF] Rutschmann | Climate dependent heating efficiency in the common lizard[END_REF].

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 [START_REF] Beebee | Buccal swabbing as a source of DNA from squamate reptiles[END_REF], 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 = 0.626) and a similar index provided a reliable proxy of absolute density in another study (Dupoué, Rutschmann, Le Galliard, Clobert et al., 2017). 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). This included the boggy islets in the rear edge (n = 55 lizards among 4 populations), the forest enclaves of the admixture zone (137 lizards among 9 populations), the mixed habitats of the continuous range (242 lizards among 16 populations), and the highland meadows in the leading edge (198 lizards among 13 populations).

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 [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF][START_REF] Karger | High resolution climate data for Europe[END_REF] 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 (Fig. S1A) and annual cumulative precipitation flux corrected by orographic effects (Fig. 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 (Fig. S1C), derived from many studies (http://www.globalforestwatch.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 (Fig. S1D), calculated from pedological and hydrogeological information [START_REF] Berthier | Enveloppes des milieux potentiellement humides de la France métropolitaine -notice d'accompagnement. Programme de modélisation des milieux potentiellement humides de France[END_REF]. This index provides a continuous measure of the hydric environment in each population, which fitted well with our direct observation in the field (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-10ng of template DNA. PCR products were run on an ABI 3730 DNA Analyser (Applied Biosystems) with Genescan-600Liz 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 was tested using GENEPOP 4.7.0 [START_REF] Rousset | GENEPOP'007: A complete re-implementation of the GENEPOP software for Windows and Linux[END_REF] 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 [START_REF] Jones | COLONY: A program for parentage and sibship inference from multilocus genotype data[END_REF] to identify and remove relatives in each population (Table S4) given that they may alter subsequent analyses [START_REF] Rodríguez-Ramilo | The effect of close relatives on unsupervised Bayesian clustering algorithms in population genetic structure analysis[END_REF]. The final dataset included n = 599 individuals.

Statistical analyses

Genetic metrics

We used R software (R Development Core Team 2015) and the DiveRsity R-package [START_REF] Keenan | DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors[END_REF] to estimate genetic diversity in each population including expected heterozygosity (He), observed heterozygosity (Ho), deviation test to Hardy-Weinberg equilibrium (HWE), and allelic richness (Ar). We used the same package to estimate heterozygote deficit within a population (FIS) and genetic differentiation between populations (FST) after 1,000 bootstraps. In all these analyses, we excluded one population (Lagunas, 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 [START_REF] Kamvar | Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality[END_REF] to analyse the molecular variance (AMOVA) and check how variation in genetic diversity differed i) among ecological units, ii) among or within populations, and iii) within individuals. We used the program BOTTLENECK [START_REF] Piry | BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data[END_REF] to check if any populations experienced a recent bottleneck. Shortly after a bottleneck event, heterozygosity increases more rapidly than allelic diversity [START_REF] Piry | BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data[END_REF]. The program tests the heterozygosity excess and compares Ho and He (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 [START_REF] Piry | BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data[END_REF].

Genetic structure

We used the Bayesian clustering analysis of the program STRUCTURE 2.3.3 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] 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 [START_REF] Earl | STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF]. 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 program, illustrating relatively few admixtures between population (alpha < 1) and good model convergence (Fig. S2). We estimated structuration of each population as the highest probability (maximal Q-value) for an individual to be affiliated to a cluster (Fig. 1).

We complemented this analysis by using another method for the examination of population structure. We considered a method based on Discriminant Analysis of Principal Component (DAPC) proposed by [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF] and implemented in the adegenet R package [START_REF] Jombart | adegenet: A R package for the multivariate analysis of genetic markers[END_REF]. This multivariate method generates models for a range of K clusters and calculate Bayesian Information Criterion (BIC) after 10,000 iterations [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF]. Optimal K is chosen following a stepping stones approach, for the minimum K when BIC decreases by a negligible amount (Fig. 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 (DC) between populations defined by [START_REF] Cavalli-Sforza | Phylogenetic analysis models and estimation procedures[END_REF], since these are more likely to retrieve correct relations among populations [START_REF] Takezaki | Genetic distances and reconstruction of phylogenetic trees from microsatellite DNA[END_REF]. 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 [START_REF] Barton | MuMIn: Multi-Model Inference[END_REF] 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 [START_REF] Keenan | DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors[END_REF]]. We ran four independent sets of 15 linear mixed model [R package nlme [START_REF] Pinheiro | nlme: Linear and Nonlinear Mixed Effects Models (R package version 3[END_REF]] comparisons, using either demographic (lizard abundance) or genetic markers (Ar, FIS and FST) 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 [START_REF] Dormann | Collinearity: A review of methods to deal with it and a simulation study evaluating their performance[END_REF]. 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 [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]], and calculated the conditional model average to test for effect size and significance of fixed co-variates. We verified that residuals of all selected models (Table 1) respected normality conditions (Shapiro test: all p > 0.104).

Pyrenean colonisation scenarios

We used Approximate Bayesian Computation (ABC) to test different Pyrenean colonisation scenarios using DIYABC 2.1.0 [START_REF] Cornuet | Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0)[END_REF]. First, we determined if variation in population genetics within and between ecological units coincided to an ancient or a recent Pyrenean colonisation. We followed the same methodology as described in [START_REF] Ursenbacher | Postglacial recolonization in a cold climate specialist in western Europe: patterns of genetic diversity in the adder (Vipera berus) support the central-marginal hypothesis[END_REF] and generated 3•10 6 simulations to compare the posterior probabilities (logistic regression) of three scenarios (Fig. 2) including: a simultaneous split among the four ecological units (scenario 1), an altitudinal progressive colonisation (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 (Fig. 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 > t1 as conditional priors (Fig. 2) and we considered the mean number of alleles, the mean gene diversity, the mean size variance, and, between two samples, the FST, 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 [START_REF] Cornuet | Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0)[END_REF]. We used a generation time of 3 years based on the equation T =  + [s/1(1-s)] [START_REF] Lande | Stochastic population dynamics in ecology and conservation[END_REF] where  is the age at maturity and s is the annual survival rate. Survival and age at maturity are likely to vary along altitude [START_REF] Heulin | Demography of a bimodal reproductive species of lizard (Lacerta vivipara): survival and density characteristics of oviparous populations[END_REF], so we considered adult survival rate of 0.51 (low and intermediate altitudes) or 0.32 (high altitudes) and a recruitment in reproduction at one year old (low and intermediate altitudes) or 2 years old (high altitudes).

Results

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: F3,37 = 3.7, p = 0.020; FIS: F3,37 = 4.6, p = 0.008; FST: F3,37 = 4.2, p = 0.012). Ar was the highest in the admixture zone and FIS were highest in the rear edge and decreased along altitudinal clines (Fig. 3a & 3b). FST was higher in both the rear and leading edge but the lowest in the admixture zone (Fig. 3c). We did not find any evidence for a recent bottleneck (p > 0.216), except in one population (BTM, p = 0.002).

Genetic structure

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

Genetic isolation

Isolation by distance (IBD) was highly significant globally (p = 0.001, Fig. 4a). Locally however, IBD differed among ecological units and according to isolation level. IBD was nonsignificant at both the rear edge (p = 0.276, Fig. 4b) and the admixture zone (p = 0.320, Fig. 4c), where all pairwise geographic distances were below 60 km. Instead, significant IBD was found in both the continuous range (p = 0.006, Fig. 4d) and the leading edge (p = 0.003, Fig. 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, Fig. 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, Fig. 5b). Variation in both FIS and FST were better explained by additive constraints in temperature and forest cover (Table 1, Fig. 5c &5d). That is, inbreeding and differentiation is higher in warm and open habitats (Table 1).

Pyrenean colonisation scenarios

The most likely scenario was the colonization scenario 2 with a posterior probability of 0.86 ± 0.04 (Fig. 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,500 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).

Discussion

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 [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. 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 [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF][START_REF] Nadeau | Eco-evolution on the edge during climate change[END_REF][START_REF] Waters | Founder takes all: Density-dependent processes structure biodiversity[END_REF]. Lower inbreeding levels in these colonized areas might for instance reflect a purge of inbred individuals during recruitment [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF], since this species can optimise mate choice for that purpose [START_REF] Richard | Optimal level of inbreeding in the common lizard[END_REF].

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;[START_REF] Horreo | Population structure of the oviparous South-West European common lizard[END_REF][START_REF] Milá | Multilocus phylogeography of the common lizard Zootoca vivipara at the Ibero-Pyrenean suture zone reveals lowland barriers and high-elevation introgression[END_REF]. 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 (Horreo, Breedveld, et al., 2019;[START_REF] Milá | Multilocus phylogeography of the common lizard Zootoca vivipara at the Ibero-Pyrenean suture zone reveals lowland barriers and high-elevation introgression[END_REF]. 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 [START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF]. 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 [START_REF] Päckert | Genetic admixture despite ecological segregation in a North African sparrow hybrid zone (Aves, Passeriformes, Passer domesticus × Passer hispaniolensis)[END_REF]. Surprisingly however, we found little evidence for bottlenecks (only one population likely due to a local isolation by a road), contrary to [START_REF] Horreo | Population structure of the oviparous South-West European common lizard[END_REF] 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 [START_REF] Horreo | Population structure of the oviparous South-West European common lizard[END_REF]] 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 [START_REF] Horreo | Population structure of the oviparous South-West European common lizard[END_REF]. 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 [START_REF] Cornetti | Genetic and ecological data reveal species boundaries between viviparous and oviparous lizard lineages[END_REF][START_REF] Templeton | Disrupting evolutionary processes: The effect of habitat fragmentation on collared lizards in the Missouri Ozarks[END_REF]. 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;[START_REF] Dupoué | Reduction of baseline corticosterone secretion correlates with climate warming and drying across wild lizard populations[END_REF][START_REF] Dupoué | Water availability and temperature induce changes in oxidative status during pregnancy in a viviparous lizard[END_REF]Dupoué, Rutschmann, Le Galliard, Clobert, et al., 2017;[START_REF] Massot | Climate warming, dispersal inhibition and extinction risk[END_REF]. Furthermore, forest cover systematically exhibited an additive influence on withinpopulation 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 [START_REF] Rutschmann | Climate and habitat interact to shape the thermal reaction norms of breeding phenology across lizard populations[END_REF] 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 [START_REF] Jofré | The role of managed coniferous forest in the conservation of reptiles[END_REF]. 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 [START_REF] Carvalho | Genes on the edge: A framework to detect genetic diversity imperiled by climate change[END_REF][START_REF] Ikeda | Genetically informed ecological niche models improve climate change predictions[END_REF][START_REF] Razgour | Considering adaptive genetic variation in climate change vulnerability assessment reduces species range loss projections[END_REF].

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 [START_REF] Cornuet | Inference on population history and model checking using DNA sequence and microsatellite data with the software DIYABC (v1.0)[END_REF].

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 [START_REF] Davasse | 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[END_REF]. As unravelled by palynology, this area was mostly deforested during the Middle Ages to develop pastoralism [START_REF] Galop | 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[END_REF][START_REF] Zanon | European forest cover during the past 12,000 years: A palynological reconstruction based on modern analogs and remote sensing[END_REF], 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 [START_REF] Gillet | Evidence of finescale genetic structure for the endangered Pyrenean desman (Galemys pyrenaicus) in the French Pyrenees[END_REF].

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 [START_REF] Hampe | Efficient mitigation of founder effects during the establishment of a leading-edge oak population[END_REF][START_REF] Hampe | Conserving biodiversity under climate change: the rear edge matters[END_REF][START_REF] Nadeau | Eco-evolution on the edge during climate change[END_REF][START_REF] Razgour | The shaping of genetic variation in edge-of-range populations under past and future climate change[END_REF][START_REF] Vilà-Cabrera | Refining predictions of population decline at species' rear edges[END_REF][START_REF] Waters | Founder takes all: Density-dependent processes structure biodiversity[END_REF]. 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 [START_REF] Stevens | A comparative analysis of dispersal syndromes in terrestrial and semi-terrestrial animals[END_REF] Highlighted models (grey area) are based on ΔAICc < 2. Highlighted variables (bolded characters) are based on α < 0.05.

Table 2. Posterior estimations of demographic parameters (means and 95% CI)

following 4•10 6 simulations and splits designed under scenario 2 (see Fig. 2). Table reports the effective population size (Ne) of each ecological unit, and the generation time before sampling associated with each out three divergence times (split of the rear edge: t1, split of the admixture zone: t2, and separation of the leading edge from continuous range: t3; see Fig. 2). For each parameter, we estimated estimation accuracy by generating 500 datasets and compile the Relative Median Absolute Error (RMAE) based on their posterior distributions. The corresponding number of years were assumed considering a generation time of three years. See text for details. 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). The selected scenario associates three divergence times between ecological units across Pyrenean ascension (split of the rear edge: t1, split of the admixture zone: t2, and separation of the leading edge from continuous range: t3). q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 0.4 0.5 0.6 0.7 0 100 200 300 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Genetic distance (D C ) Geographic distance (km) a) Global q q q q q q 0.3 0.4 0.5 0.6 0.7 0 20 40 60 b) Rear edge q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 0.4 0.5 0.6 0.7 c) Admixture zone 0 20 40 60 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 0.4 0.5 0.6 0.7 d) Continuous range 0 40 80 120 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 0.4 0.5 0.6 0.7 e) Leading edge 0 40 80 120 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q a)

  , yet our study challenged this generality and suggested fast colonisation capabilities in the common lizard likely facilitated by past habitat connectivity. Together, our results therefore illustrate the importance of integrating genetic diversity, differentiation and colonisation capabilities to accurately characterise and anticipate the potential footprints of global change on wild species and ecosystems. 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. 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.

Figure captions Fig. 1 .

 captions1 Figure captions

Fig. 2 .

 2 Fig. 2. Three demographic scenarios of Pyrenean ascension in the common lizards. For

Fig. 3 .

 3 Fig. 3. Comparison of population genetics (means ± SE) among ecological units and

Fig. 4 .

 4 Fig. 4. Population connectivity tested by isolation by distance (IBD) both globally and

Fig. 5 .

 5 Fig. 5. Variation of demographic and genetic markers with climate conditions and

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)].

 1 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.

	Abundance	model	df	logLik	AICc	ΔAICc	wi	term	estimate	SE	z value	Pr(>|z|)
		temperature	4	-73.9	156.87	0	0.45	(Intercept)	-0.005	0.254	0.02	0.986
		forest + temperature	5	-73.73	159.12	2.25	0.15	temperature	-0.476	0.096	4.80	0.000	***
		temperature + wetland	5	-73.79	159.25	2.38	0.14	forest	-0.002	0.004	0.54	0.593
		precipitation + temperature	5	-73.89	159.44	2.57	0.12	wetland	0.134	0.290	0.45	0.655
		forest + temperature + wetland	6	-73.65	161.71	4.84	0.04	precipitations	0.056	0.309	0.18	0.860
		forest + precipitation + temperature	6	-73.72	161.83	4.96	0.04
		precipitation + temperature + wetland	6	-73.75	161.9	5.02	0.04
		forest + precipitation + temperature + wetland	7	-73.62	164.53	7.66	0.01
		(Null)	3	-79.21	165.06	8.19	0.01
		precipitation	4	-78.75	166.59	9.72	0
		wetland	4	-78.96	167	10.13	0
		forest	4	-79.03	167.14	10.27	0
		precipitation + wetland	5	-78.35	168.36	11.49	0
		forest + precipitation	5	-78.71	169.09	12.22	0
		forest + wetland	5	-78.74	169.16	12.28	0

730

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

731

Acknowledgements

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.

Funding

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

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

Data accessibility

The data supporting this article are available at the Zenodo Repository at http://doi.org/10.5281/zenodo.4266021

Ethics

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 (#2017-s-02) relating to an authorization of capture, marking, transport, detention, use and release of protected lizard species.

Author contributions

AD and AT conceived the study together with JC and FA. AD, MR, MS, JTQ, CV, CR, MauB, MatB, OL and PB captured lizards. MR and MS performed molecular assays. AD, AT