Effects of Temperature Rise on Multi-Taxa Distributions in Mountain Ecosystems

Abstract

Mountain biodiversity is associated with rare and fragile biota that are highly sensitive to climate change. To estimate the vulnerability of biodiversity to temperature rise, long-term field data are crucial. Species distribution models are an essential tool, in particular for invertebrates, for which detailed information on spatial and temporal distributions is largely missing. We applied presence-only distribution models to field data obtained from a systematic survey of 5 taxa (birds, butterflies, carabids, spiders, staphylinids), monitored in the northwestern Italian Alps. We estimated the effects of a moderate temperature increase on the multi-taxa distributions. Only small changes in the overall biodiversity patterns emerged, but we observed significant differences between groups of species and along the altitudinal gradient. The effects of temperature increase could be more pronounced for spiders and butterflies, and particularly detrimental for high-altitude species. We observed significant changes in community composition and species richness, especially in the alpine belt, but a clear separation between vegetation levels was retained also in the warming scenarios. Our conservative approach suggests that even a moderate temperature increase (about 1 °C) could influence animal biodiversity in mountain ecosystems: only long-term field data can provide the information to improve quantitative predictions, allowing us to readily identify the most informative signals of forthcoming changes.

1. Introduction

Climate change is driving community composition rearrangements and biodiversity losses worldwide, currently representing a crucial theme in conservation biology (e.g., [1,2,3]). Understanding and predicting global warming impacts on biodiversity is essential to develop effective conservation strategies (e.g., [4,5]).

Not all ecosystems display the same degree of exposure, sensitivity, and vulnerability [3,6]. Mountain ecosystems host a high level of biodiversity and are especially threatened (e.g., [7]), owing to the presence of species with small distribution areas, low dispersal ability, and high levels of physiological and ecological specialization [8,9,10]. In the European Alps, a large number of endemic species is present, and during the last decades, many habitat types and animal and plant populations have undergone exceptional decline [11,12,13]. Many mountain ecosystems reported warming higher than the global average [14,15,16,17]: in the Italian Alps, the measured air temperature rise was of the order of 0.2–0.4 °C per decade since 1950 [18,19]. For the Alps, climate projections suggest future warming of 0.5–0.7 °C per decade [20], again higher than the global average [21,22].

Alpine biodiversity has already responded to the temperature rise. For alpine flora, common effects are upward displacements of high altitude plants, treeline advance, the decline of arctic-alpine species, changes in community composition, and species richness (e.g., [11,13,23,24,25]). For alpine fauna, similar responses have been measured, such as upward shifts of single species and changes at assemblage level [26,27,28,29], but the number of studies is markedly lower.

Individual species are expected to respond differently to climate warming, and the most important effects will presumably be at the community level [1,30]. Even though the temperature rise has been reported to impact community composition, the number of long-term standardized datasets with a sufficiently high temporal and spatial resolution is still limited, in particular for multi-taxa analyses [31,32,33].

Understanding how multiple species and communities will respond to climate change is essential for management, from global to local scale [34,35]. Indeed, the effectiveness of protected areas and of the Natura 2000 network under climate change has been questioned (e.g., [3,36,37]). It has been argued that nature reserves should improve their role from pure conservation to the stewardship of the changes, representing key sites for the study of climate change responses [38,39]. Conservation actions should consequently include the management of “smoother transitions” to new conditions, which necessarily requires the understanding of current dynamics and future states [39].

Species distribution models (SDM) are simple tools to estimate taxa response to climate warming, assessing the potential vulnerability of individual species and communities as a whole [40,41,42,43]. SDM could have large uncertainties in their outcomes, as they use a correlative approach based on the observed, realized species niche. SDM usually assume niche conservatism and the absence of evolutionary processes, ignoring biological interactions and physiological mechanisms [44,45]. More sophisticated, mechanistic models are currently available, but they have high data demands and thus cannot be used for any of the studied taxa, owing to the lack of data (e.g., [46,47]). Consequently, simple SDM are often considered as the main tool to estimate species and community vulnerability to climate change [40,48] and fill knowledge gaps in understudied taxa and habitats (e.g., invertebrates, [31,49]; altitudinal gradients and mountain ecosystems, [50,51]). SDM allows us to explore different scenarios and quantify biodiversity changes for a wide range of species at the same time, obtaining useful indications for protected areas’ management [5,39,40].

To obtain quantitative information on biodiversity along altitudinal gradients on spatial and temporal scales suitable for site management, a multi-taxa monitoring project was started in 2007 in the northwestern Italian Alps [52]. Five taxa (Coleoptera Carabidae, Coleoptera Staphylinidae, Araneae, Lepidoptera Papilionoidea and Hesperiioidea, Aves) representing a wide array of ecological and evolutionary traits, have been monitored in three separated mountain protected areas. Climatic and vegetation data were collected on the same spatial scales. The results of the analyses have shown that invertebrate species’ richness and community composition are related to temperature and altitude, highlighting the potential vulnerability of the alpine communities to climate change [52,53].

The present work uses the above dataset to assess the risk of short-term changes in Alpine animal biodiversity under climate warming. The coarse spatial resolution of current (global and regional) climate projections, however, does not allow to assess the response of mountain ecosystems, that are characterized by high topographic heterogeneity and small-scale microhabitat and microclimate variability [54,55]. Short-term, fine-grained climate predictions would certainly be needed by site managers [56], but there is still a generalized paucity of climate information for mountain areas [57].

On the other hand, it is important to provide informative projections, even within simplified modeling exercises, to investigate the relationships between land-use and climate change [58,59]. Our study, based on local-scale biodiversity and environmental data, is a contribution to address the challenge of estimating the response of multi-taxa distributions. We applied “what if” scenarios, increasing the maximum, minimum, and mean temperatures with respect to the increase observed in the Alps during the last decades. To assess the potential role of land cover, we compared projections with and without vegetation constraints.

Many works indicate differences in the ability of different taxa to adapt to climate change, but most of the studies are based on meta-analysis or data coming from independent studies [33,56,60]. Our approach compared data coming from different taxa collected with the same methodologies and in the same monitoring framework. Our main objective was thus to estimate patterns of congruence in multi-taxa response to climate warming. In particular, we asked whether: (i) the number of species changing their distribution differs among taxa, degree of vulnerability and scenarios and (ii) species richness and community composition significantly change along the altitudinal gradient. This allowed us to obtain indications on both indicator selection and habitat protection priorities.

2. Materials and Methods

2.1. Biodiversity Inventory in the Northwestern Italian Alps: Data Sources

In 2007, 3 protected areas in the NW Italian Alps (Gran Paradiso National Park, Orsiera Rocciavré Natural Park, Veglia Devero Natural Park) started monitoring animal biodiversity along altitudinal gradients. Twelve altitudinal transects have been identified (overall covered elevational range 550–2700 m), each characterized by 5–7 sampling units spaced apart by an altitudinal range of 200 m, for a total of 69 plots (circular areas with 100 m radius). Plot characteristics and their spatial locations are detailed in the Supplementary Material (Table S1).

Five taxa have been monitored in 2007, using semi-quantitative methods. Carabids (Coleoptera Carabidae), staphylinids (Coleoptera Staphylinidae), and spiders (Arachnida Araneae) were monitored through pitfall traps (5 per plot, filled by 10 cc of white vinegar, checked every 15 days from May to September). Butterflies (Lepidoptera Papilionoidea and Hesperiioidea) were monitored through linear transects (200 m along one of the diameters of each plot, once per month from May to September). Birds (Aves) were monitored through point counts (lasting 20 min, twice per plot between April and July, choosing the most appropriate period depending on altitude). We refer to [52] for a detailed description of the sampling methodology. Butterflies were mainly identified in the field, and only a few specimens were collected for subsequent identification in the laboratory. Epigeic invertebrates (Carabidae, Staphylinidae, Araneae) were identified in the laboratory by expert taxonomists and are currently stored at the park’s headquarters.

These taxa were selected because they are good candidates for estimating biodiversity patterns along altitudinal gradients [61,62,63,64]. In particular, they are: (i) well represented in mountain ecosystems, both in terms of species richness and number of individuals, and (ii) characterized by species with different ecological needs and different levels of specialization. Altogether, the selected taxa include different trophic levels and different taxonomic relatedness. Moreover, they can all be sampled using easy-to-apply, cheap, standardized, and well-established techniques, that allow monitoring repeatability through space and time.

To characterize microclimatic conditions at the plot scale, we installed one data-logger (iButton, DS1922, Maxim, Sunnyvale, CA, USA) at the center of each plot, at about 1 m above ground, recording hourly air temperature from June to September. We calculated mean, maximum, minimum, and standard deviation of daily temperature at each plot and averaged them to obtain seasonal values. Each plot was described by: (i) mean altitude; (ii) geographic position (a categorical variable indicating which protected area they belong to); (iii) vegetation belt (Montane, Subalpine, Alpine, indicating altitudinal sections characterized by given vegetation and climate, [65,66]); (iv) structural diversity. To obtain structural diversity, we empirically estimated during field surveys, using 5% classes and checking vegetation maps, the percentage of ground covered by different structural layers (herbaceous layer, low shrubs < 1 m, tall shrubs between 1–5 m, trees, stone, and bare ground cover). We then calculated structural diversity as the Shannon index of the herbaceous layer, shrub, and tree cover. This index quantifies microhabitat heterogeneity, and it is often adopted in conservation studies (e.g., [52,67,68]).

2.2. Model Simulation: Current Conditions and Temperature Change Scenarios

We applied Species Distribution Models to 304 species (45 carabids, 40 staphylinids, 99 spiders, 80 butterflies, 40 birds), using distribution data coming from 62 plots. Indeed, 7 plots have been discarded owing to missing temperature records and 359 species were present in less than 4 plots (selected as a threshold for presence accuracy). Each of the 304 species, as shown in Table S2a, was characterized by: (i) the taxon they belong to; (ii) the degree of vulnerability (a dichotomous variable indicating if the species is restricted to high altitude); (iii) the level of endemism (a dichotomous variable indicating if the species is endemic of the Alpine biogeographical region).

We modeled each species individually and subsequently estimated species richness and community metrics on model outputs that were converted to presence/absence data, following the approach “predict-first, assemble later”, or stacked species distribution models [34]. This approach implies unsaturation at the community level, which is currently considered a very likely assumption [69], and allowed us to create unconstrained species richness and community composition maps.

To create species distribution models, we used Maxent, a machine-learning approach based on maximum entropy (Maxent software 3.3), developed by S. Phillips and colleagues ([70], freely available at http://www.cs.princeton.edu/~schapire/maxent).

We selected Maxent because it is a high-performing and robust bioclimatic approach, widely used in recent years [71,72,73]. In particular, Maxent (i) can be run with both continuous and categorical variables [71]; (ii) it is stable also with correlated predictors [71], and consequently we were allowed to simultaneously model the effects of temperature and altitude; and (iii) it is highly trustworthy with a reduced number of presences [74,75,76], consequently allowing modeling also of rare and endemic species [77,78,79].

The species distribution is modeled by comparing environmental conditions at plots where the species is present with the conditions encountered across the study area. These latter are defined by a set of background points that identify the available environment. For each plot, Maxent estimates the probability distribution that best fits the environmental conditions at the plot, remaining as close as possible to a uniform distribution (maximum entropy principle). For each plot, a logistic output is then generated. This can be interpreted as an estimate of the presence probability, given the environment, with values ranging from 0 (lowest probability) to 1 (highest probability) (see [70,71,80] for a detailed explanation of methodology). We used the default parameterization of Maxent regarding feature types, regularization, and prevalence [74,81].

We selected for both the sample plots and the background points the same 62 monitored plots, to limit our predictions to a set of sampling units for which all data were collected at the same scale. Several works indicate in fact that macro- (coarse-scale air conditions) and micro-climate (experienced by organisms) can be significantly different, in particular in mountain ecosystems. Consequently, it is microclimate that should be used in SDM [82,83]. The role of microhabitat associations in determining current and future distributions of living organisms has also been emphasized [84]. As environmental predictors, we thus chose the variables measured in each plot that define local microclimate, altitude, geographical location, and vegetation cover (microhabitat).

We combined the environmental predictors in different ways, obtaining three model classes with an increasing number of environmental constraints (Figure 1):

• Temperature (T), which considers only temperature-derived variables (seasonal mean, maximum, minimum temperature, and standard deviation, in °C) and altitude to model species distribution;
• Temperature+Region (TR), which considers temperature-derived variables, altitude, and geographical location;
• Temperature+Region+Vegetation Structure (TRV), which considers temperature-derived variables, altitude, geographical location, and vegetation cover.

To project species distributions for each model class (T, TR, TRV), we applied three different “what if” scenarios of temperature change:

• 1Degree (d), in which minimum, mean, and maximum temperature are all equally increased by 1 °C;
• 1.5Min (min), in which minimum temperature is increased by 1.5 °C, mean temperature by 1 °C, maximum temperature by 0.5 °C;
• 1.5Max (max), in which minimum temperature is increased by 0.5 °C, mean temperature by 1 °C, maximum temperature by 1.5 °C.

We based our choices on results from the analysis of temperature trends in the Alps, especially on [14] and [85], which reported a larger increase for minimum and maximum temperature respectively, and on [20] for future scenarios. Elevation-dependent warming is of relevance for most mountain regions [17], but the temperature increase can be taken as relatively homogeneous in the range of elevations considered here, even in the case of the Alpine belt [86]. “What-if” scenarios based on assuming a given temperature (and precipitation) change are a relatively recent way to look at climate change impacts (e.g., [87,88]). This approach is adopted here, in a simplified form, owing to the difficulty of regional climate models to properly represent future climatic changes in mountain areas of a limited extent [89].

We determined model accuracy by evaluating the area under the receiver operating characteristics curve (AUC), using the same dataset employed for training. This provided an optimistic measure of prediction success [90], which was anyway useful to assess warming effects on small datasets. Few species had an AUC < 0.6 (Table S2a–c; Figure S1) and we considered them to be stable under warming scenarios.

2.3. Analysis of Model Outputs

We transformed the model output into binary maps (presence/absence), using as threshold the minimum predicted value for the training sites, also termed “lowest presence threshold” [70,91].

This threshold identified as “species presence sites” all sites that are at least as suitable as those where the species is known to occur. After running the models for each individual species, thus allowing for different responses of different species, we calculated the α- and ß- diversity parameters to estimate changes from the current situation (defined as Maxent modelization before temperature increase scenarios), following the recommendations of [92].

2.3.1. Species Distribution

For each model class, we compared each warming scenario with the current distribution in terms of the number of occupied plots for each species.

We identified the species that significantly changed their distribution under warming (varying species), using a two-tailed binomial test. The number of successes was estimated from the number of occupied plots following the temperature increase, which was then compared to the number of plots currently occupied. To verify whether the three scenarios and the three model classes significantly differed in the number of varying species, we calculated the χ² statistics on the contingency tables.

2.3.2. Species Richness

For each model class, we compared the distribution obtained in each warming scenario with the current distribution, in terms of species richness per plot. The differences were tested using the t-test for paired samples (significance level assessed after 999 randomizations, following [93]).

We quantified the changes in species richness per plot, for each model class and warming scenario, as:

ES = (Sc − Sp)/(0.5 * (Sc + Sp))

(1)

where ES is the effect size, Sc is the current species richness, and Sp is the projected species richness. The increases and decreases of ES are symmetrical and range from −2 to +2 [94].

To identify taxa that are more sensitive to the temperature increase, we analyzed the effect size as a function of the taxonomic group (considering as baseline the ensemble of all taxa). To this end, we used linear mixed-effects models with plot identity as a random effect.

We carried out an exploratory analysis using model classes, scenarios, and taxa and their interactions as explanatory variables. We selected the best model based on the Akaike’s Information Criterion corrected for small samples (AICc), computed with the “MuMIn” package [95]. Since the best model included no effects of warming scenarios, identifying instead a significant interaction of taxa and model class (Table S3), we then applied the linear mixed-effects models to each taxon separately, with the model class as an explanatory variable.

To understand how effect size changed as a function of elevation, we first graphically analyzed its variation along the altitudinal gradient and then focused on the differences between vegetation belts using again the linear-mixed effects models (with plot identity as a random effect).

Linear mixed-effects models have been implemented with the “lme4” package [96].

2.3.3. Community Composition

We represented changes in assemblage composition under warming scenarios by applying correspondence analysis (CA) to each taxon and model class. We focused on the first two CA-axes, whose explained variances were always significant (based on 999 randomizations obtained by changing species presence across plots while keeping prevalence constant). The cumulative explained variances ranged from 32.26 to 63.80 (Table S4a). The first axis was determined by altitude, minimum, and average daily temperature, while for the second axis we did not identify a clear pattern (Table S4b). To understand whether warming changed community composition as a whole and across vegetation belts, we applied the Wilcoxon Rank Sum Test on plot scores of current and projected distributions and compared the shift in plot scores between vegetation belts, using the Kruskal Wallis Test. Community homogenization under temperature increase was tested, comparing the distance from the distribution centroid for current and projected assemblages. Significance was assessed by the t-test for paired samples (using 999 randomizations [93]).

The variation in community composition of single plots under climate warming scenarios was quantified by the Jaccard Index and compared across model classes, scenarios, and taxa by the Friedman Test [97]. We also partitioned the overall dissimilarity in its turnover (β_sim) and nestedness (β_sne) components, calculated using the betapart package [98]. This allowed us to quantify whether communities exposed to temperature increase were characterized by the substitution of some species by others (turnover) or whether one of the two assemblages was a subset of the other (nestedness).

In the following, all values are represented as mean ± standard error; statistical analyses were performed using R 3. 5.2 [99].

3. Results

3.1. Species Distribution

The comparison between the results provided by each scenario and the corresponding baseline showed high variability in the species response to the temperature increase.

The majority of species showed no variation (42.1 ± 1.0%), some others displayed an increase (30.7 ± 0.7%), and others a decrease in the number of occupied plots (27.1 ± 0.9%). Considering only endemic and vulnerable species, the percentage of species with decreasing occupancy was higher (39.8 ± 1.1% for endemics; 57.4 ± 2.5% for vulnerable species) and in the case of the vulnerable species, the percentage of increasing species also became lower (7.5 ± 0.7%). Considering only varying species, both the endemic and the vulnerable species showed a higher percentage of decreasing occupancies and a lower percentage of increasing species, when compared to the whole set of species (Figure 2).

For the species that showed a variation, the amount of change was usually low (1st quartile = −4, median = 1, 3rd quartile = +4), even if the values ranged between −28 and +21 plots.

Comparing the three model classes, we found that the number of species with changing distributions was lower with a larger number of environmental constraints. Indeed, varying species in each scenario significantly differed between model classes (TRV = 25.7 ± 4.4, TR = 43.7 ± 4.4, T = 67.7 ± 8.4; 20.64 < χ2 < 27.82, p < 0.001).

Comparing the three scenarios, we observed a higher number of varying species in the Min scenario (Max = 37.3 ± 10.7, Min = 56.7 ± 14.6, D = 43.0 ± 11.3) in each model class, but the differences were significant only in the T model class (χ² = 8.56, p = 0.02).

3.2. Species Richness

Changes in species richness showed considerable variations across taxonomic groups, warming scenarios, and vegetation belts, even if the amount of change was generally low, as exemplified by the effect size, ranging from −0.4 to 0.3 (Figure 3).

Butterflies displayed an increase in species richness per plot for almost all models and warming scenarios (Figure 3e); the impact of temperature increase is significantly stronger on this taxon than for all taxa pooled together (

Table 1. Response of different taxa. Results of effect size analysis across taxa by linear mixed effect models (DF = 3281, the plot identity is taken as a random effect). The baseline is provided by all taxa pooled together. The estimates (and the standard errors) of the fixed effects are reported, followed by t-values and p-values. The ANOVA test of the global effect of the factor variable “taxon” is significant (F = 7.768, p < 0.0001). Significant differences of single taxa are shown in bold.

Response of Different Taxa	Estimate	t-Value	p-Value
Intercept	0.004 (0.005)	0.926	0.354
Birds	0.008 (0.003)	2.292	0.022
Butterflies	0.027 (0.003)	7.778	0.000
Carabids	−0.001 (0.003)	−0.201	0.841
Spiders	−0.028 (0.003)	−7.965	0.000
Staphylinids	0.005 (0.003)	1.420	0.156

).

In the case of spiders, the results showed the opposite trend, with a decrease in species richness in all cases (Figure 3d) and significantly lower effect size than for all taxa pooled together (Table 1). Carabids showed no significant differences in any of the scenarios (Figure 3b). Staphylinids (Figure 3c) and birds (Figure 3f) showed slightly different results when considering different model classes.

The overall species richness showed a significant increase considering only the TRV model class.

For all taxa except carabids, we observed significant differences in the species richness changes predicted by the different model classes (

Table 2. Differences across model classes. Results of effect size analysis across model classes by linear mixed effect models, for all taxa pooled together and for each taxon separately (DF = 494, plot identity as a random effect). The baseline is the T model. The estimates (and the standard errors) of the fixed effects are reported, followed by t-values and p-values. The ANOVA test of the global effect of the factor variable “model class” is reported in the lower two rows (F-value and p-value). Significant values are indicated in bold.

Differences Across Model Classes	All Taxa	Carabids	Staphylinids	Spiders	Butterflies	Birds
Intercept	0.003 (0.004)	0.012 (0.008)	0.000 (0.007)	−0.030 (0.005)	0.047 (0.007)	−0.009 (0.007)
t-value	0.673	1.5	0.082	−5.750	6.299	−1.291
p-value	0.501	0.134	0.935	<0.0001	<0.0001	0.197
TR	−0.004 (0.002)	−0.012 (0.007)	−0.000 (0.006)	0.006 (0.003)	−0.033 (0.005)	0.024 (0.005)
t-value	−1.768	−1.629	−0.023	1.729	−6.865	4.723
p-value	0.078	0.104	0.981	0.084	<0.0001	<0.0001
TRV	0.009 (0.002)	−0.013 (0.007)	0.026 (0.006)	0.012 (0.003)	−0.014 (0.005)	0.039 (0.005)
t-value	3.622	−1.736	4.429	3.213	−3.009	7.573
p-value	0.0003	0.083	<0.0001	0.001	0.003	<0.0001
F-value	15.099	1.893	13.15	5.173	23.684	29.259
p-value	<0.0001	0.152	<0.0001	0.006	<0.0001	<0.0001

). For butterflies, T models indicate an increase in species richness, which is less pronounced for TR and TRV models. For spiders, T models indicate a more pronounced decrease in species richness than TRV models. For all taxa pooled together, staphylinids and birds, TRV models indicate a significantly larger effect than T and TR models.

The analysis of changes in species richness along the altitudinal gradient indicated a more pronounced effect at higher altitude, in particular above 2000 m (Figure 4). In the Alpine belt, we observed a significant effect for all taxa pooled together, for butterflies, for birds and for spiders (

Table 3. Differences across vegetation belts. Results of the analysis of effect size across vegetation belts by linear mixed effect models, for all taxa pooled together and for each taxon separately (DF = 496/59/59, plot identity included as a random effect). The baseline is the montane belt. The estimates (and the standard errors) of the fixed effects are reported, followed by t-values and p-values. The ANOVA test of the global effect of the variable “model class” is also reported in the lower two rows (F-value and p-value). Significant values are indicated in bold.

Differences Across Vegetation Belts	All taxa	Carabids	Staphylinids	Spiders	Butterflies	Birds
Intercept	−0.007 (0.006)	0.008 (0.012)	−0.001 (0.010)	−0.031 (0.008)	0.007 (0.010)	−0.007 (0.011)
t-value	−1.076	0.682	−0.139	−3.862	0.661	−0.701
p-value	0.283	0.496	0.890	0.0001	0.509	0.484
Subalpine	−0.003 (0.009)	−0.021 (0.017)	0.008 (0.015)	−0.007 (0.011)	−0.002 (0.015)	0.018 (0.015)
t-value	−0.292	−1.212	0.555	−0.649	−0.128	1.172
p-value	0.771	0.230	0.581	0.519	0.899	0.246
Alpine	0.029 (0.008)	0.004 (0.016)	0.020 (0.014)	0.023 (0.011)	0.062 (0.014)	0.036 (0.014)
t-value	3.503	0.228	1.505	2.224	4.428	2.564
p-value	0.001	0.821	0.138	0.030	<0.0001	0.013
F-value	9.285	1.252	1.170	4.759	13.972	3.308
p-value	0.0003	0.293	0.317	0.012	<0.0001	0.043

). In the case of spiders, this implies that the decrease in species richness observed in the montane and the subalpine belt is less pronounced at higher altitudes (Figure 4; Table 3).

3.3. Community Composition

In almost all cases, warming scenarios determined significant changes in community composition along the first CA axis. The Wilcoxon test was applied to assess changes in plot scores from current to warmer conditions; it is significant in 53 over 54 cases, p < 0.01, the only exception are staphylinids in the TRV model class, max scenario. The pattern along the second axis is less coherent (Wilcoxon test with p < 0.01 in 35 cases over 54).

For both the current conditions and the warming scenarios, plot community composition significantly changed from the montane to the alpine belt (Figure 5a; current, KW Test, df = 2, all p < 0.0001, Figure 5b; projected, KW Test, df = 2, all p < 0.0001; Figure 5c).

Plot scores were significantly correlated in all cases (Spearman’s rank correlation; first axis, 0.978 < ρ < 0.998; second axis, 0.894 < ρ < 0.997), indicating that species composition changed gradually and coherently along the altitudinal gradient, retaining similar patterns of species dissimilarity.

Differences in the first axis between pairs of plots (current minus projected), if significant, showed that the amount of change in community composition was more pronounced for the montane belt compared to the others (Figure 5d). Interestingly, we observed significant differences for all taxa pooled together (for all model classes and scenarios), and in many cases also for carabids (T and TR, all scenarios), staphylinids (only T, Min and D scenarios), and butterflies (T, all scenarios; TR, Min and D scenarios; TRV, Max scenario), always showing the lowest values (around 0) for the montane belt.

As an estimate of community homogenization, we estimated the mean Euclidean distance from the distribution centroid for current conditions and for warming scenarios, showing always lower values in the warming case (which were also significant in 29 to 54 cases, p < 0.05).

The Jaccard Index values were close to zero for all model classes and scenarios. Considering all taxa pooled together, we observed significant differences across scenarios (p < 0.0001), with the lowest values in Max, and across model classes (p < 0.0001), with the highest values in T.

For all taxa pooled together (Figure 6a), the variation in community composition in the warming scenarios is dominated by turnover (β_sim = 0.066 ± 0.005). The nestedness component is lower (β_sne = 0.026 ± 0.002), indicating no clear pattern of species loss or gain. This global pattern has been observed also for carabids and spiders (Figure 6b,d), while in the case of butterflies and birds, differences are less pronounced and there is a tendency toward higher values of nestedness (Figure 6e,f). Staphylinids show no clear pattern (Figure 6c).

4. Discussion

Disentangling the role of temperature increase in shaping biodiversity patterns, especially in mountain ecosystems, is a fundamental challenge of conservation biology [2,3,8]. Species respond at individual level [1,2] and there is a current knowledge gap for many taxonomic groups [31,32,33,49]. Projections of potential future species distributions along elevational gradients need high-resolution spatial data [83,100]. In our work, the Maxent approach allowed to analyze data collected at a fine spatial scale (plots with 100 m radius) and to take into account the response of rare and localized species. Reliable projections should also consider that ecotypes from different geographical origins can have different responses to the same climatic variations (e.g., [101]). Our approach reduced this bias in two ways. First, we used the same set of plots to derive the empirical relationships between species and climatic/environmental variables and to project future potential distributions. Second, the limited geographic extent (northwestern Italian Alps) of the sampled plots assures that different populations of the same species will respond coherently.

In general, our simulations indicated that even for a moderate temperature increase, species richness and community composition can display subtle but significant alterations. Since different studies suggested different roles of minimum and maximum temperature increase [14,85], we simulated three scenarios having different relative importance of the minimum and maximum daily temperature rise. We also adopted three different empirical biodiversity model classes, characterized by the presence of temperature alone or of other environmental variables (vegetation and geographical location) in addition to temperature. Differences between scenarios were not as marked as those emerging from the different model classes.

We detected a high variability in the response of the different species, depending on the position and breadth of the climatic niche and the species range. The majority of modeled species did not show much response to the temperature increase. On the other hand, high-altitude species were the most affected, consistently with the fact that habitat specialists are negatively influenced by environmental change (e.g., [7]). The probability of extinction under climate change reflects the species ability to shift with suitable habitats; a scarce ability to withstand climate change impacts could depend also on a reduced niche breadth [7,102,103].

In our analysis, endemic species appeared as less affected than vulnerable ones. In the altitudinal gradient explored here, the spatial distribution of endemic species is spread over the whole altitudinal and temperature range, owing to the presence of species that are characteristic of the montane belt. Our results indicate that many of the endemic species studied here are not restricted to high altitudes and can shift their occupancy area along the altitudinal gradient. These results are not in disagreement with the general vulnerability of endemic species to climate change (e.g., [7,104]): thanks to the small scale of the analysis, our study indicates how local altitudinal gradients can offer a possibility of survival for otherwise declining species.

Significant differences between taxonomic groups were also observed.

Carabids did not show any relationship with temperature [52] and, as a consequence, they did not display any significant change in species richness under temperature increase scenarios. On the other hand, they displayed the highest “temporal dissimilarity” in community composition (even if always low, with an overall Jaccard index < 0.15). This behavior was mainly determined by turnover, i.e., by the substitution of groups of species under climate warming.

In the case of staphylinids, the pattern of change was also not marked: few significant differences were observed, and they were not coherent across model classes and scenarios, hampering the possibility of drawing general considerations.

Butterflies are heliothermic and can be considered climatically sensitive [33,105]. In our analysis, for many species the distribution was mainly determined by temperature. This taxon displayed the highest increase in species richness under temperature increase. This pattern is exacerbated by the model classes using temperature only as a constraint, and it is particularly strong in the alpine belt. These results indicate that butterflies can respond quickly to warming, determining new species assemblages, in particular at high altitudes and that the vegetation structure has the potential to buffer this response. Similar results were found analyzing observed temporal changes in butterfly distribution and community composition [106,107], also during a short time frame [108] and in the same mountain ranges considered here [53].

On the opposite side, spiders strongly suffered in temperature increase scenarios, showing a decrease in species richness. In our monitoring program, spiders were highly localized and many species were present only in a small number of plots. Therefore, variations of microclimatic conditions can effectively reduce their areas of occurrence, determining a general decrease in species richness. Models using only temperature provide the most marked changes, indicating again that the vegetation structure can buffer the effects of warming. For spiders, the highest decrease in species richness was observed in the montane belt. Our results suggest that the montane belt will become too warm for many species. Indeed, spider species are usually influenced by the small-scale vegetation structure and are adapted to the local environmental conditions [109,110]. This leads to a high level of turnover along altitudinal gradients and consequently to assemblages that are sensitive to small changes in local climatic conditions, as observed here.

Birds display a significant increase in species richness for many of the tested models and scenarios and a general positive effect of warming. Temporal analysis of bird and butterfly communities showed that both taxa respond to warming, with a faster response of butterflies [106]. For birds, changes in community composition were dominated by the nestedness component, indicating that species losses and gains represent the principal pattern of change.

Considering “temporal” dissimilarity as a whole, we observed significant changes for all cases (model classes, scenarios, taxa). Plot assemblage composition coherently changed in warming scenarios, showing a tendency to homogenization as observed in other studies [111,112,113].

As in other studies focused on plants [11], our warming scenarios indicated an upward shift from the lower to the higher belts, implying an increase in species richness at the highest altitude. Alpine meadows had a high probability of experiencing the largest modifications in community composition, as we coherently observed for almost all taxa. Such changes in richness and composition, if effectively realized, will probably be detrimental to biodiversity: the colonization of the alpine belt by species from lower altitudes could increase competition, with negative effects for localized and specialized taxa [13,114,115]. Consequently, the increase in species richness could be transitory, followed by a decline of strictly alpine species stressed by the growing competition and on the verge of going beyond their tolerance breadth [24,116].

The montane belt, on the opposite, is the lower limit of the altitudinal gradient considered here and we cannot account for the colonization of species coming from still lower elevations. We classified plots along the vegetation belts considering both altitude and potential vegetation [65,66]: consequently, some montane plots display temperature values that are lower than for subalpine ones, allowing for a partial increase in species richness under climate warming. In any case, the montane belt community composition was the most stable, showing lower dissimilarity than the upper two belts.

Among the different warming scenarios, the one with a larger increase of minimum temperature displayed the largest number of varying species, indicating that the minimum daily temperature can be an important limiting factor for species distribution [65,66]. Considering that the minimum temperature has increased at a faster rate than the maximum temperature during the latter half of the 20th century [117], variations in biodiversity patterns larger than the ones simulated here could occur.

Even if uncertainty is the rule in species distribution models [40,45], projections are useful to explore biodiversity patterns, identifying specific areas or groups of species that should be monitored. In this regard, the different responses displayed by different taxa confirm the importance of using a multi-taxa approach to estimate climate change effects on animal biodiversity. It is important to include taxa with potentially opposite responses to climate warming, such as butterflies and spiders. Targeted long-term field data are then essential to revisit and fine-tune estimates of biodiversity responses, comparing them with real changes in species responses, and adopting timely conservation strategies [118].

Even if the alpine belt displayed the highest projected change in species richness, monitoring focused only on high altitudes could miss some of the important changes along the altitude gradient. The high level of turnover due to species shifts from the montane and to the alpine belts could happen also in the subalpine belt, and the montane belt can be influenced by the arrival of species from the surrounding lowlands.

In conclusion, our analysis suggested that moderate warming has the potential to change biodiversity patterns in mountain ecosystems, with significant differences between taxa and along the altitudinal gradient. Such changes could determine new ecological relationships, which could in turn influence ecosystem processes in an unpredictable way.