Anthropogenic and natural barriers affect genetic connectivity in an Alpine butterfly

Dispersal is a key biological process serving several functions including connectivity among populations. Habitat fragmentation caused by natural or anthropogenic structures may hamper dispersal, thereby disrupting genetic connectivity. Investigating factors affecting dispersal and gene flow is important in the current era of anthropogenic global change, as dispersal comprises a vital part of a species’ resilience to environmental change. Using finescale landscape genomics, we investigated gene flow and genetic structure of the Sooty Copper butterfly (Lycaena tityrus) in the Alpine Ötz valley system in Austria. We found surprisingly high levels of gene flow in L. tityrus across the region. Nevertheless, ravines, forests, and roads had effects on genetic structure, while rivers did not. The latter is surprising as roads and rivers have a similar width and run largely in parallel in our study area, pointing towards a higher impact of anthropogenic compared with natural linear structures. Additionally, we detected eleven loci potentially under thermal selection, including ones related to membranes, metabolism, and immune function. This study demonstrates the usefulness of molecular approaches in obtaining estimates of dispersal and population processes in the wild. Our results suggest that, despite high gene flow in the Alpine valley system investigated, L. tityrus nevertheless seems to be vulnerable to anthropogenically‐driven habitat fragmentation. With anthropogenic rather than natural linear structures affecting gene flow, this may have important consequences for the persistence of species such as the butterfly studied here in altered landscapes.

| 115 suitable habitat, increasing time, energy, and survival costs (Baguette et al., 2013;Bonte et al., 2012). Also, the ability to track shifting climatic niches is important in times of global climate change (Senner et al., 2018;Stevens et al., 2012). Consequently, dispersal comprises a major part of a species' resilience to environmental change (Senner et al., 2018), with poorly dispersing species predicted to suffer the most .
Constraints on dispersal may arise from intrinsic and extrinsic factors such as barriers. In human-modified landscapes, dispersal barriers include both anthropogenic and natural features. Several studies have found that natural barriers such as forests, rifts and rivers can structure populations through hampering dispersal (e.g., Heidinger et al., 2013;Muñoz-Mendoza et al., 2017;Omondi et al., 2019), while others found no such effects (e.g., Lada et al., 2008;Link et al., 2015).
Here, we examined the population genomic structure of the Sooty Copper butterfly Lycaena tityrus (Poda, 1761) in the European Alps, with the aim of understanding the effects of natural and anthropogenic barriers on this structure. Species of the genus Lycaena have (i) limited dispersal ability, as evidenced by mark-recapture studies; (ii) form closed population structures with individuals moving mainly within patches; and (iii) which may consequently result in strong genetic population structuring (Finger et al., 2009;Fischer et al., 1999;Fischer & Fiedler, 2001;Ricketts, 2001). Indeed, Lycaena species are among the least mobile butterfly species (Ricketts, 2001).
Lycaena tityrus is declining in large parts of central Europe Van Swaay & Warren, 1999), while the alpine subspecies L. tityrus subalpinus (Speyer, 1851) still occurs in high densities.
Alpine environments are characterized by an abundance of natural habitat with widespread natural barriers strongly structuring landscapes. Unlike in most other studies, this provides an opportunity to investigate gene flow in a largely intact landscape (Leidner & Haddad, 2010;Ugelvig et al., 2012).
Apart from considering dispersal, we also examine altitudinal adaptation, taking advantage of steep gradients over short distances in the Alpine habitat of this species. In L. tityrus, altitude has been linked with pronounced differences in several life history traits and with variation at the PGI locus (Karl et al., 2008). However, little is known about the molecular and physiological mechanisms underlying altitudinal variation in this and other species (Karl et al., 2008;Keller et al., 2013). Mechanisms may include proteins regulating membrane fluidity; for instance cold tolerance may be increased by changes in the ratio of saturated versus unsaturated fatty acids (Haubert et al., 2008;Ohtsu et al., 1998). Furthermore, cold adaptation in high habitats may be based on changes in various metabolic pathways (Cheviron & Brumfield, 2012).
Until recently, dispersal investigations have been hampered by the lack of appropriate methods to measure long-distance dispersal and gene flow, i.e., successful dispersal (Fountain et al., 2018;Kim & Sappington, 2013). Traditional approaches such as mark-release-recapture or individual tracking can be difficult, and can underestimate dispersal rates and fail to reveal past dispersal events (Sielezniew et al., 2011;Ugelvig et al., 2012). New molecular genetic techniques as used in landscape genetics allow dispersal to be analysed more precisely and provide additional information on gene flow (Fountain et al., 2018;Nève, 2009;Ugelvig et al., 2012), including estimates of relatedness that can capture dispersal events covering more than one generation (Jasper et al., 2019;Schmidt et al., 2018;Yang et al., 2020).
Here, we applied double digest restriction-site associated DNA sequencing (ddRADseq) to characterize single nucleotide polymorphisms (SNPs) for exploring population genomic structure in Alpine populations of L. tityrus. Specifically, we tested three hypotheses. (i) As the species is considered to be dispersal-limited, a strong population genetic structure is expected, with individuals caught in close vicinity being more closely related to each other than to individuals caught further apart (isolation by distance). (ii) Natural (forests, ravines, rivers) and anthropogenic barriers (roads) should constrain dispersal in this species and therefore affect genetic structure. (iii) Butterflies from different altitudes are expected to be adapted to local conditions, leading to the detection of outlier loci with links to relevant biochemical pathways.

| Study organism
We investigated genetic connectivity in the Sooty Copper (Lycaena tityrus, Alpine spp. subalpinus; Lepidoptera: Lycaenidae) in a valley system in the European Alps. Lycaena tityrus is a widespread temperate-zone butterfly 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 (Brunzel et al., 2008;Settele et al., 2008). The larvae feed mainly on Rumex acetosa L., but also utilize some congeneric plant species (e.g., R. acetosella L., R. scutatus L.; Ebert & Rennwald, 1991;Settele et al., 2008;Tolman & Lewington, 2008). The species has 1-3 generations per year and hibernates as half-grown larvae. The alpine subspecies L. t. subalpinus is confined to the higher altitudes of the Alps and some other European mountain ranges, where it is relatively widespread. It has only one generation a year (Tolman & Lewington, 2008). Its altitudinal distribution ranges from 1200 to 2500 m a.s.l. (Tolman & Lewington, 2008).

| Population sampling and study area
The study area is situated in the Central Alps, Austria, spanning an altitude from ca. 1200 to 2200 m a.s.l. (Figure 1). Specifically, we sampled male butterflies along the Ötz valley, which is one of the longest lateral valleys in the Austrian Central Alps. It starts near the village Ötz at an altitude of 785 m a.s.l. After around 35 km, near Zwieselstein (1470 m a.s.l.), the valley splits into two tributaries, the Gurgl and the Vent valley. Both tributary valleys are separated by mountains up to 3500 m a.s.l. . We will hereafter refer to the lowest section until Zwieselstein as "Ötz", and to the two tributaries as "Gurgl" and "Vent". Sampling sites ranged from 1259 to 2074 m a.s.l., the maximum distance between sites was 19.5 km.
Furthermore, we sampled five more isolated sites (hereafter "outside groups" to avoid confusion with phylogenetic outgroups) in the surroundings of the Ötz valley system (Figure 1b)

| ddRAD library preparation
From each individual, we used the thorax and head for the extraction of genomic DNA with the DNeasy Blood & Tissue Kit (Qiagen, Germany). We followed the manufacturer's instructions but including an additional step of RNase A treatment. We applied double digest restriction-site associated DNA sequencing (ddRADseq) following Rašić et al. (2014) with some modifications. We used 120 ng DNA per individual for digestion in a 40 µl reaction solution with 0.8 µl (10 units) each of NlaIII and MluCI restriction enzymes, 4 µl CutSmart Buffer (New England BioLabs Inc, USA), and nuclease-free water. Digestion reactions were incubated at 37°C for 3 h. After digestion, the products were cleaned with 60 µl of Sera-Mag beads F I G U R E 1 Location of (a) the study area, indicated by a black square, within Europe, (b) the Ötz valley system (black square) and the five outside groups (black dots) within the Austrian and Italian Alps, and (c) the sampling locations within the Ötz valley system on an aerial map. The yellow stars (Ötz), red triangles (Vent) and green squares (Gurgl) show the location of individuals within the three subvalleys, and the white dot illustrates the village Zwieselstein. The orange lines indicate ravines or steep wooded slopes, putatively comprising barriers to dispersal and thus resulting in nine spatial sections. The blue lines indicate rivers and the black lines roads. The underlying river and road networks were taken from the Austrian's Roads and Waterways shapefiles, available on MapCruzin.com [Colour figure can be viewed at wileyonlinelibrary.com] (1.5× volume; GE Healthcare and Thermofisher, Australia), and ligated to specific combinations of P1 and P2 adapters. Sixteen P1 and four P2 adapters were used for ligation in a 45 µl reaction solution, including 35 µl cleaned DNA (1 ng/µl), 4.5 µl 10× T4 Buffer, 0.55 µl (1000 units) T4 DNA ligase (New England BioLabs Inc, USA), nuclease-free water, and 2 µl each of P1 and P2 adapters. The ligation reactions were incubated at 16°C for 12 h, then heat-inactivated at 65°C for 10 min, cooled down to 22°C at a rate of 2°C per 90 s, and held for 5 min at 22°C.
We pooled individuals into three libraries, each containing DNA fragments from individuals with different adapter pairs.
These were pooled and cleaned with 1.5× volume of Sera-Mag beads. Size selection of 250-400 base pairs (bp) fragments was performed using a 2% gel cassette, and analysed by pippin-prep

| SNP calling
SNPs were called using a de novo pipeline without a reference genome. We used the process_radtags program within stacks version 2.3 (Catchen et al., 2011(Catchen et al., , 2013 to demultiplex the sequence data. This separated the individuals by their unique barcode, removed adapter sequences, trimmed the reads to 120 bp length, and discarded low-quality reads (Phred score < 20). We obtained an average of ~6,163,048 reads per individual, which covered 2,003,979 loci across 26,411 RADtags. We then ran the stacks denovo_map pipeline using one individual per outside group and 26 individuals from all over the Ötz valley system to generate a catalogue of loci (treating all individuals within the Ötz valley system as one population and each outside group as separate population). Each of the 186 individuals was then matched against this catalogue and SNPs were subsequently called. We filtered the data to retain SNPs which were present in more than 50% of the 186 individuals, and then retained the first SNP in each ddRAD locus. 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 mean depth of 20-45. For filtering, all individuals were treated as belonging to one population. The final data set contained 12,806 SNPs for analyses of relatedness, genetic structure, and outlier loci (average read depth was 30.86).

| Analysis of relatedness
To avoid biased population structure analyses, we calculated pairwise relatedness among individuals using the Triadic Likelihood Estimator (TrioML, which gives a relatedness coefficient) in coancestry version 1.0.1.9 (Wang, 2007(Wang, , 2011. Individuals with a TrioML estimator of ≥ 0.250 are probably full-siblings, between 0.125 and 0.249 probably half-siblings, and between 0.060 and 0.124 probably cousins. Note that TrioML can only assess genetic relatedness but not the real ancestry. We found several related individuals including one probable pair of full-siblings, of which one individual was excluded from downstream analyses. We additionally used the Specific Hypothesis Test function with 9,999 permutations implemented in the program ml-relate (Kalinowski et al., 2006) to corroborate putative full-and half-siblings, although this program cannot be used to test cousin relationships.

| Isolation by distance, resistance, and environment
To test for isolation by distance (IBD)and by landscape features, we calculated pairwise genetic distances on an individual basis according to the Bray Curtis approach, using the package vegan (Oksanen et al., 2019) in r version 3.5.2 (Bray & Curtis, 1957;R Core Team, 2015). Pairwise, individual-based geographic distances were calculated as (i) Euclidean distances (with and without outside groups) with the function PointDistance (plane) in the package raster (Hijmans & van Etten, 2016) in r, based on the WGS84[EPSG:32632]/ UTM zone 32N; (ii) distance along rivers/valley floors as likely dispersal routes (only without outside groups); and (iii) altitudinal distances (with and without outside groups). We tested for correlations between genetic and geographic distances (isolation by distance) using Mantel tests (using the Mantel function from the vegan package in r; Mantel, 1967) with 9,999 permutations.
To investigate the spatial arrangement of dispersal corridors and barriers, we additionally used an estimated effective migration surface (eems; Petkova et al., 2016) algorithm. This analysis shows the spatial patterns of genetic diversity and gene flow, and thus visualizes geographic regions that deviate from isolation by distance due to increased or reduced gene flow (please note that this analysis is not testing the influence of landscape variables). We tested the deme numbers 200 and 1,000, using RunEEMS_SNPS in eems with three independent runs of 200,000 Markov Chain Monte-Carlo (MCMC) iterations, each with a burnin of 100,000 and a thinning interval of 9,999.
To subsequently analyse isolation by resistance (IBR) and by environment (IBE, see below), we first modelled the distribution of L. tityrus within the Ötz valley system with maxent version 3.4.1 (Phillips et al., 2017) including bioclimatic variables as well as landscape features. To avoid effects of multicollinearity (|r| > .7), we performed pairwise Pearson correlations for the 19 bioclimatic variables provided by Worldclim (worldclim.org). Accordingly, we only used three bioclimatic variables, namely annual mean temperature, annual precipitation, and precipitation seasonality.
The landscape features included were altitude, water bodies, tree cover, grass cover, land cover, aspect (i.e., how much sunlight a surface receives), hillshade (intensity of light on a surface), and slope.
The raster file data for altitude were taken from Natural Earth (naturalearthdata.com), for water bodies, tree cover, and grass cover from Copernicus (land.copernicus.eu), and for land cover from the European Space Agency (ESA) GlobCover (due.esrin.esa.int). We used the digital surface model from the Open Data Portal Austria (data.opendataportal.at) to generate aspect, hillshade, and slope using the function Surface Analysis in arcgis version 10.3.1 (ESRI).
To combine bioclimatic variables and landscape features, which need to be identical in resolution, extent, and projection, we used the packages raster and rgdal in r. For maxent distribution modelling, we selected default settings but changed the random test percentage to 25. We performed jackknife analyses to determine the landscape features which reduce the model reliability when omitted. For subsequently analysing IBE, we calculated the matrices of pairwise least cost path distances based on the reverse suitability values of the maxent estimates for potentially isolating features, i.e., altitude, water bodies, and tree cover, using the function costdistance in arcgis (ESRI). For IBR, we used the reverse suitability values of maxent estimates for the isolating features, i.e., altitude, water bodies and tree cover and ran circuitscape version 4. 0 (McRae et al., 2013) to receive resistance surfaces and distances for these features. Afterwards, we followed Cushman et al. (2006) to analyse 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. For each landscape feature, we calculated eight partial Mantel tests to analyse the relationships between genetic and landscape distance matrices, partialling out the effect of other landscape distance matrices (IBD, IBR, IBE).
To account for the nonindependency of records from one individual, we additionally run maximum-likelihood population-effects (MLPE) following Clarke et al. (2002) and Peterman (2018). As the MLPE models for isolation by environment and resistance showed that most of the variation in genetic distances was found among individuals (ca. 50%) while environmental factors explained < 0.01% throughout, these results are not presented.

| Analysis of molecular variance (AMOVA)
We used an hierarchical AMOVA in arlequin version 3.5 (Excoffier & Lischer, 2010) to assess the genetic variance among (i) the three valleys (Ötz, Vent, Gurgl); (ii) groups being isolated by potential dispersal barriers, namely ravines separating the valley into different sections or populations located on plateaus above the valley separated by steep forested slopes (ravines/isolated groups; Figure 1c); (iii) among individuals; and (iv) within individuals. We used a priori grouping to test how well those putative clusters are supported.

| Analysis of genetic structure
For identifying genetic clusters, we performed a structure analysis within the Ötz valley system using the function snmf in the package lea (Frichot & François, 2015) in r and in structure (Pritchard et al., 2000). These analyses revealed that all individuals belonged to one cluster (results not shown). Therefore, we grouped the individuals collected within the Ötz valley system according to different factors which may affect genetic connectivity. The factors considered were (i) sample site (i.e., groups of 4-7 individuals, which were caught in close proximity, typically on the same meadow; mean distance 0.42 km); (ii) altitude (m a.s.l., here we used the altitude of each individual, without grouping them); (iii) ravine/isolated group, and side of the (iv) river; and (v) road within the Ötz, Vent, and Gurgl valleys. The river and road networks were obtained from the Austrian's Roads and Waterways shapefiles in MapCruzin.com.
For analysing the influence of the above factors, we calculated pairwise genetic distances on an individual basis (Bray & Curtis, 1957) to run distance-based redundancy analyses (dbRDA; Legendre & Anderson, 1999) with the r package vegan (Oksanen et al., 2007). dbRDAs were run for each factor separately owing to strong intercorrelations (Spearman correlation coefficients between .26 and .81), and p-values were corrected following Benjamini-Hochberg (Benjamini & Hochberg, 1995). We used the functions Capscale and Vif.cca for dbRDA as well as the ANOVA function from the package stats in r.
We used the packages adegenet (Jombart, 2008;Jombart & Ahmed, 2011) and hierfstat (Goudet, 2005) in r to calculate observed heterozygosity, expected heterozygosity, inbreeding coefficient, and global F ST for the different valleys, altitudes, ravines/isolated groups, sides of rivers and roads, and sample sites.

| Outlier loci analysis
To identify potential outlier loci at different altitudes, we clustered v. 1.1 using samtools version 1.2 (Li, 2011;Li et al., 2009). To further characterize the loci under selection, we used blast2go version 5.2.5 (Götz et al., 2008) for functional analysis of the identified proteins.

| Isolation by distance, resistance, and environment
Genetic distances among all individuals (Bray-Curtis) ranged from 0.122 to 0.161, and the geographic (Euclidean) distance within the the subvalleys was not significant throughout. The eems analysis showed isolation by distance for individuals from Ötz and Gurgl, and a moderate migration rate for individuals in Vent ( Figure 5).
However, no major dispersal barrier was found within the Ötz valley system.
The receiver operating characteristic (ROC) analysis in the maxent model for the distribution of L. tityrus indicated a high model accuracy with an AUC value of 0.992. Annual precipitation showed the highest contribution to the model fit with 79.6%, followed by land cover with 6.6%, annual mean temperature with 6.5%, and tree cover with 2% ( Figure 6). The results of the concomitant IBE and IBR analyses showed that IBE (altitude and tree cover) affects genetic structure in the Ötz valley system, while IBR was not significant throughout (Table 1).

| AMOVA
AMOVA revealed that the highest amount of genetic variation occurred within individuals (Table 2), while much less structure existed among individuals, ravines/isolated groups, and the three valleys.
Nevertheless, genetic structure was statistically significant at each hierarchical level.

| Spatial genetic structure
Sample site, altitude, ravines/isolated groups, side of road but not side of river had a significant impact on genetic structure within the Ötz valley system (Table 3). For basic population genetic information and F ST values for the different groups tested see Appendix Tables S1 and S2. In general, individuals from each of the 25 different sample sites were clustered to some extent, indicating that the individuals caught within sample sites tended to be genetically more similar than those from other samples (Figure 7a). This suggests some structure despite the low genetic distances among samples. The sample site depicted by black triangles in Figure 7a, which includes part of the individuals of section 1 in Figure 1 (i.e., individuals from the lowest altitude), showed the strongest divergence from other groups, followed by the sample site visualized by brown circles, which is identical with section 8 in Figure 1 and thus involves individuals from a plateau. The nine sections defined by putative barriers to dispersal (ravines/isolated groups) were clearly genetically different, with overlap occurring only among sections 2, 3, and 4, which all belong to the Venter valley (Figure 7b, cf. Figure 1). Interestingly, the were most strongly differentiated from all others.

| Outlier loci according to altitude
In total, 312 SNPs were identified as potential outliers according to altitude, being found in all three outlier tests. baypass, FDIST2, and H-FDIST2 detected 639, 381, and 382 significant outlier SNPs, respectively ( Figure S1). Fourteen out of the 312 putative outlier SNPs had a match in blast2go (Table 4). Of these, five could be assigned to membrane proteins, two to transcription, and one each to RNA bind-

| Population genetic structure
Our results point to high levels of dispersal and gene flow within the Note: Factors are sample site (i.e., groups of 4-7 individuals caught in close proximity), altitude, ravines/isolated groups (i.e., groups isolated by putative dispersal barriers), left and right side of river and road, respectively. All variables were analysed separately. Effect sizes are given as partial Eta squared (η 2 ), and significant p-values are given in bold. Accounting for multiple testing using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995) confirmed all significant factors.
genetic variance among valleys and ravines/isolated groups (AMOVA results); (v) no structuring effect of rivers; and (vi) surprisingly large geographic distances between putatively related individuals. Regarding the latter, we found two putative pairs of cousins with one individual within the Ötz valley system and the other within an outside group.
This suggests that at least some individuals are highly mobile and have the potential to disperse over long distances, although Lycaena species are among the least mobile butterfly species (Ricketts, 2001). Note that we exclusively sampled males, while females are often the more dispersive sex in butterflies (Fischer et al., 1999;Hanski et al., 2004;Reim et al., 2018). This, however, should not affect results based on genetic relatedness in a sexually reproducing organism.  The relatively high levels of gene flow we have detected in a species formerly considered to be sedentary highlights the power of high density molecular markers to provide new insights into dispersal in the wild (Escoda et al., 2017;Maas et al., 2018;Schmidt et al., 2018). In contrast to genetic methods, direct estimates of dispersal such as obtained through mark-release-recapture studies systematically underestimate dispersal due to the difficulty of detecting rare long-distance movements (Sielezniew et al., 2011;Ugelvig et al., 2012). Genetic approaches do not suffer from such constraints, and moreover provide estimates of realised dispersal based on individuals that successfully arrive and reproduce in new patches (Fountain et al., 2018;Nève, 2009;Sielezniew et al., 2011;Ugelvig et al., 2012).
Despite high levels of gene flow and only weak evidence for isolation by distance within the Ötz valley system, genetic structuring was observed across sample sites, ravines/isolated groups, roads, and altitudes. The strongest genetic differentiation between sample sites involved individuals from the lowest altitudes (black triangles, forests surrounding habitat patches (Crawford & Keyghobadi, 2018;Haddad, 1999;Heidinger et al., 2013;Schmitt et al., 2000;Tewksbury et al., 2002). Likewise, human settlements may hamper dispersal (Miles et al., 2019;Vandergast et al., 2009;Yang et al., 2011). Both factors may have contributed to genetic differentiation in the Ötz valley system, as is evidenced by our maxent analyses.
Overall, natural features (altitude, precipitation, ravines, and forests) seemed to have a larger impact on the genetic structure of L. tityrus than anthropogenic ones. The evidence of existing IBE in altitude and tree cover within the Ötz valley system needs to be interpreted with caution because simple and partial Mantel tests suffer from high type I error rates (Cushman et al., 2013;Guillot & Rousset, 2013). Ravines essentially provide no suitable habitat for L. tityrus, and separated different genetic groups except within the Vent valley. In this area, two relevant ravines are shorter and shallower than others, providing at least some adequate habitat for the butterfly. The occurrence of such "stepping stones" may explain higher levels of gene flow in Vent. In some species, stepping stones may increase dispersal rates and maintain population connectivity and gene flow (Gilpin, 1980;Saura et al., 2014;Ugelvig et al., 2012).
In general, both species-specific dispersal ability and the spe-   (Douwes, 1976;Mattila, 2015). For instance, surface temperatures may reach 55°C on roads, while the maximum temperature of Alpine rivers ranges between 10 and 20°C (Linden & Petty, 2008;Webb & Nobilis, 1995). Perhaps behavioural factors are involved, in that butterflies are less reluctant to cross familiar natural structures than novel roads. This speculation generates testable hypotheses with regard to edge-following and edge-crossing behaviour (Fernández et al., 2016;Nguyen & Nansen, 2018), which could be considered in future research. An alternative explanation is that traffic and associated road kills rather than roads per se affect gene flow across roads. Several studies have revealed that roads may comprise major dispersal barriers for insects (Bhattacharya et al., 2003;Keller & Largiadèr, 2003;Miles et al., 2019;Muñoz et al., 2015;Schmidt et al., 2018). Small insects were found to be overrepresented in road kills (Muñoz et al., 2015;Skórka et al., 2013;Soluk et al., 2011). We speculate that road kill is a major threat to dispersing butterflies in our study area, evidenced by dead individuals that can be found along the roads (own observation). Either way, anthropogenic habitat fragmentation through roads seems to be a greater threat than natural rivers to genetic connectedness for L. tityrus in the Ötz valley system.

| Altitudinal adaptation
Overall, 312 outlier loci out of the 12,806 SNPs scored were linked to altitudinal differences. This is relatively low compared to similar studies on other species (cf. Jordan et al., 2017;Sarkar et al., 2019;Woodings et al., 2018), which may reflect limited adaptation in the Ötz valley system but also relatively high gene flow making weak selection difficult to detect. As no annotated genome for L. tityrus exists, only eleven out of the 312 outlier loci could be assigned a putative function. The detection of five membrane-related proteins may suggest selection on membrane properties along the altitudinal gradient, potentially caused by constraints on membrane fluidity.
Note that the fluidity of membranes may affect cold tolerance in ectotherms (Hazel, 1995;Hochachka & Somero, 2002). Unsaturated fatty acids in the membrane can increase membrane fluidity, causing higher cold tolerance (Brown et al., 2019;Haubert et al., 2008;Ohtsu et al., 1998). Furthermore, in two studies on honey bees and bumble bees (Jackson et al., 2020;Montero-Mendieta et al., 2019), other potassium channels were detected as potential outlier loci in relation to altitude. Four out of the remaining six loci with known function encode housekeeping proteins (Table 4). Another protein, Kinesinlike protein, is associated with ATP binding (Wang et al., 1996), which may indicate climatic adaptations related to metabolism. Indeed, metabolic cold adaptation, with metabolism increasing under thermally challenging conditions, is a widespread feature in ectotherms (Addo-Bediako et al., 2002;Shik et al., 2019). The last protein, prophenoloxidase, has a gene product involved in oxidation reduction.
It inhibits proteins and recognises bacterial and fungal components, and is therefore an important part of the insect immune defence (Söderhäll & Cerenius, 1998). Prophenoloxidase was also found to have an impact on local adaptation to different altitudes (González-Tokman et al., 2019;Karl et al., 2010).
In conclusion, our study showed surprisingly high levels of gene flow in L. tityrus within an Alpine valley system. Nevertheless, we found specific landscape features that limited gene flow, namely forests, ravines, and roads. These findings highlight the power of genome-wide SNPs for estimating realised dispersal and population structure in the wild. Interestingly, linear anthropogenic barriers (roads) served as genetic barriers, but linear natural barriers (rivers) did not. We speculate that this might be due to individuals being reluctant to disperse across novel, anthropogenic features.
Consequently, L. tityrus seems to be vulnerable to anthropogenic habitat fragmentation. This may, in addition to habitat loss, explain its regional declines in Central Europe Van Swaay & Warren, 1999). Additionally, we found evidence for selection on loci related to membranes, metabolism, and immune function along altitude. The functions involved suggest temperature as a selective agent, which is expected to be closely related to altitudinal variation.

ACK N OWLED G EM ENTS
Many thanks to Nancy Endersby-Harshman for helpful advice in the laboratory. This research was facilitated by use of the Nectar Research Cloud, a collaborative Australian research platform supported by the Australian Research Data Commons (ARDC) and