Large‐ and small‐scale geographic structures affecting genetic patterns across populations of an Alpine butterfly

Abstract Understanding factors influencing patterns of genetic diversity and the population genetic structure of species is of particular importance in the current era of global climate change and habitat loss. These factors include the evolutionary history of a species as well as heterogeneity in the environment it occupies, which in turn can change across time. Most studies investigating spatio‐temporal genetic patterns have focused on patterns across wide geographic areas rather than local variation, but the latter can nevertheless be important particularly in topographically complex areas. Here, we consider these issues in the Sooty Copper butterfly (Lycaena tityrus) from the European Alps, using genome‐wide SNPs identified through RADseq. We found strong genetic differentiation within the Alps with four genetic clusters, indicating western, central, and eastern refuges, and a strong reduction of genetic diversity from west to east. This reduction in diversity may suggest that the southwestern refuge was the largest one in comparison to other refuges. Also, the high genetic diversity in the west may result from (a) admixture of different western refuges, (b) more recent demographic changes, or (c) introgression of lowland L. tityrus populations. At small spatial scales, populations were structured by several landscape features and especially by high mountain ridges and large river valleys. We detected 36 outlier loci likely under altitudinal selection, including several loci related to membranes and cellular processes. We suggest that efforts to preserve alpine L. tityrus should focus on the genetically diverse populations in the western Alps, and that the dolomite populations should be treated as genetically distinct management units, since they appear to be currently more threatened than others. This study demonstrates the usefulness of SNP‐based approaches for understanding patterns of genetic diversity, gene flow, and selection in a region that is expected to be particularly vulnerable to climate change.


| INTRODUC TI ON
An important aim in contemporary evolutionary biology and ecology is to understand the dynamics of genetic diversity and structure of a range of representative species in space and time, based on demographic history and geographic factors in the environment (Després et al., 2019;Jackson et al., 2018).
Understanding such spatial and temporal variation is of particular relevance in the current era of dramatic climate change and habitat loss (Butchart et al., 2010;Lucena-Perez et al., 2020;Wilson et al., 2016). Spatio-temporal features form biogeographic patterns such as through barriers or corridors for dispersal, affecting the ability of organisms to locate suitable habitat and expand their ranges (Hewitt, 2003;McRae & Beier, 2007;Sheth et al., 2020;Zhivotovsky, 2016). However, the impacts of such features are often species-specific, depending on the ability of species to disperse and adapt, and on the scale and extent of (un-) suitable habitat and barriers (Hewitt, 1999;Keyghobadi, 2007;Sheth et al., 2020).
On a temporal scale, past climates influence the present-day distribution of species (Parmesan et al., 1999;Schmitt et al., 2006;Thomas & Lennon, 1999). In Europe, populations of many species have been strongly affected by oscillations of cold and warm periods in the Pleistocene, which have resulted in latitudinal and altitudinal range shifts (Coope, 1970;Hewitt, 1996Hewitt, , 1999Schmitt et al., 2006). During cold periods, many European species became restricted to southern refuges from which they expanded northwards in warm periods (Després et al., 2019;Strandberg et al., 2011;Tab erlet et al., 1998). On a spatial scale, environmental heterogeneity is the main factor influencing the population structure of species (Kokko & López-Sepulcre, 2006;Lowe & McPeek, 2014;Yang et al., 2017). Heterogeneous environments, characterized by topographic barriers, landscape clines, or unsuitable habitats, affect population connectivity and the distribution of species (Jackson et al., 2018;Lowe & McPeek, 2014;Manel et al., 2003).
Geographic structures at the spatial scale can be further divided into large-and small-scale structures. Small-scale structures that might influence local or regional patterns of diversity comprise natural barriers like rivers, lakes, mountain ridges, gorges, and forests, as well as anthropogenic barriers like roads, railways, agricultural plots, and settlements (e.g., Heidinger et al., 2013;Miles et al., 2019;Nilsson et al., 2013;Sheth et al., 2020;Storfer et al., 2010). Large-scale structures include mountain ranges, oceans, and deserts. In Europe, the Alps and the Pyrenees have formed, since the Last Glacial Maximum, large-scale geographic barriers that are difficult to cross for less mobile species (Hewitt, 1996(Hewitt, , 1999Tab erlet et al., 1998). These mountain ranges currently harbor hybrid zones and endemic taxa across heterogeneous landscapes composed of unsuitable and suitable habitats in close proximity, created by small-scale geographic features (Dagnino et al., 2020;Dirnböck et al., 2011;Hewitt, 2004;Martin et al., 2002). They provide a useful region for investigating the impact of geographic features on population structure and connectivity. Furthermore, alpine environments involve steep gradients in environmental conditions such as temperature, oxygen concentration, and ultraviolet radiation, affecting species assemblages (Cheviron & Brumfield, 2012;Montero-Mendieta et al., 2019) including insects and their host plants (Hodkinson, 2005;Horn et al., 2006;Montero-Mendieta et al., 2019). Thus, specific adaptations to high-altitude environments can be expected (Cheviron & Brumfield, 2012;Karl et al., 2008Karl et al., , 2009Polato et al., 2017).
Our focus here is on the way these various structures have influenced the population genetic structure and patterns of genetic variation. Most research on spatio-temporal patterns has focused on ecosystems, ecological communities, and species distributions and much more rarely on genetic diversity within species (Coates et al., 2018;Hoban et al., 2020;Miraldo et al., 2016). A consideration of genetic diversity and population structure is important for several reasons. First, reduced genetic diversity may interfere with a species' ability to respond to changing environments (Willi et al., 2006); it may also be a signature of a persistently small population size that can result in inbreeding (and hence inbreeding depression) in a population (Bijlsma & Loeschcke, 2012;Day et al., 2003;Frankham, 1995;Keller & Waller, 2002). Second, genetically distinct populations may be useful for delineating management units for conservation purposes (Coates et al., 2018;Moritz, 1994). If populations of a given species are sufficiently differentiated from other populations, they may form evolutionary significant units (ESUs; Frankham, 2010;Moritz, 1994;Ryder, 1986) that have different evolutionary trajectories, although this requires some caution because genetic drift can also result in genetically unique populations with limited evolutionary capacity (Weeks et al., 2016).
Here, we examined population genetic structure and diver- heterogeneity (Martin et al., 2002). Lycaena tityrus forms closed populations and is likely to have a limited dispersal ability which is strongly affected by dispersal barriers at small scales (Trense et al., 2021). Until recently, genetic variation in butterfly populations has been investigated through markers such as allozymes, microsatellites, or a small number of sequenced mitochondrial and nuclear genes, but inferences based on these markers have been limited by the fact that they sample only a small part of the genome (e.g., Maresova et al., 2019;Schmitt et al., 2006;Ugelvig et al., 2012). New molecular techniques now provide thousands of SNP markers across the genome allowing the identification of past refuges, genetic diversity, geographic barriers, and gene flow in butterflies with unprecedented rigor (Fountain et al., 2018;Nève, 2009).
In this study, we used genome-wide single nucleotide polymorphisms (SNPs) to investigate genetic diversity and geographic population structure of L. tityrus across the eastern European Alps. In a previous study using the same species, we found high gene flow within a single Alpine valley, but also genetic structuring linked to ravines, forests, roads, and altitude (Trense et al., 2021). Here, we analyzed a much broader spatial scale by sampling L. tityrus individuals in 30 different valleys in the European Alps, while the previous study was confined to only one valley. In the present study, we test whether (a) populations are structured by the presence of high mountain ridges and/or large river valleys potentially comprising barriers to dispersal, (b) current patterns of genetic structure are influenced by long-term processes such as postglacial range expansions, (c) different evolutionary lineages are present within this species, and also whether (d) populations from different altitudes show signatures of local adaption based on an analysis of outlier loci.

| Study organism and population sampling
Here, we used the alpine subspecies L. tityrus subalpinus (Lepidoptera: Lycaenidae). Lycaena tityrus is a widespread butterfly of the temperate zone with a range from Western Europe to central Asia (Ebert & Rennwald, 1991). It inhabits different kinds of grassland, including moist and dry meadows, sandy heathland, bogs, and open woodland (Settele et al., 2008). The principal larval host plant is Rumex acetosa L., but some congeneric plant species (e.g., Rumex acetosella L., Rumex scutatus L.; Ebert & Rennwald, 1991, Settele et al., 2008, Tolman & Lewington, 2008 are also used. Lycaena t. subalpinus is confined to higher altitudes of the European Alps and some other mountain ranges, where it is relatively widespread and has one generation a year (Tolman & Lewington, 2008). The altitudinal distribution of L. t. subalpinus ranges from 1,200 to 2,500 m a.s.l. (Tolman & Lewington, 2008). Thirty populations were sampled in Austria (Salzburg, Tyrol, Vorarlberg), Italy (South Tirol, Lombardy), and Switzerland (Graubünden), spanning an altitude from ca. 1,260 to 2,110 m a.s.l. (Figure 2, Table 1). We caught nine males per population in the summers of 2018 and 2019. All 270 individuals used in the current analysis were stored in liquid nitrogen until DNA extraction.
Here, we re-used the extracted DNA of 23 from the 186 individuals F I G U R E 2 Geographic distribution of the 30 sampling locations for Lycaena tityrus across the European Alps. Details are given in Table 1. The black lines indicate the national borders. The map was generated with QGIS version 3.14 (www.qgis.org) tested in Trense et al. (2021). However, sequencing of these samples was done again together with the new samples.

| ddRAD library preparation
From each male, we used head and thorax for the extraction of genomic DNA with the E.Z.N.A. ® Tissue DNA Kit (Omega Bio-tek).
We followed the manufacturers' instructions but included an additional step of RNase A treatment. Afterward, we applied doubledigest restriction site-associated DNA sequencing (ddRADseq) following Trense et al. (2021) and used the restriction enzymes NlaIII and MluCI. We pooled individuals into four libraries, each containing DNA fragments from individuals with different adapter pairs. Libraries were cleaned with 1.5× volume of Sera-Mag beads. We selected 250-400 base pair (bp) fragments using a 2% gel cassette and Pippin-Prep software 4.3 (Sage Science). This was followed by a polymerase chain reaction (PCR) enrichment in a 10µl reaction with 1 µl of size selected DNA, 2 µl 5× Phusion ® HF Reaction Buffer, 0.2 µl dNTPs (10 mM), 0.1 µl (100 units) Phusion® HF DNA Polymerase (New England BioLabs Inc.), nuclease-free water, and 2 µl (10 µM) each of Illumina P1 and P2 primers (Peterson et al., 2012). The profile of thermal cycling consisted of denaturation at 98°C for 30 s, followed by 12 cycles with 10 s at 98°C, 30 s at 65°C, and 70 s at 72°C, and an extension for 5 min at 72°C. For the final library, we used seven PCRs. The four libraries were sequenced on an Illumina NovaSeq 6000 platform generating 150-bp paired-end reads.

| SNP calling
For calling SNPs, we used a de novo pipeline without a reference genome, namely, the process_radtags program within StackS version 2.3 (Catchen et al., 2011(Catchen et al., , 2013. Here, we demultiplexed the sequence SNPs were subsequently called. We filtered the data to retain SNPs, which were present in more than 50% of the 270 individuals, and then retained the first SNP in each ddRAD locus. Note that some ddRAD loci can be located in the same gene. Further filtering with the program vcftoolS version 0.1.11 (Danecek et al., 2011) was done to retain loci that were at Hardy-Weinberg equilibrium, had a minor allele count of >2 (Linck & Battey, 2019), missing data of <5%, and a mean depth of 20-45. For filtering, individuals of the same sampling location were treated as belonging to one population but excluding two individuals of population 9 because of incomplete sequencing.
The final dataset contained 13,455 SNPs from a total of 99,717 polymorphic RAD loci.

| Analysis of relatedness
To avoid biased population analyses, we used the Triadic Likelihood Estimator (TrioML, giving a relatedness coefficient) in coanceStry version 1.0.1.9 (Wang, 2007(Wang, , 2011 to calculate the pairwise relatedness among individuals. Individuals with a TrioML estimator of ≥0.250 are most likely full-siblings and between 0.125 and 0.249 probably half-siblings. Real ancestry cannot be assessed by TrioML though, but only genetic relatedness. We found several related individuals including seven probable pairs of full-siblings. One full-sibling per probable pair was excluded from downstream analyses. To confirm putative full-and half-siblings, we additionally used the Specific Hypothesis Test function with 9,999 permutations implemented in the program ml-relate (Kalinowski et al., 2006).

| Isolation by distance and resistance
To test for isolation by distance, we calculated pairwise genetic distances between individuals according to Bray Curtis (Bray & Curtis, 1957), using the package "vegan" (Oksanen et al., 2019) in r version 3.5.2 (R Core Team, 2015). Pairwise, individual-based geographic distances were calculated as (a) Euclidean distances with the function PointDistance (plane) in the package "raster" (Hijmans & van Etten, 2016) in R, based on the UTM system, and as (b) altitudinal distances by a distance matrix. We tested for correlations between genetic and geographic distance (isolation by distance) as well as altitudinal distance by using Mantel tests (using the Mantel function from the "vegan" package in r; Mantel, 1967)  to receive resistance surfaces and distances for these features and followed Cushman et al. (2006) to analyze the relationship between genetic and landscape distance matrices by using the simple and partial Mantel test in the package "ecodist" (Goslee & Urban, 2007) in R with 9,999 permutations. We performed five partial Mantel tests for each landscape feature to analyze the relationships between genetic and landscape distance matrices, partialling out the effect of other landscape distance matrices. To account for the nonindependency of records from one individual, we additionally ran maximumlikelihood population effects (MLPE) following Clarke et al. (2002) and Peterman (2018).

| Analysis of molecular variance (AMOVA)
We used an hierarchical AMOVA in arlequin version 3.5 (Excoffier & Lischer, 2010) to assess the genetic variance among the 30 populations, among individuals within these populations and within individuals.

| Analysis of genetic structure
We calculated observed and expected heterozygosity, inbreeding coefficient, and global F ST for all 30 populations using the packages "adegenet" (Jombart, 2008;Jombart & Ahmed, 2011) and "hierfstat" (Goudet, 2005) in r. The pairwise F ST values were calculated in arlequin. Furthermore, we identified genetic clusters using the function SNMF (sparse non-negative matrix factorization) in the package "LEA" (Frichot & François, 2015) in R. We tested 31 ancestral populations (K = 1-31) with 100 repetitions for each K. The minimal cross-entropy value was selected to visualize the genetic clustering. Based on the results of the eemS analysis which indicated that major river valleys and high mountain ridges may comprise genetic barriers for this species (cf.  (Bray & Curtis, 1957). To test how well the different groups were statistically supported, we ran distance-based redundancy analyses (dbRDA; Legendre & Anderson, 1999) with the r package "vegan" (Oksanen et al., 2007).
For analyzing the population demographic history, we used the approximate Bayesian computation (ABC) method as implanted in Diyabc version 2.1.0 (Cornuet et al., 2014). Based on the results of the SNMF analysis, four population clusters (western, central, eastern, and southeastern populations) were defined (Figure 5b). For the four population clusters, we ran seven scenarios ( Figure S1.1, Appendix S1). We set the conditions of the prior time distributions to t1 < t2 to avoid incongruences in the simulated genealogies. The overall performance of scenarios was assessed with a principal component analysis using 100,000 simulated datasets and the observed data. We then assessed the posterior probability of the scenarios with a logistic regression procedure based on the 1% closest simulated datasets compared to the observed data. We used the posterior distributions for the best supported scenario for simulating 1,000 pseudo-observed datasets to assess whether this model could successfully reproduce the observed data.

| Outlier loci analysis
To identify potential outlier loci linked to altitude, we used the altitude of each sampling location; that is, the individuals of one location were assigned the same altitude. Three different F ST outlier analyses were performed to minimize the risk of false positives, namely, baypaSS version 2.1 (Gautier, 2015), FDIST2 in arlequin (Excoffier & Lischer, 2010), and bayeScan version 2.1 (Foll & Gaggiotti, 2008 Note: Factors used are the sampled locations (populations) as well as different clusters as defined by major river valleys, major river valleys in combination with high mountain ridges, and high mountain ridges. Effect sizes are given as partial Eta squared (η 2 ), and significant p-values are given in bold.

| Isolation by distance and resistance
Genetic distances among all individuals (Bray-Curtis) ranged from 0.086 to 0.137, and the geographic (Euclidean) distance between sampled populations ranged from 6.2 to 254.4 km. Mantel tests for isolation by distance showed a significant correlation between genetic and Euclidean distances (r = 0.574, p < .0001; Figure 3), but not between genetic and altitudinal distances (r = 0.033, p = .170).
The eemS analysis showed several dispersal barriers (Figure 4a) The results of the IBR analyses showed that the highest correlations for the simple Mantel tests were for altitude, slope, and tree cover ( Table 3). The partial Mantel tests showed that roads and water bodies affected the genetic structure of L. tityrus in the European Alps (Table 3). The MLPE models for isolation by resistance, finally, showed that all environmental factors contributed to population genetic structure (Table 4).

| Population structure
Regarding molecular indices, populations 3 and 30 showed the highest and lowest observed heterozygosity, respectively (Table 5). The highest expected heterozygosity and inbreeding coefficient were found in population 2, while the lowest values were observed in 29.
The program arlequin provided F ST values between each population pair, and from these values, we calculated the mean F ST value for each population. Based on arlequin and bayeScan, the highest mean F ST value was found in population 29, and the lowest one in population 13 (also 10 and 12 in bayeScan; Table 5, for details see Table S1.1, Appendix S1). An AMOVA revealed that the highest amount of genetic variation occurred within individuals (94.6%), followed by a significant structuring among populations (4.5%), while variation among individuals within populations was not significant (Table 6).
Distance-based redundancy analyses (dbRDA) revealed that all a priori groupings were statistically supported (Table 2). Based on effect sizes, the grouping according to sampling location is most strongly supported (Figure 6a), followed by the grouping according to river valleys in combination with high mountain ridges (Figure 6b).
The DIYABC analysis showed the strongest support for scenario 5 with a posterior probability of 0.762 and confidence intervals of 0.667-0.860 ( Figure S1.1, Appendix S1). Scenario 5 considered that the western, central, and eastern population clusters evolved separately, and that the southeastern population cluster split from the eastern population cluster. All other scenarios received only little support with posterior probabilities below 0.132.

| Loci under selection according to altitude
In total, we identified 514 SNPs as potential outliers according to altitude with at least two analyses. baypaSS, FDIST2, and bayeScan detected 673, 1,866, and 269 significant outlier SNPs, respectively.
Thirty-six out of the 514 putative outlier SNPs had a match in omicSbox (Table 7). Of these, ten could be assigned to membrane transport, Note: Here, we presented only the results of 6-and 5-factor models. Note that model including all six factors had the lowest AIC. Models included genetic distance (GD) according to Bray-Curtis (Bray & Curtis, 1957) as response variable, altitude, grass cover, road, slope, tree cover, and/or water bodies as fixed factors, and Lycaena tityrus individuals (n = 261) as a random factor. We used the Akaike information criterion (AIC) as an indicator of model quality. The marginal R 2 (R 2 m) gives the proportion of variance explained by fixed factors, and the conditional R 2 (R 2 c) by fixed and random factors.
four to membrane/peptidases, three to DNA binding, two each to metal binding, transcription and one each to post-translational modification, nucleotide binding, phosphatase inhibitor, kinase activity, translation, tyrosinase, lipid binding, microtubule binding, ATP binding, protein phosphatase regulator activity, nucleosome assembly, protein dimerization activity, and protease. In LFMM, we detected 16 and 13 outlier SNPs being related to altitude and temperature variables, respectively (data not shown). Nine outlier SNPs were found to be linked to both altitude and different temperature variables.  Figure 5).  (Trense et al., 2021). Note here that the results of simple and partial Mantel tests need to be interpreted with caution because they suffer from high Type I error rates (Cushman et al., 2013;Guillot & Rousset, 2013). However, the MLPE models indicated that all factors contributed to the population genetic structure of L. tityrus in the European Alps. In particular, high mountains and large river valleys seem to represent important barriers hampering gene flow, which is supported by our dbRDA and eemS results (  Figures 4 and 6). High mountains above ca. 2,300-2,500 m a. s. l. do not comprise suitable habitats for L. t. subalpinus, especially when covered by glaciers. Accordingly, mountain ranges may act as effective barriers for dispersal and thus gene flow among bee, butterfly, and grasshopper populations at large scales (Britten et al., 1995;Després et al., 2019;Hewitt, 1996;Jaffé et al., 2019).

| D ISCUSS I ON
Our study shows that this applies also to smaller spatial scales within mountain ranges. Likewise, large rivers such as Danube, Isère, Rhine, and Rhone are known to affect large-scale genetic patterns, indicated by genetic differentiation of insect populations from either side of rivers (Cupedo & Doorenweerd, 2020;Mardulyn, 2001;Schmitt et al., 2007). However, smaller rivers were not found to comprise a barrier for L. t. subalpinus in one Austrian valley (Trense et al., 2021).
Hence, whether a river valley displays a barrier to gene flow depends Note: Given are the percentages of variation among populations, among individuals within populations, and within individuals. Significant p-values are given in bold.

TA B L E 6
Results of an AMOVA for 30 populations of Lycaena tityrus F I G U R E 5 Genetic clusters as computed with a sparse non-negative matrix factorization (SNMF), based on 13 455 SNPs and 30 Lycaena tityrus populations. Four genetic clusters were identified. In (a), each vertical line represents one individual, and its likely assignment to a specific genetic cluster encoded by different colors. In (b), each circle represents one population (cf. Figure 2), with the percentage of individuals assigned to a specific genetic cluster on its size, elevation, and climatic conditions (Ćosić et al., 2013;Link et al., 2015;Muñoz-Mendoza et al., 2017). In L. t. subalpinus, only large river valleys at low altitudes seem to act as barriers. These currently do not offer suitable habitat for the species, mainly due to agricultural intensification, although they may have historically.
Cluster analyses revealed a western, eastern, southeastern, and a central cluster, which may indicate different refuges from which the Alps were recolonized independently after the last glacial period. A similar pattern was found in Erebia butterflies (Schmitt et al., 2016), and an eastern and western genetic cluster occurred also in alpine Erebia alberganus  and in Drusus discolor (Pauls et al., 2006). Accordingly, the Diyabc analysis suggested the occurrence of three glacial refuges in the western, central (perhaps along the river Etsch), and eastern Alps, with the southeastern population cluster originating from the eastern refuge. Interestingly, the different genetic clusters in our study are characterized by striking differences in genetic diversity, being high in the western, intermediate in the central, and low in the eastern and especially in the southeastern cluster. These data suggest that the largest glacial refuge of L. t. subalpinus was located in southwestern Europe or at the southwestern edge of the Alps. This refuge may have included the lowlands between the western Alps and the Pyrenees, as the alpine subspecies occurs in both mountain areas, which may thus indicate relict populations (Tolman & Lewington, 2008). In addition, admixture of different western (and southwestern) refuges or hybridization with low-altitude populations, that is, L. tityrus tityrus, may have contributed to the high genetic diversity in the western Alps. Generally, the southwestern Alps are known to be important evolutionary centers for arctic-alpine species Schmitt, 2009;Schönswetter et al., 2005). The lower genetic diversity in the central and especially eastern genetic clusters suggests a smaller refuge size (or rather smaller effective population sizes) as compared with the supposed western refuge, in combination with founder effects for the southeastern cluster (Austerlitz et al., 1997;Fayard et al., 2009 decrease in genetic diversity from the core to the edges of a species' range (Brown et al., 1995;Eckert et al., 2008). We can also not rule out that our eastern sampling points were further away from the eastern refuge, for instance located in the Slovenian Alps, explaining the decreased genetic diversity.
In addition to effects of the last glaciation, the patterns found could be also influenced by more recent developments. In the dolomites, for instance, L. t. subalpinus is much rarer than further north, which is most likely due to differences in geology. The dolomites with their karst landscapes have much dryer soils strongly reducing the abundance of Rumex plants compared with the central Alps with their crystalline bedrocks. Concomitantly, population size and connectivity are likely to be low in the dolomites, which may also explain the low genetic diversity of these populations. Whether this may also apply to (parts of) the eastern cluster is currently unclear.
Thus, while the large-scale population structure of L. t. subalpinus is likely shaped by the last glacial period, local and regional factors may also be important.
With regard to the conservation status of the taxon, L. t. subalpinus represents an evolutionary significant unit for conservation, based on substantial morphological, ecological, and genetic differences compared to the lowland form L. t. tityrus (Karl et al., 2008(Karl et al., , 2009Tolman & Lewington, 2008). Whether this also applies to the inner-Alpine genetic clusters found here is open to debate. We suggest that the dolomite populations should be considered as their own entities, as these show strong genetic divergence in the SNMF analysis and are likely more threatened than the central Alpine populations. For the protection of the taxon L. t. subalpinus though, conservation efforts should be concentrated in the western Alps with their genetically diverse populations and the importance of genetic diversity in future adaptation.
As predicted, we found footprints of selection to the unique alpine environments including differences in temperature, oxygen, ultraviolet radiation, and food availability (Cheviron & Brumfield, 2012;Dillon et al., 2006;Montero-Mendieta et al., 2019). Overall, 514 outlier loci out of 13 455 SNPs were linked to altitudinal differences. Since no annotated genome is available for L. tityrus, we found putative functions in only 36 outlier loci, 14 out of which were associated with membranerelated proteins or functions. Similarly, five out of 11 outlier loci with putative functions were found to be membrane-related proteins in a previous study (Trense et al., 2021). Nevertheless, no common outlier SNPs were found in both the current and the previous study. The missing overlap of outlier SNPs may result from selection and other processes differing between the scales investigated, leading to different outliers being detected. In the previous study, we sampled two contiguous subvalleys, with accordingly little overall genetic divergence such that outliers may represent more recent selection events (Trense et al., 2021).
In the present study, we sampled several different sites in the European Alps, which are genetically clearly differentiated, thus representing a longer history of genetic differentiation. It is also likely that the outlier loci detected in the current study are associated with environmental variables that vary substantially at a broader scale. This may include selection on membrane features along the altitudinal gradient, probably in relation to membrane fluidity, which affects cold tolerance in ectotherms (Hazel, 1995;Hochachka & Somero, 2002). Higher proportions of unsaturated fatty acids in the membrane increase the membrane fluidity and thus cold tolerance (Brown et al., 2019;Haubert et al., 2008;Ohtsu et al., 1998). We were also able to show here that some outlier SNPs linked to altitude were also linked to temperature, suggesting thermal selection. Interestingly, five other outlier loci, the serine/threonine protein and phosphatidylinositol phosphatase/kinase, the zinc-finger protein, and other potassium channels and phenoloxidases were also detected in other studies on insects as potential outlier loci in relation to altitude Montero-Mendieta et al., 2019;Trense et al., 2021;Waldvogel et al., 2018). These outlier loci are involved in several cellular processes (Cassandri et al., 2017;MacKinnon, 2003;Wera & Hemmings, 1995).

| CON CLUS IONS
This study shows high genetic differentiation between L. t. subalpinus populations in the European Alps. Large-scale patterns were likely shaped by the last glacial period and the location of refuges in the western, central, and eastern Alps. At a smaller scale, high mountain ridges and large river valleys limit gene flow and thus structure populations. We suggest that conservation efforts should focus on the western Alps based on the high genetic diversity found there.
Furthermore, the dolomite populations could be treated as separate management unit. Our findings demonstrate the usefulness of genome-wide SNPs for estimating population structure, constraints on dispersal, and selection in the wild. In times of global climate change, it is important to better understand the population genetic structure of alpine species, because they seem to be particularly vulnerable to global warming (Engler et al., 2011;Hoffmann, 2010;Schmitt et al., 2014).

ACK N OWLED G M ENTS
We Furthermore, we thank the Bundesamt für Umwelt (BAFU) in Switzerland for the permission using genetic resources. We are grateful for the constructive comments provided by three anonymous reviewers. Open Access funding enabled and organized by ProjektDEAL.

CO N FLI C T O F I NTE R E S T
None declared.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found in the online version of the article at the publisher's website.