Population decline at distribution margins: Assessing extinction risk in the last glacial relictual but still functional metapopulation of a European butterfly

To determine the interplay between climate and land use changes in driving population dynamics in a butterfly species, Coenonympha hero, at the southern limit of its distribution.


| INTRODUC TI ON
Many species are declining in Europe as the result of climate warming and land use changes. The conversion of natural habitats by human activities is the threatening factor influencing the largest number of species, while climate change is expected to have a growing role in the next decades. Climate change affects the distribution and abundance of species, with increases in abundance in the cooler parts of species' ranges and declines in warmer parts (Fox et al., 2014;Lehikoinen et al., 2013;Parmesan et al., 1999). Changing temperatures can have a particularly marked effect on butterflies, with species at the edge of their distribution showing the most dramatic range shifts (Menéndez et al., 2006;Parmesan et al., 1999;Warren et al., 2001).
The conjunction of climate warming and land use changes is likely to particularly affect cold-adapted boreal species at the southern margin of their distribution and habitat-specialist species that occupy semi-open habitats such as forest glades. Populations in these temporary habitat patches are often small and geographically isolated from each other and thus may strongly suffer under demographic and stochastic environmental fluctuations, making them particularly prone to local extinctions (Melbourne & Hastings, 2008). Furthermore, such isolated populations may express inbreeding depression causing reduced fitness and lower adaptive potential to environmental changes (Habel & Schmitt, 2018;Nonaka et al., 2019). These species persist over the landscape through a dynamic process of extinction/recolonization of local habitat patches with a metapopulation functioning (Hanski & Thomas, 1994). Dispersal is a central process in metapopulation dynamics, yet the effective connectivity between favourable patches is not easy to assess using direct methods such as mark-release-recapture that do not account for (1) rare long-distance dispersal events and (2) effective dispersal, that is the actual success in establishment after dispersion. Genetic methods can give better estimates of effective dispersal than direct methods with less effort. The large number of single nucleotide polymorphisms (SNPs) available from non-lethal sampling together with the development of computationally efficient inference methods such as approximate Bayesian computation (ABC) provides unprecedented amount of information about past and present population effective size, population fragmentation history and contemporary effective gene flow, even for endangered species persisting in small and fragmented populations at their range limits.
Although the recent decade has seen an explosion of species distribution modelling (SDM) based on climate variables to predict future species ranges (reviewed in Zurell et al., 2020), there is increased concern that these global climate predictors might not capture the relevant habitat conditions encountered by individuals, because both topography and land cover can strongly influence local microclimate and population persistence, especially at range margins (Suggitt et al., 2011). For example, there is increasing evidence that climate warming is promoting altitudinal range shifts, which could drive many populations to the extinction in the future (Rödder et al., 2021), and that habitat fragmentation may present a major impediment to species poleward range shifts (Fourcade et al., 2021). Demographic genetic inferences can help improving the performance of SDM to capture past and current habitat suitability of species at their distribution margins, by choosing the environmental predictors best accounting for the population size fluctuations over time and space at fine-spatial scale. Here, we show the complementarity between demographic genetic inferences and species distribution models based on multiple predictors to understand past and current population functioning at species distribution margin, key information for conservation management.
The Scarce Heath Coenonympha hero (Linnaeus, 1761) (Lepidoptera, Nymphalidae) is a butterfly species occurring very locally from eastern France across central Europe and southern Scandinavia to temperate Asia. It is among the most threatened butterfly in Europe and has already become extinct in the Czech suitable habitat patches (semi-open areas) maintaining a dynamical non-equilibrium population despite temporary local extinctions (Cassel-Lundhagen & Sjögren-Gulve, 2007), with large populations acting as sources for small satellite populations (Cassel-Lundhagen et al., 2008). Based on current pattern of genetic variation, we inferred the demographic history of the French population, using ABC to test for alternative demographic scenarios, involving ancient (e.g. several thousand years) or recent (historical) population decline and fragmentation. Based on European occurrences available for C. hero, we modelled the potential current suitable climatic conditions for the species and projected its putative distribution into the Mid-Holocene (MID, 6 kya), Last Glacial Maximum (LGM, 20 kya) and Last Interglacial (LIG, in order to link the demographic history of the species with broad-scale climatic changes. Because local conditions are likely to be important for this low mobile species, we tested whether adding more local environmental factors (global aridity index, altitude, land cover) improved model performance. We then focused on France to construct a habitat model taking into account occurrence status (still present, no-longer present = extinction) and their dates of report in order to evaluate which landscape features are the most important for the species persistence and to guide land use management decisions in a conservation perspective.

| Species distribution models
In order to determine the best environmental predictors of C. hero current distribution in France and Europe and identify changes in climatically suitable areas over time in Europe, we constructed a series of SDMs at global (Europe) and regional (France) scales (Table 1), based on geo-referenced occurrences available for C. hero in the Global Biodiversity Information Facility database (GBIF, https:// www.gbif.org/). We retained field observations recorded after 1900 and records from distribution atlas (unknown date; Kudrna et al., 2011) all with coordinate precision <5 km (Appendix S1) and added 31 new occurrences for the present study samples, resulting in 4615 occurrences in Europe (Table S1). Sampling biases due to unequal sampling effort (e.g. monitoring projects, sampling interests, accessibility) can affect the accuracy of SDM predictions by overrepresenting the ecological importance of the more sampled region (Fourcade et al., 2014). To reduce this possible bias, we randomly selected one occurrence when several points fell within the same raster cell of 1 km 2 (Appendix S1), resulting in 599 occurrences distributed in 12 countries ( Figure 1; Table S1). Depending on the model, environmental predictors included climatic, topographic and land use variables. We collected or calculated a total of 31 environmental variables (Table S2.1: Appendix S2) and retained non-collinear variables ( Figures S2.1 and S2.2). We first constructed a climate-only model (SDM1) using six BIOCLIM variables relevant with respect to the ecology of the species. Indeed, butterflies are highly sensitive to temperature fluctuations, and the main period of activity for C. hero (larval development and imago reproduction) occurs during the warmest months. We thus selected the following: mean annual temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), mean temperature of warmest quarter (BIO8), annual precipitation (BIO12) and precipitation seasonality (BIO15) from WorldClim v2 (Fick & Hijmans, 2017). To determine current and past climatically suitable areas for C. hero in Europe and assess species distribution changes over time, we projected this model in the current , Mid-Holocene (MID), Last Glacial Maximum (LGM) and Last Inter-Glacial (LIG) periods. To determine whether a more detailed characterization of the physical environment improved model performance, we removed variables of low importance (BIO12) and added the Global Aridity Index (GAI) from the CGIAR-CSI v2 database (Trabucco & Zomer, 2019) and elevation (DEM) (European Environment Agency) to construct SDM2. Two other climate-related variables, the growing degree days (GDD) and soil moisture (SMC), were not included because they were strongly correlated to BIO1 and GAI, respectively (Appendix S2).
We then constructed a habitat-only model using the CORINE Land Cover (CLC) map of Europe level 3 (European Environment Agency) to characterize current land use (CLC2018). We reclassified the original 44 classes into 9 levels (Table S2.2), including artificial (urban) and natural (crops, bocage, grass, shrub, forest, water) surfaces, and calculated the percentage of each natural habitat type per km 2 . To represent anthropization, we used the human footprint (HF) index (Venter et al., 2018) from the NASA Socioeconomic Data and Applications Center (SEDAC). We also included in this model the net primary productivity (NPP) as a proxy for global resource availability (Imhoff et al., 2004). Finally, we combined climatic, topographic and habitat predictors retained in SDM2 and SDM3 to construct SDM4 (Table 1; Appendix S2). For these four models, we used 599 presences and three datasets of 10,000 random pseudo-absences were selected using a surface range envelope model.
Among the 599 presences, some are historical records, where the species has not been reported since 1980 (Table S1). Since land use has dramatically changed in Europe during the last decades, considering extinct localities as presence data can potentially skew distribution models. While we have access to detailed information about the status (year of last record) for occurrences in France, distinguishing the areas where the species is extinct and the areas not surveyed since the 1980s for all European countries is difficult ( Figure S2.3). We thus focused on France to build a fifth model (SDM5) based on 135 occurrences: 41 presences, 44 known extinctions (true absences) and 50 absences randomly sampled in areas where the species is absent without occurrences records (Figure 1; Appendix S2) to refine potentially suitable habitats for the species in France. For this model, we used the same climatic and topographic variables as SDM4 but the local habitat was characterized at the time of occurrence report using the last available CLC layer (CLC1992, CLC2000, CLC2006, CLC2012 and CLC2018) ( data  and species occurrence data (records between 1900 and 1970) on the mapped potential distributions (Appendix S2). SDM results including or not these 1900-1970 records are not significantly different so we included them in SDM1 to SDM4 models because excluding them would lead to underestimate the width of the species climatic niche and erroneous past (LIG, LGM, Mid-Holocene) distribution projections.
For each analysis, an ensemble of five statistical models was used, including GLM, GAM, GBM, MAXENT and RF (Araújo & New, 2007;Marmion et al., 2009;Appendix S2). Models were calibrated using 70% of observations randomly sampled from the initial data and evaluated against the remaining 30% data using the true skill statistic (TSS, Allouche et al., 2006) and the area under the curve (ROC, Swets, 1988

| Ordination-based method
In addition to SDM, we performed principal component analyses (PCA) using the ade4 R package (Chessel et al., 2004). First, we com-   Figure 1). Coenonympha hero is mainly found in the Doubs department, in the Plateau du Doubs (DP) and the nearby wetland (DW), and in Jura department (JU).
Some populations have recently been reported at the limit between Doubs and Jura in a centre position (CE). The species was thought to be extinct in the Saône-et-Loire (SAO) department but was recently rediscovered. We collected 149 specimens from 31 localities in the five cited areas with three to six individuals per locality ( Figure 1; Table S2). Individuals were captured and immediately released after one leg removal (non-lethal sampling). SNP genotyping-Among a total of 82 million reads obtained, we trimmed adapters using the BBduk tool of the BBmap package (Bushnell, 2014), removed reads with length <120 nucleotides and cut all reads to this value. We retained 91% high-quality (HQ) reads after removing reads of low quality or with uncalled bases (Table   S3). We used the denovo_map.pl pipeline in STACKS v2 (Rochette et al., 2019) to reconstruct loci and call genotypes, using a maximum of 7 mismatches to merge two stacks into a polymorphic locus (-M; ustacks function). Highly repetitive stacks were dropped using the "Removal" (-r) option. A catalog of the loci from all the individuals was built, with a maximum of 10 mismatches for merging two individual loci (-n; cstacks function). Loci within each individual were searched against the catalog (sstacks function). A SNP data set was produced with the genotype of each individual for every polymorphic position (populations function). We removed called genotypes from per-sample read depth <5. We required a locus to have <20% missing data, and individuals with >25% missing data were discarded

| Population structure and diversity
Population genetic structure and diversity was investigated using the 817 SNP data set. We first performed a PCA in the pcadapt R package (Luu et al., 2017). We then used ADMIXTURE (Alexander et al., 2009) for inferring admixture coefficients for all individuals. The optimal number of genetic clusters, K (range 1-10), was determined based on a 10-fold cross-validation. Genetic diversity was estimated using observed heterozygosity (H o ), population genetic diversity (=expected heterozygosity, H e ), allelic richness (A R ), inbreeding coefficient (F IS ) and population-specific F ST using the hierfstat R package (Goudet, 2005). We also calculated pairwise relatedness coefficients among individuals (Wang, 2002) (Goslee & Urban, 2007).
Spatial patterns of genetic variation among populations were further investigated by computing genetic dissimilarities on a spatial grid of 500 demes among the three French departments using the EEMS programme (Petkova et al., 2016) ( Figure S3.4). We ran three independent analyses with a burn-in of 2,000,000 and MCMC length of 10,000,000 with a thinning interval of 10,000. The results of the three analyses were combined to produce maps for within-demes (diversity rates) and between-demes (migration rates) estimates using the rEEMSplots R package (Petkova et al., 2016) ( Figure S3.5).

| Demographic inferences
To reconstruct the demographic history of C. hero, we used two methods: the stairway plot method v2 based on the site frequency spectrum (SFS) (Liu & Fu, 2015) that reconstructs continuous changes in effective population size over time and the approximate Bayesian computation approach (ABC) implemented in DIYABC v2 (Cornuet et al., 2014) that allows to compare alternative scenarios and to estimate past and current population sizes and discrete times for size changes. Based on population genetic structure results, we defined three populations corresponding to three geographical regions (Appendix S4).
For the stairway plot method, the SFS of each population (DPW, JU, CE) was computed from imputed genotypes of full SNP data sets (2373 SNPs). Missing data were imputed via matrix completion (Chi et al., 2013) from allele frequencies and ancestry coefficients considering three genetic groups. A total sequence length of 254 kb was used, including monomorphic and polymorphic sites among polymorphic loci, as well as monomorphic loci (817 polymorphic loci among 2118 loci) with a mutation rate of 2 × 10 −9 (Appendix S4).
Effective population sizes through time for each population were estimated using the median and 95% confidence intervals of 100 bootstrap replicates.
For the ABC method, we used the 817 SNP data set. We first compared four alternative scenarios of population demographic history: constant population size (Scenario 1), recent decline (Scenario 2), old decline (Scenario 3) and two declines (Scenario 4) ( Figure S4.3). After identifying the best scenario for each population separately, we determined the divergence time between the three populations by comparing three alternative scenarios: divergence before old decline (Scenario 1), divergence due to old decline (Scenario 2) or divergence after recent decline (Scenario 3) ( Figure S4.5). All priors were set to uniform distribution. Prior intervals were defined based on the stairway plot results for ABC reconstruction of single-population demographic history, and using the prior distribution of parameters from ABC reconstruction of single-population demographic history for divergence scenarios (Appendix S4). For summary statistics, we used the within-population genetic diversity (H e ) for single-population demographic history and added the genetic differentiation (F ST ) between populations for the divergence scenarios.
We generated 1,000,000 data sets for each scenario and compared alternative scenarios by calculating their posterior probabilities using a logistic regression on the 1% simulated data closest to the observed data (Cornuet et al., 2010). The goodness of fit of the best scenario was assessed using a PCA, by checking the accordance between observed and simulated data sets from prior distributions.
We evaluated confidence in scenario choice by generating 100 newly simulated data sets from priors to compute type I and type II errors (Cornuet et al., 2010). Once the most likely scenario was identified, the posterior distributions of parameters were estimated by computing a local linear regression on the 1% of the simulated data closest to observed data set, after applying a logit transformation to the parameters value (Cornuet et al., 2010).

| Current and past distribution in Western Europe
The projected potential distribution in Europe of C. hero based on current BIOCLIM data (SDM1) predicts a high suitability around the Baltic Sea (between latitudes 52°N and 60°N) with continu-   Figure 3c).

| Current distribution in France
The comparison of environmental characteristics experienced by  (Table   S1; Figure S2.3) all correspond to known areas of extinction. It is worth noting that these occurrences are located in the warmer climatic conditions (high values for BIO1; Figure 4a) where human pressure is high  We thus performed a French model (SDM5) using this knowledge on the presence and absence (extinctions) of the species (Figure 1). This model revealed that the potential geographical distribution of C. hero in France is mainly reduced to the Jura massif. In accordance with the PCA on climatic variation, GAI is the best environmental predictor of C.
hero distribution in France (69%; Figure 4d).  Investigating the determinants of spatial population structure (distance, barriers), we found significant IBD among all the populations (Mantel, R = .586, p < .001), and similar trends within each region although linear regressions were not significant for JU and CE ( Figure S3.3). We further tested whether populations were more genetically distinct than expected under IBD (i.e. higher genetic drift in each region than gene flow between regions). Based on EEMS migration surfaces, we found genetic differences exceeding IBD suggesting barriers to gene flow between regions ( Figure 6).

| Demographic inferences
Based on population genetic results, we inferred the demographic history of three populations corresponding to three geographical regions sampled (DPW: DP-DW, CE and JU; Appendix S4). In the three regions, there is a deficit in rare alleles compared with neutral expectations, which is suggestive of population decline ( Figure S4.1).
In order to reconstruct the demographic history of C. hero populations in France, we first used the stairway plot method to reconstruct continuous fluctuations in effective population size over time.
We then used the stairway plot results to design alternative scenario topologies and define prior intervals for ABC analysis to estimate the most likely effective population size and population decline timings (Appendix S4). The two methods are complementary, and each can be used to inform the other.
The three populations show similar demographic histories ( Figure 7). For the stairway plot method, we found two events of population decline for DPW and CE while JU population show a strong old population decline followed by a progressive decline over time. For the ABC analysis, we found Scenario 4 as the best scenario: two demographic declines ( Figure S4.3 and Table S4.1), which is congruent with the stairway plot results (Figure 7). Although the posterior probabilities of best scenarios range between 0.51 and 0.68, 95% confidence intervals of posterior probabilities do not overlap with intervals of other scenarios (Table S4.1). However, except for CE, type I errors were rather high (~0.30), suggesting that Scenarios 2 and 3 (one decline) are difficult to tell apart from Scenario 4 (two distinct declines) as also suggested by type II errors (Table S4. In accordance with genetic diversity indices, estimations of effective population size support a south-north gradient with current effective population size estimates twice higher in CE as compared to JU, and five times higher in DPW for both stairway plot and ABC methods (Figure 7; Table S4.2). Timings for old and recent population declines are similar for the three populations (Figure 7).

F I G U R E 7
Comparison of changes in effective population sizes over time obtained by two different methods. The stairway plot method is based on the SFS and requires to fix the spontaneous mutation rate (here 2.0 × 10 −9 ; Appendix S4) and reconstructs continuous fluctuations over time. The ABC method does not require information about the mutation rate but reconstructs discrete events of population decline and requires defining priors. Each population was analysed separately Ancestral effective population sizes: reported for European species include Mediterranean refuges (the Iberian, the Apennine and the Balkan peninsulas) and Central European refuges (Schmitt & Varga, 2012;Taberlet et al., 1998;Wendt et al., 2021).

For boreal European species, extra-Europe
Siberian refuges were also proposed (Schmitt & Varga, 2012). While most of the literature relies on mtDNA variation and haplotype geographical distribution to infer putative Pleistocene refuges (e.g. Dapporto et al., 2019;Ehl et al., 2021), with all the biases inherent to the particular transmission mode of mtDNA (Després, 2019), demographic inference based on multiple nuclear markers allows to go a step further in reconstructing in detail the demographic history of population expansion after the last glacial event. Combining demographic reconstructions and palaeoclimatic SDM further facilitates inferring the location of putative refuges, by identifying stable climatically suitable areas allowing populations to persist during glacial/interglacial oscillations (Ashcroft, 2010;Stewart et al., 2010).
SDM fitted with palaeoclimatic variables (Mid-Holocene, LGM, LIG) revealed that climatically suitable areas for C. hero during the LIG were mostly found in Western Europe from Spain to Norway across a large latitudinal range (Figure 2). During the LGM, the foothill of the Jura massif constituted one of the few identified climatically suitable areas together with the Balkans, but this later area was unsuitable during the LIG. This points the Jura region as the only climatically highly suitable area during the LGM, but also before, and after, during the mid-Holocene and current climates. SDM also provided insights into the factors shaping population demographic history. The two demographic inference approaches were congruent and identified a population decline in two steps: a first old decline before the LGM and a second more recent decline, which resulted in the fragmentation of the Jura population into three subpopulations (Figures 7 and 8). The first demographic decline corresponds to the transition from a warm and wet climate during the LIG to cold and dry climate with ice-sheet covering most northern Europe at LGM, supporting that C. hero persisted during the LGM in the Jura. Current populations in the Jura region being at the southernmost margin of C. hero distribution range, they constitute glacial relict populations.
During post-glacial warming, the ice-free area increased towards North where C. hero is currently found. Yet, despite evidence for a geographical northward expansion, we found no evidence for a population size expansion after the LGM (Figure 7). The climate during mid-Holocene was even less suitable than current climate, with only a few highly suitable patchy areas in north-western Europe ( Figure 2). This mid-Holocene species-hostile picture despite evidence for successful species colonization northwards was possibly buffered by the presence of large deciduous forests that offer pockets of cooler microenvironment during hot and dry summers for climate-sensitive species (Bladon et al., 2020;Suggitt et al., 2011).
Indeed, C. hero is currently mainly found in forested areas throughout Europe, especially in central and eastern Europe although it is an open-canopy butterfly (Figure 3). In accordance with a unique refuge with a limited population size, only a single mitochondrial haplotype (CO1 608 bp, BOLD) is found throughout Europe (Appendix S5), which is much lower than haplotype diversity generally observed in other European butterflies (Cassel-Lundhagen et al., 2020;Després et al., 2019;Dincă et al., 2021).
Using ecological modelling and genetic demographic inferences, we found evidence for a unique Western glacial refuge for C. hero in Europe. Our finding is therefore a contribution to the ever-growing literature showing that many species, including butterflies, survived the LGM in Western Europe outside the Mediterranean peninsulas and that boreal species did not necessarily had a far-east Siberian refuge Minter et al., 2020). A second larger LGM refuge outside Europe for C. hero is pointed out by the higher mtDNA haplotype diversity found in eastern Asia, but this refuge did not contribute to Europe post-glacial recolonization (Appendix S5).  Table 2) suggesting allelic loss through genetic drift. Genetic drift could result from a lack of connectivity between these sites and the core population due to their north-eastern isolation. Indeed, the core population is >15 km, far more than the dispersal capacity of this small butterfly which has been estimated by mark-releaserecapture (MRR) to a few hundred metres (Cassel-Lundhagen & Sjögren-Gulve, 2007). It could also potentially be due to low local habitat quality and small local population, as suggested by SDM4

| Metapopulation functioning
with a lower proportion of high-quality habitats in the southern regions ( Figure 4). The lowest genetic diversity was found in the two southern clusters, CE and JU, and especially in one of the southernmost localities (CER, H e = 0.222). This is unexpected based on main phylogeographic paradigms, which hypothesize lower genetic diversity to the north (e.g. Dincă et al., 2021). The unexpected northern richness and southern lower genetic diversity could be due to rather recent (thousands of years) factors (increased aridity, habitat fragmentation) than to the northward expansion that occurred at the end of the Pleistocene.
Genetic diversity (H e ) within each sampled site ranged from 0.222 to 0.282 which is comparable to the diversity observed in other Coenonympha European butterfly populations using the same ddRAD procedure (Capblancq et al., 2015;Després et al., 2019), and nearly twice that of a currently expanding but young cropfeeding butterfly, Pieris rapae (Ryan et al., 2019). Most sites were at demographic equilibrium (F IS non-significantly different from zero; Table 2), which suggests, together with the weak genetic structure observed at the regional scale, an overall good connectivity between sites and a dynamical metapopulation functioning. This was unex- not related, suggesting that the founding event from the core population is not very recent (more than 3 generations ago) and that the local population is not restricted to the descendants of a single founding female. Furthermore, the overall genetic diversity in this site is not lower than in other more connected localities, yet the high population-specific F ST is a signature of founding effect and further drift due to limited connectivity with the source population.

| Climate and land use constrain C. hero persistence
As all other Satyrinae, C. hero larvae feed on various types of widespread grass and are not specialized on a particular host plant (Tiitsaar et al., 2016). Neither do they depend on a very particular biotope such as two highly endangered wetland Coenonympha sp. only found in low altitude alkaline marshes (C. oedippus), or in cold raised bogs (C. tullia). These two habitat-specialist species are particularly affected by land drainage and are threatened throughout Europe (van Swaay et al., 2006). Depending on the country in Europe, C.
hero can be found in more of less forested steppes, meadows and grasslands, and in semi-open deciduous or mixed forests, fens or transition mires (Sielezniew & Nowicki, 2017;van Swaay et al., 2006 Figure 4). Furthermore, records from Switzerland, Netherlands and Germany fall into this same ecological space, in accordance with the lack of recent observations of the species and suspected extinction in these countries (Kudrna et al., 2011;Lütolf et al., 2006;van Swaay & Warren, 1999). This local extinction pattern supports range contractions at the warm edges in response to climate and habitat changes (Cahill et al., 2014). The mean annual temperature in the Jura region ranging between 7°C and 10°C, these glacial relict populations are at major extinction risk given the unprecedented rate of contemporary climate change (Bennett et al., 2021). Interestingly, despite evidence for a dramatic recent decline especially in the southernmost populations, they retain high genetic diversity, meaning the whole metapopulation might have some evolutionary potential to respond to environmental change. This adaptive potential of C. hero French populations is suggested by their marginal position in the environmental space of the whole species in Europe ( Figure 4) and could potentially contribute to evolutionary rescue northern European populations that appear to be already suffering from genetic erosion (Cassel & Tammaru, 2003) and inbreeding depression (Cassel et al., 2001).
In conjunction with warming, the reduction of suitable areas due to anthropization, as observed for C. hero in France, can have dramatic impacts on the most isolated populations, resulting in between-population genetic differentiation and withinpopulation genetic erosion (Frankham et al., 2014;Hampe & Petit, 2005). Some 2000 years ago, genetic inferences identify a second population decline of similar (Doubs) or much stronger (Jura) amplitude as the first decline (Figure 7), although climatic projections suggest more suitable current conditions especially in northern Europe than during Holocene (Figure 2). This second decline is also supported by records of local extinctions during the last century in the most populated European countries (high human footprint index; Figure 4), suggesting that increasing anthropization and land use changes prevented C. hero population expansion and precipitated its decline. Indeed, human population in Europe has been multiplied by more than 1000 ( Figure S4.7: Appendix S4) during the more recent C. hero population decline. Despite these threats on C. hero, we found no evidence for increased consanguinity nor genetic erosion within sampled sites, which appear to maintain good connectivity with neighbouring sites and to function as a metapopulation, with high genetic diversity maintained at the regional scale. The current effective population sizes estimated, namely ~34,000, 6000 and 5000 individuals for the DPW, CE and JU populations, are compatible with the population size estimates by MRR in other studies (e.g. in Poland: ~200 individuals per ha, Sielezniew & Nowicki, 2017; in Sweden: ~50 individual per ha, Cassel et al., 2001) given the quantity of suitable habitats available in the Jura predicted by SDM.
The climatic niche of C. hero in France is best captured by the global aridity index (GAI) than by any other bioclimatic variable ( Figure 4), a local index combining precipitation, potential evapotranspiration and temperature (Trabucco & Zomer, 2019). throughout the species distribution range as these glacial relicts are expected to retain more genetic diversity and therefore more evolutionary potential than younger northern populations (Taberlet et al., 1998). In face of global warming and rapidly increasing temperature throughout Europe (Squintu et al., 2021), local adaptations selected in southernmost populations might prove to be useful to evolutionary rescue populations not yet exposed to extreme climatic conditions (Razgour et al., 2019). Coenonympha hero appears to maintain over time through small interconnected populations able to buffer climate change by shifting to cooler habitats (more forested, altitude). Nonetheless, the less diverse and most declining population is the southernmost population (JU) and particular attention is needed to maintain gene flow with the northern core population (DPW). Both climate and land use variables are needed to accurately predict the distribution of C. hero and ongoing changes (warming and increased human footprint) combine as factors of extinction risk, while we found no evidence for C. arcania introgression by mapping C. hero ddRAD data set to C. arcania genome (data not shown), indicating that hybridization with other Coenonympha species is not a driver of C. hero decline. Conservation actions should aim at maintaining and increasing patches of semiopen habitats and forest glades to help the species to maintain in its relictual historical refuge in Jura, but also to favour the connectivity between these marginal southernmost populations and all the other populations in Europe, as it will favour the spreading of warmadapted genes elsewhere in Europe.

ACK N OWLED G EM ENTS
We thank François Dehondt and Conservatoire Botanique National

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw sequences (fastq) analysed in the present study are available at the European Nucleotide Archive repository (http:// www.ebi.ac.uk/ena) and accessible under study accession number PRJEB48550 (individual accession numbers in Supporting Information). The two SNP data sets (in VCF format) used for population genetics and demographic inferences, and the occurrence data sets used for each SDM analysis are provided in Supporting Information. In accordance with the French regulation regarding endangered species (article L.124-4 du code de l'environnement), the precision of coordinates of the localities sampled in the Jura massif was degraded to 10 × 10 km precision but are available from authors upon request.