This study was conducted in the Swiss National Park (SNP) in the Central Alps. The SNP is located at the south-eastern edge of Switzerland. Established in 1914, the SNP is aiming for a strict preservation of natural processes. It has a surface area of 170.3 km², including various habitats in an elevational range from 1’380 m a.s.l. up to 3’173 m a.s.l. (Schlüchter et al. 2021). Adjacent areas, some of which are difficult to access due to their location in high alpine terrain, are used extensively for livestock farming and forestry. Compared to other regions in the Alps, the climate in the region of the SNP is dry and continental with a pronounced annual cycle. In winter, there is a continuous blanket of snow for months. Due to the alpine terrain, temperatures are mainly determined by elevation and exposition (Haller et al. 2014; Gubler and Robinson 2025). In Buffalora at the eastern edge of the SNP (1’971 m a.s.l.), mean annual temperature (MAT) is 1.1 °C with a mean annual precipitation of 936 mm (MeteoSwiss 2025). Since 1917, MAT has increased here by 2.1 °C with warming accelerating since the late 1980s, while mean annual precipitation has been decreasing since 2010 after a long period of stability (Gubler and Robinson 2025).

Our study sites are located inside or in the immediate surroundings of the SNP. The butterfly surveys conducted in 1998, 2001, and 2004 in the SNP region used 20 study plots with an area of 0.25 ha each which are scattered across the SNP on elevational gradients at three distinct sites (Pasche et al. 2007; Gonseth et al. 2013). The first six plots surveyed in 1998 are in the catchment area of the Pass dal Fuorn at the center, eight plots surveyed in 2001 are in Val Trupchun to the south, and the last six plots surveyed in 2004 are in Val Mingér to the north of the SNP (Fig. 1). Most of the plots used in this study are located on dolomite bedrock, except in Val Trupchun, which is located on marl bedrock (Haller et al. 2014; Swisstopo 2025). The plots are either south- (Val Trupchun, Pass dal Fuorn) or north-exposed (Val Minger, Pass dal Fuorn). Even though the plots are predominantly covered by alpine pastures and grasslands, the immediate surroundings of the plots vary and encompass a wide range of different habitats. Typically, plots below the forest line are located at small to medium sized clearings, often in close proximity to mountain streams and the gravel banks formed by them. Different types of alpine grasslands, scree fields, avalanche tunnels and bare rock comprise the surroundings of the plots above the forest line. Most of the plots are located in areas with a slight to moderate slope, but the surrounding terrain can be very steep (Pasche 2005; Haller et al. 2014). The plots are not distributed entirely evenly along the covered elevational gradient from 1.672 m a.s.l. up to 2.431 m a.s.l., with a higher density in the middle to lower parts of the range (Fig. 2).

Figure 1. 

Butterfly research in and around the Swiss National Park: map of the study area in the Engadin, Canton of Grisons.

Figure 2. 

Distribution of sampling sites over elevation within the overlapping elevational range. The histograms of both sampling periods are overlaid and visualise the different survey designs, namely the plot-based approach of the old data and the continuous transects of the recent data. For the recent data, the distribution of transect metres indicates how many metres of both transects are within a vertical section of 25 metres.

In 2025, two transects were established, both following hiking trails: one of them (Transect 1) is in the valley Val dal Botsch at the centre of the national park (Fig. 1). The majority of the transect follows the south-facing valley and climbs its south-exposed slopes on porous dolomite (Swisstopo 2025). Transect 1 covers an elevational range of 883 metres from 1’794 m a.s.l. to 2’677 m a.s.l. on 5.8 kilometres horizontal distance. The second transect (Transect 2) is located further west at the high alpine lake plateau of Macun, an exclave of the national park (Fig. 1), and the adjoining valley Val Zeznina. This second transect covers an elevational range of 1’125 m from 1’730 m a.s.l. up to 2’855 m a.s.l. on 7.6 kilometres in north- and north-west facing slopes. In contrast to transect 1 in Val dal Botsch, transect 2 in Val Zeznina and on Macun is located on crystalline and thus silicate bedrock such as gneisses and amphibolites (Gubler and Robinson 2025; Swisstopo 2025). The transects surveyed in 2025 cover a wide range of habitats and slopes. Below the forest line, the transects transit a mosaic of subalpine forests and clearings interspersed with gravel banks and scree fields. Above the tree line, the surrounding of the transects is dominated by alpine grasslands, but scree fields, bare rock and avalanche tunnels are frequent habitat features. Although the variation and number of habitats is exceeding the habitat diversity in the survey plots of the old data, conditions are equivalent across the larger part of the transects. The transects cover the elevational gradient evenly. Only the lowest part of the elevation range has a disproportionately high coverage due to the gentler slope within this area (Fig. 2).

For the old data, within each plot, a transect of 320 or 350 m was established, where butterflies were counted. Butterflies passing within a semicircle with a radius of 5 metres were recorded. Burnet moths were neglected. The transects were walked at a steady pace at a speed of approximately 2 to 3 km/h. Butterfly surveys were only conducted between 11 am and 14 pm if weather conditions were suitable for butterfly activity. Each plot was surveyed five to six times per season, while every site was examined in one season only. The plots in the Fuorn- catchment were investigated in 1998, the plots in Val Trupchun in 2001 and the plots in Val Mingér in 2004. Each transect was surveyed five to six times from mid-June to early September at regular intervals (Pasche 2005).

For the transect survey conducted in 2025, the methodology was adapted to the continuous transects over the elevational gradient and more closely aligned with the methodology used in the Swiss Biodiversity Monitoring (BDM, Biodiversitäts-Monitoring Schweiz 2021). In the three months from June 10^th to September 10^th, each transect was surveyed twelve times in a weekly interval (3 to 11 days in between each transect walk, average 6.5 days). As in the old survey, all butterflies (Lepidoptera: Papilionoidea) were counted in a distance of 5 metres. For the transect sections on forest roads, the semicircle was split up, and butterflies were counted in 5 metres distance from the edges of the road. The target speed for walking the transects was adopted from the methodology used in the old data, but the time interval was extended due to the length of the transects. Thus, transects were walked between 10 am and 17 pm (10.30 am to 17 pm after August 10^th). Weather conditions considered as suitable for transect surveys were specified based on the guidelines given by the BDM (Biodiversitäts-Monitoring Schweiz 2021). The observations were recorded using the Webfauna app (Gerber et al. 2025), which located all individuals precisely via GPS coordinates. To avoid biases by diurnal activity patterns, the start time of the transect walks was varied and the surveys were started in different sections: at the upper end, the lower end and the middle of the transect.

Whenever possible, butterflies were identified at the species level by visual inspection. If this was not possible due to distance or butterflies were only flying past, these individuals were caught with a net for identification. In cases of doubt, the identification key given by Baudraz et al. (2020) was consulted. Grizzled skippers (Hesperiidae: Pyrgus sp.) that could not be identified based on external characteristics and individuals out of the Melitaea athalia – complex were collected and identified by genital dissection. The species complexes Leptidea sinapis – juvernica and Colias alfacariensis – hyale were not identified at the species level. Heavily worn individuals that could no longer be reliably identified were not included. The nomenclature was adopted from Wiemers et al. (2018).

The habitat-data used in this study was drawn from the HABITALP dataset (Lotz 2006). This dataset is based on the analysis of high-resolution colour infrared aerial images taken in August 2000. Habitat polygons were delineated and classified manually at a scale of 1:5000 resulting in an accurate representation of on-site conditions and habitat transitions (Lotz 2006). This high level of precision is the reason why the HABITALP dataset was preferred to more recent datasets with automated interpretation of aerial images, despite its age. Designated habitat types and degrees of coverage per soil coverage layers were extracted from the dataset (Suppl. material 1). It was assumed that changes to habitat characteristics caused by wind throws and slope movements would only slightly alter the shapes of habitat polygons and the overall classification. Natural succession was considered to affect all terrestrial habitat and land coverage types equally. Consequently, it was expected not to influence the subsequent analyses.

The following calculations were conducted with R version 4.3.2 (2023, https://www.r-project.org). To link habitat data with recent butterfly data, both transects surveyed in 2025 were split up in sections of 200 m distance and around these sections buffer areas of 100 m were established. This choice was based on observations that during capture-recapture surveys, the majority of individuals of examined alpine butterfly species was recaptured within a radius of 100 m from their initial capture (Junker et al. 2010; Polic et al. 2014; Ehl et al. 2019). For each buffer area, the centroid and its elevation were determined. The proportions of the different habitat types and degrees of coverage were extracted from the habitat data for each buffer area using the SF-package (Pebesma et al. 2025). Additionally, for each species, buffer areas with a 100 m radius were established around the coordinates of all observations. These buffer areas were merged, and area proportions were calculated per habitat type while degrees of coverage were neglected. As hiking paths were evenly distributed across the sections and the distribution of roads and buildings was uneven between the transects, these anthropogenic landscape features were excluded from the analyses to avoid ecologically nonsensical correlations.

The temperature adaptation of butterflies was represented using the Species Temperature Index (STI). This index is based on the average MAT across the geographic range per species and represents a relative measure of temperature adaptation of butterflies (Schweiger et al. 2014). Because STI is drawn from a large geographical range and does neither differentiate between lowland and alpine subspecies nor account for local conditions in the Central Alps, a more specific Swiss Species Temperature Index (SwSTI) was calculated in the same manner as the Europe-wide. All calculations were conducted with R version 4.3.2 (2023, https://www.r-project.org). BDM butterfly data (Weber et al. 2004) from 2013 to 2024 from all plots within the four alpine biogeographic regions of Switzerland (Northern flank of the Alps, Eastern Central Alps, Western Central Alps, Southern flank of the Alps)(Swisstopo 2025) was used. Species that were observed at less than ten plots or had less than 30 observations were excluded. Each observation within this dataset was assigned the MAT value of the nearest datapoint in a raster of MAT values of Switzerland. The raster data at a 25 × 25 m scale cover the period from 1981 to 2010 (Külling and Adde 2024). To match butterfly and temperature data, the TERRA-package (Hijmans et al. 2025) was used. The SwSTI was calculated as the mean of the MAT values of all observations per species. The fact that the periods when these two datasets were created do not overlap was disregarded, as SwSTI is a relative index.

All statistical analyses were conducted with R version 4.3.2 (2023, https://www.r-project.org). To disentangle the influence of habitat composition and elevation on the composition of the butterfly communities, a Variance Partitioning was performed. Here, only the recent butterfly data was used to take advantage of the continuous elevational gradient and to avoid possible distortion of the elevation effect by combining two different datasets. The variance partitioning was performed using the VEGAN-package (Oksanen et al. 2025). Habitat parameters and the elevation of the centroids of the buffer areas served as predictors for the composition of butterfly communities within the 200 m sections. The effects of both predictors were tested for robustness by redundancy analyses followed by ANOVA. As a dependency between habitat parameters and elevation was assumed, the habitat parameters were examined in the same way using redundancy analysis to determine their correlation with elevation.

Changes in the temperature adaptation of butterfly communities were assessed by the comparison of the Community Temperature Index (CTI) between old and recent data. For both datasets, every recorded individual was assigned the STI value corresponding to its species. Because the old data only covered the elevational range from 1’672 m to 2’431 m a.s.l., the recent data were truncated at 2450 m a.s.l. and data from below were used for the temporal comparison. The Community Temperature Index (CTI) was used as an abundance weighted mean value of STI values. Both a Generalized Linear Model (GLM) and a Generalized Additive Model (GAM) were fitted to test for differences between survey periods and to analyse the relationship between CTI and elevation. The GLM was built using the LME4-package (Bates et al. 2025) with elevation and sampling period as predictors for CTI values including the interaction term between the two factors. A Gaussian family was used. The GAM was performed using the MGCV-package (Wood 2025). For CTI as dependent variable, a smooth term over elevation per year (k = 3) and the year of the survey were included as explanatory variables. k-value was selected to avoid oversmoothing due to too low values and overfitting due to too high values taking model diagnostics into account. A Gaussian family and an identity link function was used. Both models were compared based on the Akaike Information Criterion (AIC) and the model with the lower AIC-value was selected.

The analysis performed with the Europe-wide CTI was repeated in the same manner for the SwCTI using SwSTI instead of STI. Two species were excluded from the analysis because the data used for the calculation of SwSTI was deemed insufficient (Suppl. material 2). The models were built the same way as for CTI and the model with the lower AIC-value was selected.

To detect changes in species’ elevational distribution, species with at least ten observed individuals within the overlapping elevational range in both survey periods were selected. Afterwards, the mean elevation was calculated per species and survey within the overlapping range. The two survey periods were compared by a t-test paired per species. Because of the plot-based approach used in the old data, investigating changes in leading and rear edges over time was not possible.

Temperature indices were linked to the upward displacement per species to examine differences between cold-adapted and warm-adapted species. The difference of the means between the survey periods was calculated per species and checked for a correlation with STI and SwSTI. Therefore, a linear model was built using STI and SwSTI as predictors for the differences in mean elevation.

The width of the habitat range used per species on a local scale was depicted by the Shannon diversity index of habitat types present in the surrounding of the observations within the recent data. The proportional habitat data from within the 100 m radius buffer zones around the observation points were used to calculate a Shannon diversity index of habitats per species. The index values were correlated with the difference in mean elevation by fitting a linear model.

Overall, 6030 individuals of 82 species were counted during the transect survey in 2025. Along transect 1, 58 species and along transect 2, 63 species were recorded. 39 species (47,6%) were found along both transects. In contrast, there are 14’071 recordings from 84 species in the old data. Here, the overlap between the three sites was 30 species (35,7%). Species in both datasets match by 74,2%. Species missing in 2025 were mostly found in Val Trupchun (10 out of 15) and species newly discovered in 2025 are predominantly found on transect 2 (7 out of 9).

Testing habitat and elevation for their influence on butterfly community composition, habitat was found to explain almost half of the variance (adj. R² = 0.464; F = 3.32; p < 0.001). The independent impact of elevation however appeared to be comparatively low (adj. R² = 0.015; F = 8.748; p < 0.001). The explained variance that could not be assigned to either habitat composition or elevation is summarized in the shared effect. The shared effect of elevation and habitat, meaning the availability of specific habitats at a specific elevation, also explained a considerable amount of the variance in butterfly communities (adj. R² = 0.101). Nevertheless, a considerable amount of the variance remained unexplained (adj. R² = 0.42). Yet habitat and elevation are not independent: habitat is for more than a quarter of the variance explained by elevation (Inertia = 741.2; Proportion = 0.312; adj. R2 = 0.3; F = 26.331; p < 0.001).

For evaluating the potential changes of CTI over time and elevation, the GAM was chosen based on AIC values (ΔAIC = 119,35). The GAM built for CTI revealed a clear negative contribution of elevation to STI values with index values decreasing at increasing elevation (1998–2004: edf = 1.991; F = 537.8; p < 0.001; 2025: edf = 1.921; F = 175.2; p < 0.001). This effect was non-linear and well differentiated between the two survey periods. While the slope decreased with increasing elevation in 1998–2004, it increased with elevation in 2025. On average, the CTI values in 2025 were higher than they had been within the old data (Estimate = 0.502; t = 17.21; Pr(>|t|) < 0.001). The model explains only a small part of the deviance (adj. R² = 0.074). Projecting the model on the data, the difference of CTI between sampling periods is particularly evident at medium elevations (1’900–2’300 m a.s.l.), while towards the edges of the covered range, the index values are converging (Fig. 3).

Figure 3. 

Visualisation of the model predictions for the Generalized Additive Model (k = 3) for Community Temperature Index representing temperature adaptation of butterflies in the Swiss National Park over elevation. Old (1998-2004) and recent data (2025) are compared. Error ribbons show standard error.

For the equivalent investigation conducted using SwCTI, the GLM was found to be sufficient based on AIC values (ΔAIC = 0.26). Nevertheless, results look very similar: Here too there is a clear negative effect of elevation to the Temperature Index (estimate: -0.00318, std error: 0.000006458, t value: -48.272, p < 0.001). The decrease of SwCTI with increasing elevation was found to be steeper in 2025 compared to the previous sampling period (interaction term: estimate: -0.0004872, std error: 0.0001171, t-value: -4.160, p < 0.001) (Fig. 4). Not only CTI values but also SwCTI values were higher in the recent data compared to the old data (Estimate: 1.324, std error: 0.2444, t-value: 5.417, p < 0.001). Compared to the CTI model, the SwCTI model explains a greater proportion of the variance (adj. R² = 0.165).

Figure 4. 

The Swiss Community Temperature Index (SwCTI) model for the temperature adaptation of butterflies in the Swiss National Park over elevation is shown. Old (1998-2004) and recent data (2025) are compared. Error ribbons show standard error.

32 species met the requirements for comparing elevational distributions. Mean difference in elevational distribution between the two sampling periods ranged from -78 m (Thymelicus lineola Ochsenheimer, 1808) to 345 m (Erebia gorge Hübner, 1804). T. lineola was the only species where a negative change meaning a downward shift was detected. Overall, a clear upward shift was found (mean difference = 134.195; t = 9.021; df = 31; p < 0.001) (Fig. 5). A summary of the elevational shift for all species is provided in the Suppl. material 3).

Figure 5. 

Mean elevation of occurrences [m a.s.l.] is compared for butterfly species in the Swiss National Park between the old (1998-2004) and recent (2025) butterfly surveys. A paired approach per species was used for assessing the difference between the survey periods.

Neither for STI values (estimate = 3.914; se = 9.441; t = 0.415; Pr(>|t|) = 0.68) nor for SwSTI values (estimate = -8.175; se = 8.838; t = -0.925; Pr(>|t|) = 0.362) a correlation with the difference in mean elevation was detected (Fig. 6). However, the magnitude of the upslope shift was found to be correlated with the diversity of habitats in the surrounding of the recordings per species. The upward displacement increased with decreasing habitat diversity (estimate = -215.98; se = 74.39; t = -2.903; Pr(>|t|) > 0.01) (Fig. 7).

Figure 6. 

Mean upward displacement in absolute values [m] over 20 years is shown in relation to temperature adaptation represented by SwSTI. Standout species are named and their data points are marked with a circle. The error ribbon shows standard error.

Figure 7. 

The correlation of Shannon Habitat Diversity and mean upward displacement [m] of butterfly species in the SNP over 20 years is illustrated. The Shannon Habitat Diversity was calculated based on the area percentages of different habitat types in 100 m buffer zones around the observation points per species. The error ribbon shows standard error.

We were able to show that habitats are more important in explaining the composition of alpine butterfly communities than elevation as an independent factor. The influence of elevation is more apparent in combination with habitats, i.e. in the availability of specific habitats at a certain elevation. Hypothesis 1 can therefore be accepted, although the influence differs greatly between the two factors. Mobile and migratory species found over a wide range are a possible reason for the low explanatory power of elevation. At the same time mobile and generalist species are most likely also reducing the significance of the habitat factor. However, most alpine species are rather sedentary (Polic et al. 2014; Klaiber et al. 2017). The reason for the major influence of habitats on alpine butterfly assemblages is the close association of most butterflies with a limited set of host plants and, especially in the larval stage, a narrow microclimatic niche (Schweizerischer Bund für Naturschutz 1987; Reinhardt et al. 2021). The high proportion of 42% unexplained residual variance must be considered: Characteristics of butterfly microhabitats including e.g. the presence of host plants, moisture, and high-resolution microclimate (Schweizerischer Bund für Naturschutz 1987; Weidemann 1995) are only superficially captured by the habitat dataset (Lotz 2006). Additionally, only adult mobility and distribution is reflected by the results, while other life stages which are either sedentary or have limited mobility are neglected. These issues represent limitations to this study.

Since the occurrence of habitats also largely depends on elevation, it follows that climate change can be expected to act not only by the direct effect on butterfly communities by increasing average temperatures. Instead, also an indirect effect caused by the resulting changes in habitats must be assumed. Nevertheless, this indirect effect acts rather slow compared to the direct effect. Plants are of particular importance in terms of habitat changes as vegetation is one of the main characteristics of habitat types (Lotz 2006). However, plants, especially those inhabiting subalpine and alpine habitats, have been shown to lag considerably behind warming in their advance to higher elevations (Roth et al. 2014; Alexander et al. 2018; Rumpf et al. 2019). For butterflies, shifting to higher elevations faster than plants (Roth et al. 2014; Kerner et al. 2023), this poses a risk: Due to the different paces in shifting to higher elevations, the availability of host plants and suitable habitats can limit the further upslope displacement of butterflies (Schweiger et al. 2008; Kerner et al. 2023). Specialist species with a narrow temperature niche are especially vulnerable in this context because they are less able to track their temperature niche due to their limited set of host plants and microhabitats. Similarly, as upper elevational limits of high-alpine plant species are shifting particularly slow (Rumpf et al. 2018), high-alpine and cold-adapted species tend to be highly affected.

As long as suitable habitats are present, they possibly can buffer against range shifts and contractions. For example, several alpine species are known to colonise areas well below their main elevational distribution, if their microhabitat structures are present, even at a small scale (Schweizerischer Bund für Naturschutz 1987; Sonderegger 2005). Their temperature niches might thus be wider than generally thought and their distribution rather driven by habitat availability. In the long run however, extinction debts of plant species, which are frequently found in the context of current range shifts (Rumpf et al. 2019), can distort this picture and lead to an overestimation of long-term habitat potential.

In line with theoretical expectations, assuming CTI and SwCTI reflecting the temperature adaptation of butterfly communities, CTI and SwCTI values increased with climate warming confirming hypothesis 2 and decreased with elevation. Warming-induced adaptations of species communities reflected by increased thermophily have been reported in other studies for various groups of organisms, including butterflies (Gottfried et al. 2012; Roth et al. 2014; Bonelli et al. 2022; Khaliq et al. 2024; Bruni et al. 2025). The results presented here thus fit well into the general picture of climate-induced community change in mountain ecosystems.

Although the results were generally consistent for both indices, marked differences with lower absolute values and steeper slopes for SwCTI were noticeable despite both indices being derived using the same methodological approach. These differences are most probably attributable to the underlying spatial resolutions of the temperature indices. The STI was calculated using MAT data from a 50 km × 50 km grid with data from across Europe (Schweiger et al. 2014). Because of the large scale, climatic heterogeneity typical for mountainous regions is not fully captured by STI, resulting in very similar MAT values for species occurring in topographically heterogeneous regions and a limited sensitivity to local elevational gradients. Conversely, SwSTI was derived from 25 m × 25 m raster data (Külling and Adde 2024), thereby better reflecting the fine-scale climatic variability of the Alps with lower overall values. However, the thermal niche is only partially covered for species not restricted to the Alps or montane regions.

As CTI and SwCTI are relative indices, differences in index values between old and recent data cannot be directly interpreted as indicating a lag, keeping pace or exceeding of local warming. It is striking that for both CTI and SwCTI, indices of the two time periods are converging at high elevations. Possible explanations for this finding are the lower density of survey plots in the old data towards the upper edge of the elevational range or differences in detectability among species. Depending on species, abundance and weather conditions, detection probability can vary, especially at high elevations, and might thus potentially bias the estimates of CTI and SwCTI. Also, the convergence may reflect non-uniform upward-shifts among species. Communities in the lower and medium elevations of the examined elevational range are transforming faster than those in the higher elevations. Here, communities still appear to be dominated by cold adapted alpine species. Assuming this explanation to be true, it would implicate that cold-adapted high-alpine specialists are receding to higher elevations at an above-average rate, whilst the rapid colonisation by lowland generalist is further increasing the CTI at the lower edge of the elevational range. Instances for such range shifts are Polygonia c-album Linnaeus, 1758 and Pararge aegeria Linnaeus, 1758, which have colonised the SNP from the lower valley areas in recent decades (Pasche et al. 2007). At the same time, high alpine specialists as Euphydryas aurinia ssp. debilis Oberthür, 1909 and Pontia callidice Hübner, [1800], which were abundant above 1’800 m a.s.l. in the first half of the 20^th century (Pictet, 1942), advanced their elevational distribution significantly until the turn of the millennium (Pasche et al. 2007) and are restricted to elevations above 2’250 m a.s.l today.

The clear upward displacement in elevational distributions of butterflies revealed by the temporal comparison confirms hypothesis 3.1. This study therefore adds to the large body of evidence demonstrating range shifts across multiple organism groups in association with climate change, in the SNP region (Baur and Baur 2013; Wipf et al. 2013; Gilgado et al. 2022) and mountain ranges in general (Parmesan and Yohe 2003; Chen et al. 2009; Menéndez et al. 2014; Marshall et al. 2020). However, the magnitude of this response to global warming is difficult to interpret because the examined elevational range was small. The analysis was restricted to an elevational range of about 750 metres constraining the detectable upwards movement. Consequently, the average upward shift of 135 metres in the last 25 years may represent an underestimate. Nevertheless, similar magnitudes of upward movement have been reported for butterflies in other parts of the Alps (Roth et al. 2014; Rödder et al. 2021; Habel et al. 2023; Kerner et al. 2023), supporting both the generality and plausibility of the identified alteration. The observed mean upward displacement of 135 m corresponds closely to the expected upward shift of 140 to 160 m, suggesting that butterflies in the SNP region are largely keeping pace with local warming. It is likely that the leading and, where applicable, rear edges have also shifted upwards along with the mean elevation. Unfortunately, the data structure prevents further insights in this regard.

Considering the finding that community composition of alpine butterflies is mostly driven by habitats and the fact that vegetation as main characteristic is lagging warming temperatures, the result of butterflies keeping pace with local temperature changes might seem contradictory. This discrepancy may be resolved by the fact that only mean elevation and thus main elevations of occurrence were compared. It should be noted, however, that the magnitude and dynamics of the upward shift may differ significantly between the edges and the mean of the elevational distribution (Rumpf et al. 2018). We therefore suspect that the mean elevations shift in line with local warming, whilst the edges of the distributions move upslope delayed corresponding to the rate of habitat change. Conclusions regarding any potential butterfly-host plant mismatches and habitat buffering, as discussed previously, cannot be drawn without data on the edges of the distribution.

The reliability of the detected upward shift of butterfly species is limited by the fact that, although both butterfly surveys were conducted in the same narrow geographical context of the SNP, the study sites were not identical. Nevertheless, despite the different study sites, a high degree of similarity was found between the butterfly communities recorded in the surveys from 1998–2004 and in 2025. This finding is particularly noteworthy given the significant differences in community composition across sites within each dataset. The Val Trupchun and Val Zeznina – Macun areas, where most of the species unique for one survey were found, had no geologically equivalent areas included in the respective other survey. Thus, reasons for differences in species composition between survey periods are most likely due to differences in bedrock (Schweizerischer Bund für Naturschutz 1987; Sonderegger 2005). Apart from local environmental conditions, discrepancies in abundances may be explained by the potentially large annual fluctuations in abundances, biannual life cycles, and weather conditions during the life cycle (Schweizerischer Bund für Naturschutz 1987; Roy et al. 2001; Reinhardt et al. 2021). The recent data were collected over a single season only and may therefore not be representative of a typical season, due to interannual variation in butterfly abundance and phenology. This fact and resulting differences in detectability may have influenced the detected upward shift. In general, we assumed that the mean elevations we found would not be greatly affected by variations in abundance. Furthermore, the distribution of data across the elevation gradient with a discontinuous pattern and the highest density of survey plots in the middle of the elevation range within the old data, and a continuous pattern in the recent data may have further affected the observed shift in mean elevation. Despite these limitations, the temporal comparison of elevational distributions of butterfly species was deemed to be valid. This was due to the high similarity of the species communities found and the fact that study areas were comparable regarding habitat features and situated at a short distance from one another.

As we did not find a statistically reliable relationship between species’ temperature adaptation and their shift in mean elevation, our results do not explicitly support previous studies that reported stronger upward shifts in cold-adapted species (Neff et al. 2022; Kerner et al. 2023). Nevertheless, such a relationship may still exist but remain obscured by high interspecific variability and limited statistical power. Besides temperature adaptation, ecological specialisation appears to be a potential factor modifying the upward displacement. Their joint influence together with a limited sample size and a few potentially influential species may prevent a clear statistical separation of these drivers.

For example, the species with strongest upward shifts within the group of species with medium to high SwSTI values were Phengaris arion (Linnaeus, 1758) and Plebejus idas (Linnaeus, 1761). Both species depend on living in association with ants during their larval development, albeit to different degrees. Therefore, they can be regarded as highly specialised, even though P. idas is not strictly restricted in terms of habitat types and host-plant use (Schweizerischer Bund für Naturschutz 1987; Reinhardt et al. 2021). In contrast, T. lineola, the only species for which a downslope shift has been observed, is a generalist species inhabiting a wide range of different habitats (Reinhardt et al. 2021). A similar pattern was observed among cold-adapted species. Boloria pales (Denis & Schiffermüller, 1775), the species with the weakest upward shift, is a polyphagous species occurring across a broad elevational range in various alpine habitats (Schweizerischer Bund für Naturschutz 1987). Erebia gorge, which exhibites the strongest detected upslope displacement, is restricted to sparsely vegetated scree areas with Festuca spp. and similar microhabitats (Schweizerischer Bund für Naturschutz 1987; Sonderegger 2005).

While these observations remain descriptive and cannot be statistically confirmed, it is reasonable to assume that temperature preference influences the response to climate change, not least on the basis of studies with similar findings (Neff et al. 2022; Kerner et al. 2023). For future research, larger datasets and trait-based approaches will be necessary to in-depth examine the relative importance of temperature preference and ecological specialisation in shaping species’ responses to climate change.

The observed pattern that species occupying a narrower range of habitats within the study area show stronger responses to climate change points in a similar direction. Such species are typically habitat specialists, or species found at the edges of their elevational range or of their ecological niche, where tolerance limits are likely to be approached. Previous studies have reported similarly higher sensitivity to climate change in habitat specialists (Engelhardt et al. 2022; Neff et al. 2022; Geppert et al. 2023), while there is also contrasting evidence. Rödder et al. (2021) found stronger upward shifts for generalist species. Although these findings may seem contradictory, they might reflect different but complementary ecological processes. Generalist species, given their low habitat specificity, can rapidly expand from lower elevations into higher elevations with warming temperatures. However, due to the restricted elevational range of this study, the potential expansion dynamics of lowland generalists is not fully captured which is why this effect is not apparent. In contrast, we found specialist species and species locally restricted to a few habitat types to respond more sensitively to climate change. These species may be less able to adapt to changes in environmental conditions by habitat expansions or moving to other microhabitats with a more suitable microclimate. Therefore, they track suitable thermal conditions through elevational shifts rather than through ecological flexibility. This interpretation is consistent with Hypothesis 1 and suggests that habitat availability and specialisation determine not only the composition of butterfly communities but also their response to climate change in alpine environments.

Although there are, of course, exceptions, the species that are closely tied to just a few habitats include a particularly large number of alpine and high-alpine species (Schweizerischer Bund für Naturschutz 1987; Sonderegger 2005). Considering the other results, a consistent ecological picture emerges: high-alpine cold-adapted species with a strong habitat specificity are particularly affected by and vulnerable to climate change. Conversely, the vulnerability decreases with a lower degree of cold adaptation and lower habitat specificity. In addition to rising temperatures, alpine habitats are undergoing particularly drastic changes as a result of global warming (Wipf et al. 2013; Biskaborn et al. 2019; Rounce et al. 2023; Körner and Hiltbrunner 2024), posing an additional threat for the species that are dependent on these habitats. At the local scale, our findings therefore particularly highlight the vulnerability of small and isolated populations and those occurring at the margins of their ecological niche or geographical range. Overall, our results underline that habitat availability represents a key limiting factor for alpine butterfly species under ongoing climate warming and suggest that species with narrow ecological niches may face disproportionate risks. This highlights the importance of considering habitat specificity and landscape structure alongside climatic variables when assessing climate change vulnerability in mountain ecosystems.

A main implication of our results for conservation is to ensure that high-quality and diverse habitats are preserved and, where possible, restored. Grazing by livestock and wild ungulates can mediate the upward shift of the tree line and contributes to the preservation of original alpine vegetation while at the same time enhancing plant diversity (Speed et al. 2010; Kaufmann et al. 2021). Although tree growth and advancing tree lines can enhance habitat diversity at a local scale, this process poses the risk of losing alpine grasslands of high value for conservation on a larger scale. Together with landscape dynamics such as landslides and avalanches with following pioneer communities, grazing ungulates and traditional livestock farming thus contribute to the preservation and diversification of alpine habitats by counteracting the slow, but evident upward shift of elevational vegetation belts. As glaciers melt and permafrost thaws, newly colonised habitats emerge (Parolo and Rossi 2008; Biskaborn et al. 2019; Rounce et al. 2023), providing potential refugia for some alpine species (Sonderegger 2005). Overall, the great variety of climate change impacts on alpine ecosystems poses additional uncertainty for future butterfly assemblages.

If rear edges shift upwards to a similar extent as mean elevations of occurrences, first local extinctions in the SNP region are expected by the start of the next century. By then, species as P. callidice or E. aurinia debilis will have shifted their elevational distribution entirely above 3000 m a.s.l., the elevation of most mountain summits in the region. However, with regional warming accelerating (Gubler and Robinson 2025) and mountaintop habitats being small and isolated, extinction events might occur even earlier (Reinhardt et al. 2005). These findings underscore the importance of climate refugia and high diversity in microclimates because of their ability to buffer against warming (Turlure et al. 2010; Suggitt et al. 2018). Even if local extinctions may not be averted (Wilkes et al. 2023; Biella et al. 2024), they offer temporary preservation of alpine butterfly assemblages on a larger scale.

To keep track of future changes in alpine butterfly communities, previous and recent butterfly surveys in the SNP region should be repeated at regular intervals, e.g. every ten years. This regularity would also allow for detecting a potential acceleration in butterfly responses to climate change, that can be expected due to the acceleration in warming. Also, recent studies on the responses of alpine butterflies to climate change have exclusively surveyed imagoes with high mobility, while the preimaginal development actually has the greatest influence on species distributions and community composition (Reinhardt et al. 2021). Thus, in future studies it might be helpful to include oviposition observations and caterpillar surveys to get more thorough insights to how alpine butterflies respond to climate change.