Butterfly distribution along altitudinal gradients: temporal changes over a short time period

Mountain ecosystems are particularly sensitive to changes in climate and land cover, but at the same time, they can offer important refuges for species on the opposite of the more altered lowlands. To explore the potential role of mountain ecosystems in butterfly conservation and to assess the vulnerability of the alpine species, we analyzed the short-term changes (2006–2008 vs. 2012–2013) of butterflies’ distribution along altitudinal gradients in the NW Italian Alps. We sampled butterfly communities once a month (62 sampling stations, 3 seasonal replicates per year, from June to August) by semi-quantitative sampling techniques. The monitored gradient ranges from the montane to the alpine belt (600–2700 m a.s.l.) within three protected areas: Gran Paradiso National Park (LTER, Sitecode: LTER_EU_IT_109), Orsiera Rocciavrè Natural Park and Veglia Devero Natural Park. We investigated butterflies’ temporal changes in accordance with a hierarchical approach to assess potential relationships between species and community level. As a first step, we characterized each species in terms of habitat requirements, elevational range and temperature preferences and we compared plot occupancy and altitudinal range changes between time periods (2006–2008 vs. 2012–2013). Secondly, we focused on community level, analyzing species richness and community composition temporal changes. The species level analysis highlighted a general increase in mean occupancy level and significant changes at both altitudinal boundaries. Looking at the ecological groups, we observed an increase of generalist and highly mobile species at the expense of the specialist and less mobile ones. For the community level, we noticed a significant increase in species richness, in the community temperature index and a tendency towards homogenization within communities. Besides the short time period considered, butterflies species distribution and communities changed considerably. In light of these results, it is fundamental to continue monitoring activities to understand if we are facing transient changes or first signals of an imminent trend.


Introduction
Global warming and land use changes are considered among the main threats to biodiversity (Sala et al. 2000, Lemoine et al. 2007). Other driving factors may also interact with climate and land use changes to impact biodiversity (Brook et al. 2008); however, substantial alterations in the population and distribution of terrestrial species have already been detected worldwide, mainly in response to both of these effects (Mantyka-Pringle et al. 2012).
Mountain ranges are very sensitive to environmental changes and global warming (Huber et al. 2005, Beniston 2006). It is well established that the Alps have experienced a temperature increase over the last century (Beniston 2003, Brunetti et al. 2009) especially at higher altitudes (Acquaotta et al. 2014). In addition to global warming, the alpine chain has suffered (and is still suffering) loss of open habitats as a consequence of forest expansion (Hunziker 1995, Gellrich et al. 2007) due to traditional land use practices being abandoned (Hinojosa et al. 2016). Such a loss of open habitats is concomitant with a noticeable loss of species (Pauli et al. 2007), which are under threat because climatic and land use changes will probably continue into the future (Nogués-Bravo et al. 2007, Chamberlain et al. 2013, Pellissier et al. 2013.
The alpine biodiversity has already responded to these factors. Upward shifts of alpine plants (Walther et al. 2005, Pauli et al. 2012, Vittoz et al. 2013, invertebrates (butterflies, Wilson et al. 2007, Wilson andGutierrez 2012;carabid beetles, Pizzolotto et al. 2014) and mountain birds (Scridel et al. 2018, Rocchia et al. 2018) have been documented. In addition, community composition has changed at high alpine sites (Keller and Körner 2003), with an accelerating increase of species richness, in many cases due to an expansion of generalist species (Walther et al. 2005, Pauli et al. 2007. Although some responses are evident, few investigations have focused on alpine biodiversity temporal changes. As mentioned above, alpine environments are under threat, therefore more information is needed to understand how the main adverse factors (global warming and loss of open habitats) have affected, or are affecting, alpine biodiversity over time. Exploring temporal patterns of biodiversity is of great significance because future warming and related environmental changes are expected to cause substantial changes, not only in spatial distribution of species, but also in species turnover over time (Korhonen et al. 2010).
Long-term monitoring programs are fundamental tools to assess and monitor temporal changes of biodiversity (Morecroft et al. 2009;). However, long-term series are rarely available (Gaston et al. 2008). Therefore, it seems useful in the meantime to focus on short time-scales and investigate short-term biodiversity responses, establishing the basis for disentangling in the future if we are facing directional changes or just short-term fluctuations.
Protected areas are a key part of conservation strategies to mitigate and monitor losses of biological diversity due to changes in climate and land use (Kharouba and Kerr 2010). In those areas, the reduced direct anthropic pressure allows for more easily focusing on and detecting the effect of climate modifications, and the natural evolution of vegetation, on biodiversity. Moreover, thanks to their staff and logistic resources, protected areas have a primary role in carrying out long-term monitoring programs, with the aim of acquiring and sharing knowledge through research networks, such as the LTER network.
In 2006, a multi-taxa monitoring project in the NW Italian Alps was started (Biodiversity Monitoring Project), with the goal of analyzing long-term biodiversity changes in space and time along altitudinal transects in three protected areas (Viterbi et al. 2013).
Following the protocol developed and tested by the LTER site Gran Paradiso National Park (GPNP), seven taxa were monitored (Lepidoptera Rhopalocera, Orthoptera, Aves, Coleoptera Carabidae, Coleoptera Staphylinidae, Araneae, Hymenoptera Formicidae), using standardized, easy-to-apply and cheap sampling techniques. In 2007, the protocol was adopted by two other protected areas (Orsiera Rocciavré Natural Park, ORNP and Veglia Devero Natural Park, VDNP). The protocol has been repeated every 5 years (2 yr monitoring -4 yr stop), and as a result, data from two sampling periods are now available (1 st season, 2006-2008; 2 nd season, 2012-2013).
In this framework, we focused on butterfly data deriving from the two sampling periods (1 st , 2006-2008; 2 nd , 2012-2013) of the Biodiversity Monitoring Project. Butterflies are an ideal indicator taxon (Dennis 1993;Hellmann 2002;Parmesan and Yohe 2003), strongly influenced by the abiotic environment and climatically sensitive (Hellmann 2002;Parmesan and Yohe 2003). Their ecological sensitivity and ability to disperse (Devictor et al. 2012) make butterflies a good model system to study species responses to climate and habitat changes (Stefanescu et al. 2011, Roth et al. 2014. The main aim of our research consisted in assessing changes of butterflies' species distribution and community composition along altitudinal gradients over a short time period. We investigated butterflies' responses across multiple levels of biological organization following a hierarchical approach, in order to assess heterogeneity or consistency of changes among different levels. As a first step, we analyzed temporal changes of distribution for single species and for species classified into ecological groups (species level). Then, we explored temporal variation of species richness and of community composition (community level). We finally discussed the results in the light of the possible drivers of change within our study areas.
All areas are characterized by mountain (alpine) conditions with vegetation ranging from mixed forest to alpine prairies and glaciers. The three parks have similar mean elevations and vegetation characteristics, but differ slightly in terms of climatic regime (the highest monthly precipitation and lowest annual mean temperature are found in VDNP).
A total of 11 altitudinal transects were set (one for each valley for each of the three parks) ranging from a minimum of 600 to a maximum of 2700 m a.s.l. The transects cover different altitudinal ranges and, overall, they involve three vegetation belts (montane, subalpine, alpine; Suppl. material 1: S1). Each altitudinal transect was composed of 4-7 sampling units (plots) separated by an altitude range of ca. 200 m, to allow independence of sampled data, for a total of 62 plots. Sampling units were circular plots with a radius of 100 m, where monitoring activities were carried out to obtain semi-quantitative data.
We sampled butterflies using linear transects along the diameter of the plot (200 m in length), walked along at uniform speed. According to the protocol by Pollard (1977), we recorded all butterflies observed within an imaginary 5×5 m square. Sampling was limited to sunny conditions, under calm-to-light wind. Each plot was visited once a month from June to August. Individuals were captured and released after specific identification, except for specimens that were difficult to identify, which were retained for further determination. Sampling efforts were characterized by 2 consecutive years of sampling and we collected data during two separate time periods: 1 st sampling period (2006( -20072007-2008 in ORNP andVDNP) and 2 nd sampling period (2012-2013 for all the three parks, GPNP, ORNP, VDNP).
In each plot, we also recorded the microclimatic conditions placing data-loggers (Thermochron iButton, DS1922L, Maxim, Sunnyvale, CA, U.S.) which recorded the air temperature every hour throughout the field season. Loggers were located in the center of each plot, at least 1 m above ground and were covered with a white shield.
By using land cover data derived from aerial photos and validated in the field (Agroselviter 2009;Meloni et al. 2009), we classified each plot according to four categories of main habitat cover type: woodland, ecotone (corresponding to transition habitats), grassland (including all grassland types jointly grouped, even if above or below the tree line) and rocks.
Moreover, during each sampling session, we noted, for each plot, land management data (dominant land use) considering three categorical variables: grazing, mowing and no activities.

Species level
The first step of our analysis focused on distributional changes of butterflies' species and groups of species. Following ecological requirements, species could change distribution along altitudinal gradient, either vertically or horizontally with the colonization or the loss of sites at similar altitudes.
To obtain a comprehensive framework on butterflies species distribution changes, we firstly assessed temporal changes in plot occupancy to detect an overall signal of expansion or retraction without considering an altitudinal directionality, then, we investigated altitudinal range changes between sampling periods.

Changes in plot occupancy
We defined occupancy as the number of plots occupied by each species in each sampling session (1 st vs 2 nd ) and compared it by using a t-test for paired samples (significance level assessed after 999 randomizations, following Legendre and Legendre 2012). The randomization process allowed to perform the t-test even with non-parametric data .
To identify which group of species changed the most over time, we analyzed if the occupancy increased/decreased equally between functional groups. We compared the changes in the number of plots per species (Δ plot, 2 nd sampling session minus 1 st sampling session), between the ecological groups of conservation interest by using non-parametric tests (Kruskal-Wallis or Mann-Whitney tests, depending on the number of ecological categories). Following the classification proposed by Balletto and Kudrna (1985) and Balletto et al. (2015), we classified butterflies' species according to: feeding specialization (from polyphagous to monophagous); altitudinal range (generalist, medium altitude, specialized); high altitude species (species typically found only in the subalpine and alpine vegetation belts); light ("shade loving", "sun loving"); dispersal capacity; habitat preferences (woodland, ecotone, open areas, screes) (Suppl. material 1: S2). We selected the ecological groups according to their potential responses in the light of possible drivers of change. In particular, feeding specialization and habitat preferences groups could highlight responses in relation to micro habitat changes while altitudinal range and high altitude species groups could have changed the occupancy status because of temperature increase. Moreover, we considered dispersal capacity because we hypothesized that species with high vagility can colonize easily new areas following both drivers of change.
To assess if individual species changed their occupancy of plots over time (if the frequency of presence/absence was statistically different between the two sampling periods), we applied the McNemar's test for paired nominal data in a 2×2 table (Sokal and Rohlf 1995). We organized data about gains and losses in plots according to the four following combinations: plots with species presence (no change) or species absence in both sampling periods (no change); plots only colonized during the first sampling season (loss) or only during the second one (gains). We tested the 133 species common to both sampling periods, and we considered the species with p-values still significant after Bonferroni correction as "winners" or "losers", depending on the direction of changes.

Altitudinal range changes
We described the altitudinal range of each species with the following parameters: -altitudinal optimum (mean and median value); -higher limit (absolute maximum, 90 th percentile); -lower limit (absolute minimum, 10 th percentile).
To quantify the amount of change, we compared these parameters between sampling periods using the t-test for paired samples. To compare non-parametric data we assessed significance level after 999 randomizations . As in the case of occupancy rates, we also compared the changes in altitudinal range among ecological groups.
To test for species responses over time along the altitudinal range, we used the indicator species analysis (IndVal method) proposed by Dufrêne and Legendre (1997). The IndVal method combines the "specificity" of a species (its uniqueness in a specific group of sites) and its "fidelity" (its frequency within the same group). For each species, IndVal can range from 0 (no indication) to 1 (maximum indication). Statistical significances of IndVal are tested by means of a Monte Carlo test, based on 999 randomizations and were performed using the function multipatt of the 'indicspecies' package (De Cáceres and Legendre 2009).
For this purpose, a species matrix (150 species) was created using the average abundance data between years for each plot and sampling period. A logarithmic transformation was applied to reduce the influence of extreme abundance values .
We grouped together sampling sites into four altitudinal bands (1: 550-1200 m; 2: 1250-1650 m; 3: 1700-2150 m; 4: 2200-2700 m) and two time periods (season I, 2006season II, 2012season II, -2013. This methodology allowed us to test the association of each species with the altitudinal bands over time, resulting in 255 possible groups, but we restricted our analyses to the 62 groups that could be interpreted ecologically in terms of altitudinal constancy or shift (Suppl. material 1: S3). Depending on the degree of species association (IndVal) to altitudinal bands between sampling periods, we identified the following altitudinal changes. We classified as altitudinal expansion the case in which a species showed, in the second sampling period, new significant associations for higher altitudinal bands, at the same time maintaining the association to the altitudinal bands of the first sampling period. An altitudinal upward shift occurred when the species lost its association to the lower altitudinal bands and gained new associations to higher ones during the second sampling period. At the same time, we identified downward retraction when a species lost its association to the lower bands, without gaining new associations with new higher ones. Then, we pointed out when the species became (spread) IndVal species or lost (disappear) their indicator value for specific altitudinal bands in the second sampling period. When no changes of associations occurred, the species remained stable over time (Suppl. material 1: S3).

Community level
Species' responses may affect communities' structure and composition (Wilson and Gutierrez 2012). Therefore, after a species level analysis, it is interesting to assess changes at the community level in order to compare the pattern of changes and to detect potential relationships. In this section, we investigated temporal changes of butterflies' communities from three different perspectives: species richness, community composition and community temperature index.

Species richness
To evaluate if sampling efficiency remained comparable throughout sampling periods, and if we consequently can carry out comparisons at the assemblage level, we performed a preliminary analysis of species richness estimators and sample completeness (expressed as the proportion of the observed species richness to the averaged estimated species richness). We calculated the averaged estimated species richness, through the mean value of four abundance-based estimators (ACE Abundance-based Coverage Estimator; Chao1; Jacknife1; Jacknife2), and correlated it with the observed species richness. Correlation was high and significant in both sampling periods (Spearman test; 1 st sampling period, ρ = 0.937, p < 0.001; 2 nd sampling period, ρ = 0.918, p < 0.001). Sample completeness was also high in both seasons (mean sample completeness among plots ± standard error; 1 st = 67.59 ± 1.56; 2 nd = 72.17 ± 1.31). Consequently, we considered the sampling efficiency between sampling periods as comparable, allowing us to carry out direct data comparisons.
To analyze how species richness per plot changed over time, we compared species richness between sampling periods using the t-test for paired samples (significance level assessed after 999 randomizations, following Legendre and Legendre 2012).
To understand if changes in species richness were mainly related to characteristics of specific plots, we analyzed it as a function of the following variables: altitude, temperature, geographic location, dominant vegetation cover (habitat type), and dominant land use. We considered the rate-of-change as a dependent variable (hereinafter ROC), defined as the differences in species richness between sampling periods (1 st , 2006-2008; 2 nd , 2012-2013), divided by the species richness of the first sampling period (1 st , 2006-2008). We analyzed the ROC through linear regression, and we compared variables in a multi-model context, according to two criteria: (i) avoiding the simultaneous use in the same model of highly correlated predictors (Spearman's rs >0.5); (ii) choosing predictors to represent biologically meaningful combinations of predictive variables, and consequently avoiding data dredging. All models were compared with the null model (intercept only) and all continuous variables were standardized, to enable comparisons between variables. The multivariate model selection was performed using Akaike information criteria, in its form corrected for small samples (AICc). As a measure of goodness-of-fit, we calculated the adjusted R 2 . These analyses were performed by the R package 'MuMIn' (Barton 2018).

Community composition
We analyzed community compositions by testing for changes both in location (significant changes in community composition per site over time) and dispersion over the years (significant changes in observed differences in community composition between sites, over time). Changes in location were tested by applying non-parametric MANOVA to Bray-Curtis distance matrixes, to test if the multivariate centroids of species composition were similar in the two groups, or not (Anderson 2001;McArdle and Anderson 2001). Non-parametric MANOVA (permutational MANOVA) is an analysis of variance using distance matrixes, which partitions sums of squares of a multivariate dataset. It was performed by the function adonis of the 'vegan' package (Oksanen et al. 2012). The significance of the test was assessed by using F-tests based on sequential sums of squares obtained from permutations of the raw data (999 permutations). Since we had to keep in mind the temporal structure and the spatial dependencies of our sampling design (62 sites at two points in time), we applied a restricted randomization, which did not allow for permutations across samples. Changes in dispersion were tested by the betadisper function of the package 'vegan', a multivariate analogous of the Levene's test for comparing group variances (Anderson 2001). To test if one time period is more variable than the other, non-Euclidean distances between objects and group centroids were handled by reducing the original distances to principal coordinates. To test for significance, we applied a similar randomization approach, as previously explained.

Community Temperature Index
The term "Species Temperature Index" (STI), refers to a quantitative description of the realized climatic niche of a species (Tayleur et al. 2016). To obtain such quantification over the Italian territory, we used presence data given by the database CkMap (Ruffo and Stoch 2005). CkMap is an atlas of the Italian fauna, promoted by the Italian Ministry of the Environment, which summarizes data related to the distribution of the Italian fauna on a gridded map (10×10 km). In our case, because we required STI that referred to Alpine populations, we only focused on North Italy, considering 1396 cells. Every cell had a value of "1" (if occupied by the species) or "0" (if species presence was uncertain). Temperature data were obtained by the maps of Metz et al. (2014), using annual mean temperature (BIO1). In this way, we calculated mean temperature values for each species (realized niche optimum).
We used STI to calculate the "Community Temperature Index" (CTI), quantified as the mean STI of all the species present in a given community. We calculated the CTI for each community (plot) and sampling period, and we analyzed the changes in CTI over time.
We compared CTI values between sampling periods by using a t-test for paired samples. As for the case of ROC, we analyzed the temporal changes in CTI (Δ CTI) by a linear regression in a multi-model context regression, as a function of the same environmental variables and models.

Temperature analysis
To quantify changes in weather conditions between sampling periods, identified as a possible driving factor of changes in biodiversity, we compared monthly temperatures (from June to September), as a function of year, altitude and additive or interactive effects. Models were compared with the null one (intercept only) with Linear Mixed Effect Models (using the plot as random effect) and all continuous variables were standardized, to allow comparisons between variables. The multivariate model selection was performed using Akaike information criterion, in its form corrected for small samples (AICc). The analysis has been carried out using the R packages 'lme4' (Bates et al. 2015) and 'MuMIn' (Barton 2018).
All the analyses were carried out with the R software, version 3.3.3 (R Core Team 2017).

Results
Considering all sampling sites and time periods, we recorded 150 butterfly species, of which 133 were common to both periods, 5 were exclusive to the first sampling session and 12 were exclusive to the second sampling session (Suppl. material 1: S2).
We identified 17 species which significantly increased their distribution within the study area, of which 7 were also significant after Bonferroni correction. No species significantly decreased its occupancy level (Table 1).

Altitudinal range changes
We observed significant changes both at the minimum and at the maximum boundary of altitudinal range of the species. At the lower altitudinal limit, we observed a significant decrease in the absolute minimum value (t-test, n = 133, t = 3.03, p = 0.004, change = -96.62 m ± 31.85 m). At the higher altitudinal limit, we observed an increase in both the absolute maximum (t-test, n = 133, t = -3.08, p = 0.006, change = 75.19 m ± 24.01 m) and in the 90 th percentile (t-test, n = 133, t = -2.63, p = 0.014, change = 55.15 m ± 20.97 m).
We identified 87 species associated with one of the selected group combinations ( Fig. 1; Suppl. material 1: S4). In most cases, IndVal species resulted as stable through- Table 1. Species which significantly changed their area of occupancy, expressed as the number of plots of presence during the first (Plot1) and the second (Plot2) sampling period. The number of plots gained or lost (Δ), the p-value of the McNemar's test (p-value) and its significance level (* p < 0.05, ** p < 0.001) after Bonferroni correction (p-adj) are also expressed.

Species
Plot 1 Fig. 1, Suppl. material 1: S4), only one species showed an altitudinal upward shift, gaining the association to a higher altitudinal band and losing it into the lower one. The general increase in species occupancy also resulted in 19 species becoming indicator species during the second sampling period of specific altitudinal bands (spread). Only three species were indicators of the lowest altitudinal range only during the first time period, showing a decline during the second sampling season (disappear). For no species did we record a downward retraction.

Species richness
Species richness significantly increased from the first to the second sampling period (t-test, n = 62, t = -9.76, p = 0.001, change = 8.82 ± 0.90). The analysis of the ROC showed a significant effect of both land cover and land use: wooded habitats and managed plots increased the most (Table 2).

Community composition
Species composition significantly differed between time periods, but the obtained R 2 is extremely low, meaning a scarcely relevant change (non-parametric MANOVA, Fvalue = 5.87, r 2 = 1.91, p = 0.001). More interestingly, we observed a significant change in dispersion between sampling periods, with a lower dispersion around the median during the second sampling period (distance to centroids 1 st sampling period = 0.584, 2 nd sampling period = 0.528, p = 0.001; indicating a tendency towards homogenization; Fig. 2).

Community Temperature Index
We observed that the CTI significantly increased from the first to the second sampling period (t = -3.59; p = 0.001), indicating a common trend toward thermophily in all . Different colours represent different statuses. Species that lost their indicator value with specific altitudinal bands during the second period (disappear; black, with continuous border); species stable over time (stable; white, with continuous border); species that showed, in the second sampling period, new significant associations for higher altitudinal bands, at the same time maintaining the association to the altitudinal bands of the first sampling period (altitudinal expansion; grey, with continuous border); species that lost their association to the lower altitudinal bands and gained new associations to higher ones during the second sampling period (altitudinal upward shift; grey with dotted border); species that became IndVal of specific altitudinal bands during the second sampling period (spread; white, with dotted border).

Figure 2.
Homogenisation of community composition. Box-plot of the distances to the centroid of community composition, during each sampling period. A reduction in species heterogeneity at the community level from the first to the second period can be clearly seen. Table 2. Best linear regression model for Rate Of Change (ROC). Coefficients (± standard error) of the selected variables are indicated in the cells. Adjusted r squared is indicated as a measure of goodness of fit. alt = altitude; alt 2 = altitude (second order); park = geographic location; rme = change in mean temperature; rmi = change in minimum temperature; vegetation = dominant cover type (land cover); use = land use; Tmin = mean seasonal minimum temperature during the first season. adj r 2 = adjusted R 2 . ** p < 0.01; * p < 0.05; ° p = 0.06. the sampled areas. Moreover, interestingly, we observed that the change in CTI over time was mainly dependent on the geographic position of the sampling plots ( Fig. 3; R 2 = 14.17, p = 0.007; geographic location, p = 0.007). In particular, we recorded the significantly higher increase in CTI in the colder areas, i.e. in the plots located in VDNP, which is characterized by the lowest temperature among our study areas.

Temperature
The analysis of field-recorded temperatures indicated significant differences between sampling periods, concerning the monthly mean and minimum temperatures observed during July, August and September (Table 3). In particular, we observed a second sampling period significantly warmer than the first one with a significant increase of minimum (+1.22 °C) and mean (+0.83 °C) temperatures. The values of increased temperatures reported refer to the seasonal average (from June to September).

Discussion
Long-term monitoring, which is periodically and systematically repeated over time (usually decades), is crucial for correctly understanding the variables and mechanisms that determine species' distribution responses and patterns of community composition as a consequence of climate and land use changes , Magurran and Henderson 2010, Legendre and Gauthier 2014. Interestingly, some recent studies in the Alps have shown how significant changes can also occur over short periods, highlighting the importance of verifying distributional and community changes at temporal intervals of less than 10 years (short-term changes). Roth et al. (2014) described changes in community composition of birds, butterflies and plants in the Swiss Alps along an altitudinal gradient, over just 8 years. Erschbamer et al. (2009) showed how plant species richness in the Dolomites was significantly higher after only 5 years. In this framework, our study on butterflies composition and distribution along altitudinal gradients in the NW Italian Alps represents an important step towards a better comprehension of biodiversity patterns in mountain ecosystems, even if it was restricted to a short time-frame (2006-2007 vs 2012-2013). Indeed, understanding the spatial and temporal dynamics of species-rich communities is critical for understanding how environmental changes can affect biodiversity (McCann 2007).

Driving factors of change
Although our research is based on a short time frame, it is important to consider the role of potential driving factors of change, by analyzing climatic and land cover changes (whenever occurring) as explicative predictors.
Concerning variations in temperature between sampling periods (2006-2008 vs 2012-2013), we observed that the second sampling period was warmer than the first one. These observed differences are coherent with the trend recorded by Beniston (2006) for the Alps and are likely correlated to climate change. Indeed, Beniston (2006) observed an increase in minimum and mean temperatures over time, but no change in the maximum temperature.
Unfortunately, no punctual and detailed data on precipitations were available to introduce them into the analytical framework. However, analyzing data from the official weather stations of the Regional Meteorological Service located in the study areas, we observed an interesting pattern in snow cover, characterized by a reduction in the number of days with snow and by a seasonal shift in snow cover towards the spring, in particular since 2005. Such conditions are indicative of a warmer climate and they correspond to a common pattern for the entire Alpine chain (Scherrer et al. 2016).
Due to the short time frame of analysis, we considered the main habitat type (dominant land cover) as a constant variable because no huge structural changes have potentially occurred. However, because of the natural evolution of vegetation following climatic changes and land abandonment, minimal but effective land cover changes in terms of microhabitat scale (e.g. shrubs encroachment) can be taken into account (Zurlo 2018).
Some sampling sites (plots) were characterized by low-intensity and sustainable grazing activities without any changes between sampling periods.
No other direct land use changes and habitat alterations occurred within our study sites because of the park status and LTER accredited site of our study area.

Plot occupancy
Despite the short time-frame under analysis, butterflies showed a general increase in mean occupancy levels but species-specific plot occupancy changes differed markedly among species. However, when grouping species according to homogeneous ecological traits, we observed consistent distributional responses, in accordance with several studies that have proved the effect of ecological and life-history traits in shaping species distribution (Forister et al. 2010, Chen et al. 2012, Auer and King 2014. Few studies to date have examined the underlying reasons for distribution changes and so far results have been equivocal (Angert et al. 2011;Buckley and Kingsolver 2012).
Monophagous, altitudinal specialists and high-altitude species appeared to be more limited than other species. These categories comprise species with a high level of specialization, consequently less prone to colonize new environments, even if climatic or environmental constraints are relaxed. In particular, monophagous species are strictly limited by the presence and the quality of their single larval host plant, and are already observed and predicted to be highly vulnerable to climatic/environmental changes (Romo et al. 2014). Our results concerning butterfly specialization are interesting and reflect what has been observed in central Europe concerning habitat specialization, where a decrease of specialized and low vagile species has been observed, along with an increase of generalist species and good dispersers (Habel et al. 2016).
A significantly higher increase in plot occupancy by "shade-loving" species compared to other species can be associated with a tendency towards a higher coverage of shrubs in open areas at low and medium altitude. In the European Alps, the effect of climate change is regionally affected by human activities. Cattle grazing in the alpine pastures has been decreasing throughout the past century, allowing rapid recolonization by trees and shrubs, where the treeline had been artificially lowered (e.g., Vittoz et al. 2008a;EEA 2010;Rocchia et al. 2018).
Interestingly, no species significantly decreased its area of occupancy, but 17 species increased their occupancy area, 7 of which in a particularly marked manner, still significant after Bonferroni correction (Table 1).
Aglais urticae and Argynnis niobe showed pronounced increase in plot occupancy. Our results are in contrast with other researches that reported a decline of these two species in other parts of Europe (Kulfan et al. 1997, Saarinen et al. 2003. Those spe-cies are altitudinal generalists with a high dispersal capability like most of the other species which increased their area of occupancy within our study area (e.g. Pieris rapae and Colias alfacariensis). We therefore consider that our observations could be only transient, due to the high dispersal capacity of those species and their tendency to carry out vertical migrations. Particularly interesting was the increase in mean occupancy of C. rubi. Since its low feeding specialization (Balletto et al. 2015), this species could have increased its abundance or spread into new areas, without having limitations due to food availability. However, although we highlighted important plot occupancy changes both at ecological groups and species level, we should also consider that they could be influenced by many factors, such as species detectability, population dynamics and phenology, which could, in turn, affect plot occupancy detection. To disentangle such patterns, long-term monitoring and more detailed study on the population dynamics of target species will be essential.

Altitudinal range changes
As for plot occupancy analysis, we noticed a general process of altitudinal change, considering all butterflies species together.
Looking at the ecological groups, we observed that, even along altitude, generalist and highly mobile species were more prone than specialist and less mobile ones to change their distributional ranges. Particularly, high altitude species showed a significantly higher increase in their mean, median and 10 th percentile altitudinal parameters, corresponding to a reduction in their lower altitudinal boundaries and in the surface available. Those species are already limited in their distribution. Their presence is, in many cases, limited by minimum temperatures (Pellissier et al. 2013) and, consequently, they cannot lower their altitudinal range nor, in many cases, raise it, due to the drastic decrease of vegetation cover at higher altitudes (a high occurrence of rock cover and a strong reduction of the availability of herbs and grasses). Moreover, the reduction in permafrost, which is a relatively new and rapid phenomenon, makes high altitude rock and screes unstable, preventing colonization by plant species (Cannone et al. 2007). However, we should also consider a limit of our sampling design, as we did not monitor plots over 2700 m a.s.l., consequently reducing our possibility to observe an expansion towards higher altitudes and a colonization of new plots by high altitude species.
Concerning altitudinal patterns, we used the IndVal test to identify species associated with specific altitudinal zones. Our analysis highlighted 52 species associated with one or more altitudinal zones at both time periods: these species are widely spread within their altitudinal zone, and stable over the investigated period, consequently representing species which -in the future -could be monitored to investigate important altitudinal changes. For 12 species, we observed new significant associations for higher altitudinal bands (altitudinal expansion). We recorded a change in association from a lower to a higher altitudinal band (upward shift) only for one species, Celastrina argiolus, which was associated with the first altitudinal zone (550-1000 m) during the first sampling period, while during the second sampling period, it resulted characteristic for the second zone (1250-1650 m) only (see statistical analyses -altitudinal range changes for the description of the relationships between IndVal association and altitudinal changes). It is a species that favors woodland and shade, and probably an increase in minimum temperature with an increase in the shrub land coverage had an important effect in the shift of this species.
In contrast, other studies in European countries observed on the opposite a higher amount of species shifting their altitudinal distribution, disappearing from low elevation sites and colonizing higher altitudes (e.g., Hill et al. 2002;Konvička et al. 2003). The differences from our work, where most of the species changing their altitudinal association did not disappear from the lowest altitude sites, could be related both to the shorter temporal frame we analyzed and to the lower altitudinal gradient we monitored. Indeed, our low altitude sites are located in the montane belt, where the effects of climate and environmental changes could be less severe than in the lowland plains.
Finally for 19 species, the changes in their distribution determined new associations for specific altitudinal bands during the second sampling period (spread). These species represent important taxa to be followed over time, to understand if this pattern represents the first signal of an upward expansion.

Species richness
We observed significant changes in species richness per plot within the analyzed period. Butterfly communities are known to quickly change their arrangement due to environmental changes (Thomas 2005), and previous studies have indicated that butterflies could respond even faster than other taxa (Devictor et al. 2012).
We observed a clear and significant increase in species richness within our temporal frame, mainly related to land cover and land management. The highest rate of change was clearly evident in the wooded areas (woodland clearings), while ecotones (transitional areas dominated by shrubs, and mainly located inside the subalpine belts) showed the lowest rate of change. Other authors also observed similar results, although they were mainly related to individual species abundance. For example, Sgardeli et al. (2016) observed a higher increase in species abundance simultaneously with an increase of temperature within forest areas with respect to other areas. Indeed, on days with high temperatures and solar radiation (which are exacerbated in open areas), wooded areas can exert a buffering effect, protecting individuals from extreme temperatures and reducing temperature leaps (Oliver and Morecroft 2014). Even thermophilous and "sun-loving" butterfly species can be inhibited in their flight activity by high temperatures and solar radiation, with potentially dangerous consequences that have not yet been fully understood (Cerrato et al. 2016).
In grazed areas, we observed an increase in species richness that was twice as much as in the unmanaged area. This can be probably explained by the low-intensity, sustainable grazing of these sites. Although intense grazing activities can be detrimental for butterflies (Scalercio et al. 2014), it has been previously observed that a balanced grazing pressure can increase the presence of plant species belonging to Poaceae and Fabaceae (Fischer and Wipf 2002), which represent the most exploited plant families that host the larvae of many butterfly species. Moreover, grazing maintains woodland clearings and open herbaceous areas below the tree line, that would be rapidly colonized by shrubs and trees, without such management activities (Nagy and Grabherr 2009). Indeed, balanced grazing activities have the fundamental role of maintaining micro-habitat heterogeneity, which is, in turn, of high importance in reducing the potential detrimental effects of climate change to butterflies in mountain ecosystems (Klečková et al. 2014;Konvička et al. 2016).

Community composition
Although we observed major differences in species distribution between our sampling periods, we did not observe significant differences in community composition. These results were probably due to a consequence of the short time-frame of analysis.
Studies that have demonstrated variation in butterfly community composition clearly take into account longer time-frames (e.g., Habel et al. 2016) and, to our knowledge, no investigations have reported significant changes in community composition in such a short time-frame. Interestingly, however, we noticed a tendency towards biotic homogenization in butterfly community composition. Biotic homogenization refers to the increase in biological similarity between communities, a replacement process leading to a decrease in distinctiveness in community composition over time, as a result of the replacement of some specialist species with other generalists, which become more uniformly distributed across previously different assemblages (Olden and Rooney 2006). Indeed, species respond individually to changing environmental conditions, depending mainly on their physiological characteristics and habitat requirements (e.g., Wilson and Gutierrez 2012). This determines new species assemblages, which can only be appreciated by examining the entire communities throughout time (e.g., Wilson et al. 2007).
For example, a similar change in community composition over time, accompanied by an increase in community similarity, has been observed in the data analysis from the UK Butterfly Monitoring Scheme over a period of 20 years (Gonzàlez-Megìas et al. 2008). This tendency towards biotic homogenization has been observed in recent decades in different taxa, following land cover and climatic changes (e.g., Bühler and Roth 2011;Eskildsen et al. 2015). These phenomena often determine an increase of generalist and of highly vagile species, to the detriment of other taxa (e.g. Menéndez et al. 2006;Bonelli et al. 2011).

Community Temperature Index
The observed increase in CTI tended to derive from the increase in mean occupancy levels recorded by a large number of generalist, more termophilous species in the second sampling period compared to the first one. This phenomenon was particularly marked in the coldest areas, representing a potential threat for high altitude and microthermic species. Indeed, the generalist species successfully colonizing a higher altitude from a lower one can outcompete the more specialized and vulnerable species of the alpine belt, challenging the survival of their local populations in the long term. Our results are consistent with an increase in CTI already observed in other geographic areas (e.g., in Greece by Zografou et al. 2014 over 13 years; in the Swiss Alps by Roth et al. 2014 over 8 years; at European scale by Devictor et al. 2012 over two decades). In any case, our trend was observed over a shorter time-frame, and, if confirmed during the next monitoring sessions, could represent a warning signal for alpine butterfly fauna.

Conclusions
In recent decades, changes in species distribution and composition of communities have already been analyzed in some parts of the Alps. Nevertheless, data usually concern plants (e.g. Grabherr et al. 1994;Pauli et al. 2001;Vittoz et al. 2008aVittoz et al. , 2008b, or birds (Popy et al. 2010;Scridel et al. 2017), with very little information regarding other taxonomic groups, and mainly concerning changes in the distribution of single species (e.g., Battisti et al. 2005).
Responses to climate and habitat changes vary widely between species of the same communities (Walther 2010). Consequently, it is of crucial importance to understand the heterogeneity of species responses and the implications of the changes in community composition that follow (Wilson and Gutierrez 2012). Moreover, many previous studies that aimed at understanding changes over time in community composition have relied on the comparisons of contemporary data with historical data (atlases, collection specimens), which were often collected in a non-standardized way and/or referred to a much coarser spatial grain (Wilson and Gutierrez 2012).
Sampling units placed in well-specified areas represent a more appropriate tool for investigative purposes (e.g., Archaux 2004;Viterbi et al. 2013;Brandmayr and Pizzolotto 2016), and our data from the Italian Alps, taken from exactly the same sites at 4-year intervals, represent a first attempt to fill these gaps, focusing on changes in mountain ecosystems and mainly considering community level responses. Moreover, thanks to the LTER network, such standardized and temporally replicated data are available for broader spatial scale analysis and comparison with data coming from other sites of the same network, greatly improving our capacity to understand biodiversity patterns.
Interestingly, and somewhat alarmingly, our results suggest that, even if the timeframe under analysis is relatively short, we already observed considerable changes, in particular considering that the research was carried out in protected areas and LTER sites where habitat alteration by direct human effects is strongly reduced.
To summarize, in just 5 years, we observed: -a general increase in mean occupancy level and in species richness of butterflies; -no significant changes in the mean altitudinal optimum, but significant changes at both altitudinal limits.
Moreover, the observed changes differed across species, determining: -an increase in shared species (a tendency to homogenization) within communities, even if the overall community composition did not change; -a significant increase in Community Temperature Index (CTI).
Considering these results, it is now even more important to continue our monitoring in order to understand in the near future if our observed patterns represent only transient changes or are the first signals of an imminent trend.