Past, current, and potential future distributions of unique genetic diversity in a cold‐adapted mountain butterfly

Abstract Aim Climatic changes throughout the Pleistocene have strongly modified species distributions. We examine how these range shifts have affected the genetic diversity of a montane butterfly species and whether the genetic diversity in the extant populations is threatened by future climate change. Location Europe. Taxon Erebia epiphron Lepidoptera: Nymphalidae. Methods We analyzed mtDNA to map current genetic diversity and differentiation of E. epiphron across Europe to identify population refugia and postglacial range shifts. We used species distribution modeling (SDM) to hindcast distributions over the last 21,000 years to identify source locations of extant populations and to project distributions into the future (2070) to predict potential losses in genetic diversity. Results We found substantial genetic diversity unique to specific regions within Europe (total number of haplotypes = 31, number of unique haplotypes = 27, H d = 0.9). Genetic data and SDM hindcasting suggest long‐term separation and survival of discrete populations. Particularly, high rates of unique diversity in postglacially colonized sites in England (H d = 0.64) suggest this population was colonized from a now extinct cryptic refugium. Under future climate change, SDMs predict loss of climate suitability for E. epiphron, particularly at lower elevations (<1,000 meters above sea level) equating to 1 to 12 unique haplotypes being at risk under climate scenarios projecting 1°C and 2–3°C increases respectfully in global temperature by 2070. Main conclusions Our results suggest that historical range expansion and retraction processes by a cold‐adapted mountain species caused diversification between populations, resulting in unique genetic diversity which may be at risk if distributions of cold‐adapted species shrink in future. Assisted colonizations of individuals from at‐risk populations into climatically suitable unoccupied habitat might help conserve unique genetic diversity, and translocations into remaining populations might increase their genetic diversity and hence their ability to adapt to future climate change.

H d = 0.9). Genetic data and SDM hindcasting suggest long-term separation and survival of discrete populations. Particularly, high rates of unique diversity in postglacially colonized sites in England (H d = 0.64) suggest this population was colonized from a now extinct cryptic refugium. Under future climate change, SDMs predict loss of climate suitability for E. epiphron, particularly at lower elevations (<1,000 meters above sea level) equating to 1 to 12 unique haplotypes being at risk under climate scenarios projecting 1°C and 2-3°C increases respectfully in global temperature by 2070.

Main conclusions:
Our results suggest that historical range expansion and retraction processes by a cold-adapted mountain species caused diversification between populations, resulting in unique genetic diversity which may be at risk if distributions of cold-adapted species shrink in future. Assisted colonizations of individuals from at-risk populations into climatically suitable unoccupied habitat might help conserve unique genetic diversity, and translocations into remaining populations might increase their genetic diversity and hence their ability to adapt to future climate change.

| INTRODUC TI ON
Projecting the future geographic distribution of genetic variation within species' ranges and the potential loss of genetic variation from anthropogenic climate change requires understanding of the past, present, and future distributions of species (Wroblewska & Mirski, 2018). Geographic variation in the distribution of genes across a species' range results from a combination of historical and current conditions, which influence patterns of genetic differentiation among populations that are, or have been, geographically isolated, and from colonization bottlenecks during range shifts (Hewitt, 2004). These range shifts and their genetic consequences have primarily been driven by the fundamental niche of a species, or their "climate-envelope," and species' ranges shift to track environmental changes, altering the location of populations and their genetic structure (Hewitt, 2004;McCallum, Guerin, Breed, & Lowe, 2014;Thomas, 2010). The Earth has gone through many climate fluctuations, including glaciations in the Pleistocene and human-induced climate change in the current Anthropocene (Hewitt, 2004;Santer et al., 2019). Future anthropogenic climate warming may further impact species through distribution changes, genetic erosion, and extinctions (Botkin et al., 2007). Cold-adapted/mountain species may be especially vulnerable to future climate changes as they are already restricted to mountain ecosystems where suitable climate space is limited, and loss of genetic diversity within these range-restricted cold-adapted species may reduce their ability to adapt to future changes (Elsen & Tingley, 2015). Understanding how past climatic changes have impacted current genetic structure may allow us to make predictions for the likely extent of genetic loss under future climate change and thereby prioritize at-risk populations for conservation management (McCallum et al., 2014;Wroblewska & Mirski, 2018).
During the last ice age, ice sheets were at their greatest extension 20,000-21,000 years ago, during the Last Glacial Maximum (LGM) (Crowley & North, 1991;Ray & Adams, 2001). During the LGM, species were thought to persist where climatic conditions were buffered, at lower elevations or in more southerly regions (Dapporto et al., 2019;Morelli et al., 2016); however, some studies have shown evidence of species surviving in northern isolated refugia (Provan & Bennett, 2008;Schmitt & Varga, 2012;Stewart & Lister, 2001). Cold-adapted species which currently occur in mountain ecosystems were probably more widespread during the LGM and only became isolated in their current interglacial populations after climate-induced range retraction, although some coldadapted species were already restricted to isolated glacial refugia during the LGM (Schmitt, 2009;Schmitt, Hewitt, & Muller, 2006).
The consequences of past distribution changes will be reflected in current genetic diversity, because contractions and expansions from long-term refugia leave a genetic signature of high diversity in refugia compared to lower diversity in recently colonized populations (Hewitt, 2000;Keppel et al., 2012;Morelli et al., 2016). Thus, understanding historical interactions of cold-adapted species with climate can help us understand current genetic structure and diversity of populations.
Lepidoptera are poikilothermic and therefore sensitive to changes in climate, and those species which are cold-adapted are particularly vulnerable to warmer conditions (Deutsch et al., 2008;Elsen & Tingley, 2015). Some cold-adapted Lepidoptera are experiencing extinctions at their low latitude/elevation margins as the climate deteriorates for these species (Franco et al., 2006;Wilson, Gutierrez, Gutierrez, & Monserrat, 2007).

The Mountain Ringlet
Erebia epiphron is a butterfly found in the mountains of continental Europe and Britain, and its distribution has retracted 130-150 m uphill in Britain over the past five decades due to climate warming (Franco et al., 2006). Therefore, E. epiphron represents a good model organism to understand how past climate-induced changes have impacted current genetic structures of populations, and whether genetic diversity may be lost with further climate-induced local extinctions.
Species distribution modeling (SDMs) are commonly used to project future distributions of species under climate change scenarios (Guo et al., 2017;Urban, 2015) and to develop climate adaptation conservation strategies. These modeling approaches have also been used with paleoclimate data to hindcast past distributions and to understand how they shape current population structures (Smith, Gregory, Anderson, & Thomas, 2013). Phylogeography with genetic techniques can be used to identify divergence between populations and to infer historical distribution patterns and colonization routes (Luquet et al., 2019). Previous studies have shown how a combination of species distribution modeling and phylogeography can provide better understanding of past, present, and future distributions of species and predict the potential loss of genetic diversity resulting from climatic warming Wroblewska & Mirski, 2018;Yannic et al., 2014).
In this study, we use mtDNA sequencing to map the current distribution of genetic diversity of the cold-adapted butterfly, E. epiphron, and also use species distribution modeling to project current, past, and future distributions of the species. We use this genetic and modeling information to determine the distribution of E. epiphron in continental Europe during the Last Glacial Maximum, the locations of glacial refugia, and patterns of subsequent postglacial expansion into northerly latitudes in Britain. We identify populations with unique genetic diversity and examine potential loss of genetic diversity under future climate change scenarios in order to prioritize populations for protection. All relevant fieldwork permissions were obtained. DNA was extracted from 111 individuals with Omega Bio-tek E.Z.N.A.® DNA Isolation Kit following the manufacturer's protocol. For each individual, the head and antennae were removed and placed in 1.5 ml tubes with CLT buffer and proteinase K and homogenized with pellet pestles. A 658-bp fragment of the mitochondrial cytochrome oxidase-I (COI) gene was amplified using the primers LepF (5'-ATTCAACCAATCATAAAGATATTGG-3') and LepR (5'-TAAACTTCTGGATGTCCAAAAAATCA-3') (Hajibabaei, Janzen, Burns, Hallwachs, & Hebert, 2006). PCR amplification of individual DNA samples was carried out in 20 µl reactions which included 1.8 µl of template DNA, 1x PCR reaction buffer (Promega), 1.5 mM MgCl 2 , 0.2 mM of dNTPs, and 1U of Taq DNA polymerase (Promega GoTaq®). PCR conditions used the following profile: 94°C for 2 min (one cycle), 2 min at 94°C, 58°C for 45 s, and 72°C for 1 min (35 cycles), followed by a final extension step of 75°C for 5 min. PCR products were purified and Sanger sequenced with forward and reverse primers using © Eurofins Scientific PlateSeq service and LightRun Tube service. Chromatograms were checked visually using SeqTrace (Stucky, 2012). Additional COI sequences were obtained from a panel of 39 samples collected in England in June 2016 as a part of a whole-genome resequencing project (NERC highlight project NE/N015797/1). Briefly, the complete mitochondrial genome was assembled for each individual sample using the MitoZ toolkit (Meng, Li, Yang, & Liu, 2019) and annotated using the mitos2 web server (Bernt et al., 2013). Low coverage regions (<10) were masked to avoid introducing low-quality SNPs and the COI region was extracted for further analyses.
These 150 sequences along with 65 existing COI sequences from GenBank were combined to create a data set of 215 COI sequences from 13 mountain regions across the species' European range (for sample information see Appendix S1 and map of mountain regions see Appendix S2). These sequences were aligned with ClustalX implemented in MEGA-X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018) and the alignment checked by eye and cropped to the same length (649 bp). Haplotypes were identified and genetic diversity measures were determined using DnaSP 6 (Rozas et al., 2017). Genetic diversity measures included number of haplotypes (H n ), number of unique haplotypes (H u ), haplotype diversity (H d, the probability that two randomly sampled alleles are different), and nucleotide diversity (π, the average number of nucleotide differences per site between sequences (Nei, 1987). A TCS network (Templeton, Crandall, & Sing, 1992) of all haplotypes was constructed using PopArt (Leigh & Bryant, 2015). A CO1 phylogenetic tree was constructed in BEAST (Suchard et al., 2018) of the Erebia genus, outgroups, and the E. epiphron populations. The same methods and CO1 sequences were used from (Peña, Witthauer, Klečková, Fric, & Wahlberg, 2015) using a log-normal relaxed molecular clock, with a birth-death incomplete speciation model for the randomly generated tree prior, and then an uncorrelated log-normal relaxed molecular clock and all the programs other default settings to model the rate of evolution.
The age between Erebia and its sister taxa was set at 37.4 ± 2 Myr, (Peña et al., 2015) to estimate age in divergence between E. epiphron subpopulations.  (Smith et al., 2013). Spatial autocorrelation was tested using Moran's I in R. The butterfly distribution data were at 50-km grid resolution, but the species is likely to be restricted by local climate conditions in each grid square (Smith et al., 2013).

| Using species distribution modeling (SDMs) to map current distribution of E. epiphron
Thus, we included in models only the coldest/warmest and wettest/driest cells (4.5-km resolution) within each 50-km grid, resulting in a total of eight climatic variables being incorporated into our SDMs (see Appendix S3). 50 × 50-km grid cell resolution data are appropriate for our model building to address biogeographic questions at regional scales, because we are interested in changes in the distribution of the study species over long periods of time (i.e., millennia), rather than shorter-term changes at individual sites. This 50-km spatial resolution also ensures that the pseudoabsences (i.e., locations where E. epiphron is assumed to be absent) are more accurate representations of true absences, because these We used the mean receiver operating characteristic (ROC) value to evaluate each model, with a threshold of ROC > 0.85 for inclusion of models within the ensemble model. We restricted pseudoabsences to locations within a buffer of 250 km around presence data points to avoid placing absences in mountain systems with potentially suitable climate space that are not currently occupied by the species (e.g., Scandinavia) (Akcakaya & Atwood, 1997;Hirzel, Helfer, & Metral, 2001). Models were generated using 70% training data and 30% testing data (Franklin, 2010;Huberty, 1994).

| Hindcasting past distributions and identifying glacial refugia
We incorporated paleoclimate data into our ensemble SDM for the eight climate variables representing the coldest/hottest and driest/wettest locations within each 50-km grid square. Data for climate projections over the last 21,000 years were downloaded from PaleoView (2.5 × 2.5° (latitude/longitude) grid) (Fordham et al., 2017), and downscaled to match the resolution of the current climate data (2.5 arc minutes), using established methods (Mitasova & Mitas, 1993
Unique haplotypes were assumed to be at risk if all 50-km grid squares in one of the 13 mountain regions were predicted to become climatically unsuitable in the future (based on binary presence or absence threshold probability values from the ensemble SDM output). We set the threshold value as the probability value associated with the low elevation climatic range edge E. epiphron in its current range (low elevation range boundary in England; threshold probability = 0.49). Using this threshold, model probabilities were converted into the presence/absence to show grid squares with no change over time (i.e., population persistence), grids predicted to become climatically unsuitable (i.e., extinction), and grids predicted to become climatically suitable (i.e., colonization). Haplotype risk (H r ) was calculated as the number of unique haplotypes at risk in each of the 13 mountain regions (Figure 1a) due to projected loss of all climatically suitable areas within a region in the future.

| Current haplotype diversity across 13 mountain regions in Europe
From our 215 mtDNA samples, we identified 31 mtDNA haplotypes across Europe, including 27 haplotypes unique to a specific mountain region ( Figure 1a, Table 1). The high frequency of unique haplotypes across Europe suggests low levels of allele sharing.
There was also high genetic differentiation between populations (AMOVA, ϕ = 0.76, p < .001) and the divergence between some of these populations is dated before the Last Glacial Maximum (phy- Mountains, and has one unique haplotype (haplotype 29) which diverged from haplotype 16 by one substitution (Figure 1b).

| Hindcasting past distributions of E. epiphron and identifying areas of long-term survival
Climate suitability in the LGM (21,000 years before present) showed overlap of climatically suitable areas (at 50-km grid resolution) with many locations currently occupied by E. epiphron, as well as some southerly locations (Figure 2e). This overlap was confirmed when all 21 SDM outputs for each 1,000-year time period up to the present day were combined to show long-term climatic stability since the LGM (Figure 2f). These climate stability maps provided evidence that the locations of glacial refugia were in areas of high topographic variation within the species' current distribution in continental Europe.

| Projecting future distributions and loss of genetic diversity
As expected for a cold-adapted species, SDM outputs from both high and low future climate change scenarios project that many extant E. epiphron areas will have reduced climate suitability in the future (38%-64% loss of 7,000-km 2 occupied grids across Europe) ( Figure 3, Table 1  epiphron does not currently occur in Scandinavia, our models predict that this area will remain stable in climate suitability in the future.

| D ISCUSS I ON
By using species distribution modeling and mtDNA analyses, we explore the past, present, and potential future distributions of genetic diversity in the cold-adapted species E. epiphron. We identify high levels of genetic differentiation across Europe and found evidence of long-term climate suitability in many of these regions since the LGM, which suggests these climatically stable regions were refugial areas of long-term survival by our study species over the last 21,000 years and potentially longer-term areas of persistence over

| Limitations
This study has potential limitations, which are inherent in species distribution modeling, especially when projecting into different climates (Buisson, Thuiller, Casajus, Lek, & Grenouillet, 2010). We did not have suitable data to include sampling effort formally into our models, and so, the areas outside of the current E. epiphron distribution are considered "pseudo-absences" rather than "true" absences.
However, other butterfly species have been recorded in these squares (Lepidopterists have visited these squares) without recording E. epiphron as present, and hence, the proportion of false absences in the data is likely to be very low at the spatial ( Future work could use next-generation sequencing to further test our hypotheses and to model-specific genetic-climatic relationships in the future (see Bay et al., 2018).
Our analyses suggest that entire mountain regions of the butterfly's distribution could be lost under future climatic change, but it is possible that isolated populations could survive in particular microhabitats, at least temporarily. However, these localized populations may not contain all of the genetic variation currently present in the wider region, and overtime, these refugial populations may gradually lose genetic variation and viability (e.g., through inbreeding), and so, they may not persist in the longer term due to their isolation (metapopulation failure). A variety of processes may lead to the loss of genetic diversity following isolation, and there can sometimes be a delay in genetic loss following population decline (Kadlec, Vrba, Kepka, Schmitt, & Konvicka, 2010). For example, the sister species of E. epiphron, Erebia orientalis, is very localized and currently occurs only in the Eastern Balkans and is genetically homogeneous, potentially putting it at risk of inbreeding depression (Hinojosa, Monasterio, Escobes, Dinca, & Vila, 2019). Therefore, our model projections should be seen as representing much longer-term regional-scale expectations, rather than short-term predictions at the local population or microhabitat scale. We believe that our conclusions about the long-term (LGM to present) continental-scale dynamics of E. epiphron are robust and that this knowledge of the past helps frame future risks and provides information for conservation management.

| Long-term survival resulting in unique genetic diversity in cold-adapted species
SDM outputs provide evidence that our exemplar cold-adapted study species occurred in disjunct regions throughout the period from the LGM to the present day, based on the distribution of suitable climate; the genetic data confirm likely separation not only since the LGM, but most probably over much longer periods and successive glacial-interglacial cycles. For mountain species, limited gene flow between the disjunctive parts of their range during glacial and interglacial periods results in divergence and unique haplotypes, unlike lowland European species which colonized northwards from their glacial refugia, and where large parts of the current geographic ranges often share haplotypes (Hewitt, 2004).
LGM separation of populations has also been identified in mountain plants and other invertebrates (Bettin, Cornejo, Edwards, & Holderegger, 2007;Huck, Budel, & Schmitt, 2012;Margraf, Verdon, Rahier, & Naisbit, 2007;Pauls, Lumbsch, & Haase, 2006). The numbers of glacial-interglacial cycles over which populations have remained disjunct remain unclear, but some studies have indicated divergence dates covering several glacial-interglacial cycles or even predating the Pleistocene (Hewitt, 2000). The reality is likely to be more complex with areas of persistent separation, but with occasional links between them (i.e., rare gene flow or brief periods of connection), as indicated by the distributions and relatedness of haplotypes in Figure 1.

| Unique haplotypes in populations derived from northern cryptic refugia
Following the LGM, the ice retreated in northern Europe and many species colonized northwards, for example, via the land bridge between continental Europe and Britain, which was present until sea level rise ~ 7,000 years before present (Sturt, Garrow, & Bradley, 2013). The locations of southerly glacial refugia, which are thought to be the main sources of colonizations, have been debated extensively, with proposed glacial refugia in the Iberian Peninsula, Italy, and the Balkans (Hewitt, 2000), and this has recently been reinforced in European butterflies (Dapporto et al., 2019). However, there is also evidence for more northern cryptic refugia based on fossil, pollen, and genetic evidence (Birks & Willis, 2008;Provan & Bennett, 2008;Stewart & Lister, 2001), where species apparently persisted at higher latitudes in sheltered locations with suitable microclimates (Stewart, Lister, Barnes, & Dalen, 2010). However, most cryptic refugia described to date have been for relatively warm-adapted species.
Here, we present evidence for the existence of northern cryptic population(s) for cold-adapted species during the LGM, based on high unique genetic diversity of the present-day E. epiphron populations in England, an area that was beneath an ice sheet at the LGM (Hughes et al., 2016). since the LGM, although this seems highly unlikely given the short time period for one to three mutations to occur (Figure 1b). It is possible that these LGM populations were situated on land that is currently below sea level, at an edge of the glacier, or in sheltered low elevation microclimates on land. Multiple colonization events have also been shown in other taxa in the UK (Piertney et al., 2005), and the locations of cryptic refugia during the LGM are assumed to be ice-free areas in southern England (Bocherens, Fogel, Tuross, & Zeder, 1995, Lister, 1984, northern Scotland (Bennett, 1984), and southern Ireland (Montgomery, Provan, McCabe, & Yalden, 2014 (Weeks, Stoklosa, & Hoffmann, 2016) and thus may be particularly vulnerable to future climatic changes. Our SDMs project loss of suitable climate for E. epiphron in many locations in Europe, especially in regions with predominantly low elevation populations and few opportunities to shift uphill to high elevation, which could result in loss of genetic diversity. However, our projections of range retraction do not take into account any potential of populations to adapt to warmer temperatures in situ (Franks & Hoffmann, 2012). Future loss of genetic diversity has also been predicted in other species (Alsos et al., 2012;Beatty & Provan, 2011;Yannic et al., 2014), and rates of loss of genetic diversity in wild populations since the industrial revolution (Leigh, Hendry, Vázquez-Domínguez, & Friesen, 2019) are consistent with our projections.

| Conservation interventions to mitigate climate-driven genetic erosion
Conservation management and adaptation could protect cold-adapted populations and safeguard unique genetic diversity from climate change (Mawdsley, O'Malley, & Ojima, 2009 (Thomas, 2011).
For E. epiphron, our SDMs reveal areas in Scandinavia to be climatically suitable, although the species does not occur there, and climate is predicted to increase in suitability in future in Scandinavia for E. epiphron ( Figure 3) and for other Erebia species (Settele et al., 2008). However,  (Aitken & Whitlock, 2013). This could involve moving warm-adapted individuals into cooler populations to increase their adaptive capacity as the climate warms (Weeks et al., 2011).
However, moving locally adapted populations may result in outbreeding depression and maladaptation, negatively impacting populations (Weeks et al., 2011), although some genetic rescue interventions have resulted in increases in populations, and alleles associated with local adaptation were not lost following gene flow (Fitzpatrick et al., 2020).
Genetic conservation interventions for insects, and specifically butterflies, have been rarely implemented, although increasing habitat connectivity has led to genetic rescue of populations (Jangjoo, Matter, Roland, & Keyghobadi, 2016) and genetic data have been used to inform on reintroductions (Dinca et al., 2018). There is no evidence of attempted genetic rescue via translocations of butterflies, although translocating individuals is a genetic conservation strategy which may be important in ensuring future survival and adaptability of populations under climate change. As with translocations, these conservation options may also be controversial, but could remove the need for ongoing intervention and management at sites with declining populations (Weeks et al., 2011). We recommend that before the implementation of any climate adaptation strategy, populations are closely monitored to determine if populations are retracting and likely to become extinct in areas that are becoming too warm for the species. In addition, individual species' assessments are required to assess the genetic diversity of populations and any local adaptation, which would determine the most appropriate conservation strategy.

| CON CLUS IONS
The genetic diversification of cold-adapted mountain species, as demonstrated in our study species E. epiphron, has been shaped by Pleistocene glaciations, the locations of long-term survival of populations, and colonization patterns after the LGM, resulting in unique genetic diversity in isolated populations. Mountain and cold-adapted species are vulnerable to future climate warming, and we predict E. epiphron will lose 38%-64% of its range in the future, especially at low elevations. The uniqueness of genetic diversity contained in these populations could be at risk depending on the severity of future climate change. Conservation strategies such as translocation could ensure the survival of these cold-adapted species, but more research is needed on the likely effectiveness of such approaches.

ACK N OWLED G M ENTS
We We also thank Phil Platts for advice on SDM modeling and Cróna Mc Monagle and Ian Wilson for assistance with UK data collection.
We thank Ilik Saccheri for advice on data analysis and access to additional mtDNA data obtained via the NERC highlight project and Georgina Palmer for data collection. We thank all landowners of field sites, Natural England, and Scottish Natural Heritage for providing permission to sample. We thank the two anonymous reviews for their critical comments and suggestions for improvement in the initial draft.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
GenBank and BOLD accession numbers for each mtDNA COI sequence can be found in Appendix S1.