Greater topoclimatic control of above- versus below-ground communities.

Assessing the degree to which climate explains the spatial distributions of different taxonomic and functional groups is essential for anticipating the effects of climate change on ecosystems. Most effort so far has focused on above-ground organisms, which offer only a partial view on the response of biodiversity to environmental gradients. Here including both above- and below-ground organisms, we quantified the degree of topoclimatic control on the occurrence patterns of >1,500 taxa and phylotypes along a c. 3,000 m elevation gradient, by fitting species distribution models. Higher model performances for animals and plants than for soil microbes (fungi, bacteria and protists) suggest that the direct influence of topoclimate is stronger on above-ground species than on below-ground microorganisms. Accordingly, direct climate change effects are predicted to be stronger for above-ground than for below-ground taxa, whereas factors expressing local soil microclimate and geochemistry are likely more important to explain and forecast the occurrence patterns of soil microbiota. Detailed mapping and future scenarios of soil microclimate and microhabitats, together with comparative studies of interacting and ecologically dependent above- and below-ground biota, are thus needed to understand and realistically forecast the future distribution of ecosystems.

INTRODUCTION

Climate defines the spatial distribution of a wide range of organisms, and the need to understand this relationship is particularly revived by the current threat of climate change on ecosystems (Illan, Gutierrez, & Wilson, 2010; Lenoir & Svenning, 2015; Miller et al., 2018; Pacifici et al., 2015; Parmesan & Yohe, 2003; Scherrer & Körner, 2011; Steidinger et al., 2019; Woodward & Williams, 1987). In addition to their recognized ecological effects, the prevalent use of average annual, seasonal or monthly air temperature‐ and precipitation‐related factors in many studies is common because they are freely available (e.g. CHELSA; Karger et al., 2017), with existing robust future scenarios (IPCC, 2014). Such macroclimate data are often associated with topographic variables that are easily derived from digital elevation models (DEM; Amatulli et al., 2018) to model topoclimate and to describe topography‐related habitat characteristics concerning, for example, the distribution of energy (Wang, Qiu, Wang, Wang, & Liu, 2014) and water (Moeslund et al., 2013).

Topoclimatic factors, albeit not capturing microclimatic variations at very fine spatial resolutions (Boulangeat, Gravel, & Thuiller, 2012; Lembrechts, Lenoir, et al., 2019; Mod, Scherrer, Luoto, & Guisan, 2016; Niittynen, Heikkinen, & Luoto, 2018; Scherrer & Guisan, 2019), are important drivers of species distributions. Yet, how their importance varies among taxonomic groups is largely unknown. Most studies on the effects of topoclimate on the distributional patterns of organisms are strongly biased towards a few taxonomic groups (Table 1; Kharouba, McCune, Thuiller, & Huntley, 2013). For instance, plants and some iconic groups of vertebrates and insects have been extensively studied with varying methodologies ranging from basic observational approaches to complex frameworks (Bateman, VanDerWal, Williams, & Johnson, 2012; Bradie & Leung, 2017; Illan et al., 2010; Lawrence et al., 2014; Miller et al., 2018; Mod & Luoto, 2016; Roberts, Nielsen, & Stenhouse, 2014; Seoane, Bustamante, & Diaz‐Delgado, 2004; Staniczenko, Sivasubramaniam, Suttle, & Pearson, 2017), but comparatively fewer attempts have been made to assess the drivers of distributions of soil microorganisms (Bradie & Leung, 2017; Lenoir et al., 2020; Lenoir & Svenning, 2015; Pacifici et al., 2015)—likely due to previous methodological limitations (Riesenfeld, Schloss, & Handelsman, 2004) and, for some groups (e.g. soil protists), to a lower sampling effort (Caron, Worden, Countway, Demir, & Heidelberg, 2008; Geisen et al., 2017, 2018; Seppey et al., 2020; Wilkinson, 1998). Existing studies have evidenced the importance of soil characteristics, such as pH, nutrient content and moisture availability, for explaining the diversity, biomass and community structure of microorganisms in soils (Bahram et al., 2018; Bates et al., 2013; Fierer & Jackson, 2006; Serna‐Chavez, Fierer, & van Bodegom, 2013; Tedersoo et al., 2014). Furthermore, soil temperature has been shown to be decoupled from air temperature (e.g. Scherrer & Körner, 2011) and to better explain the fine‐scale distributions of plants, especially those with small and close‐to‐ground growth form (Lembrechts, Lenoir, et al., 2019). Soil microclimate might thus be more important than air macroclimate in explaining the distribution of organisms strongly associated to soil conditions. While a direct quantification of local soil conditions could allow to test their role, an indirect approach is to compare the ability of topoclimatic conditions to explain and predict the distribution of organisms across taxonomic groups with contrasted life histories (Kharouba et al., 2013).

TABLE 1

Number of results for ISI Web of Knowledge querya of “(climat* AND (distribution OR occurrence) AND ([Taxonomic group])”

Taxonomic group	Number of query results
plant* OR tree* OR vegetation	27,593
bird*	2,829
insect*	2,472
mammal*	2,107
lepidoptera* OR butterfl* OR moth*	1,353
reptil* OR lizard* OR turtle* OR tortoise* OR terrapin* OR snake* OR croco*	1,247
amphibia* OR frog* OR toad* OR salamander* OR newt*	1,212
fungi OR fungus	1,139
bacteria*	896
“micro*organism*” OR “micro*biota”	448
orthoptera* OR grasshopper* OR cricket*	177
bombus OR bumblebee* OR bumble*bee*	78
micro*eukaryot* OR protist*	72

Searched on 11 July 2019.

Due to the tight direct and indirect couplings among different taxonomic and functional groups, such as biotic interactions and nutrient cycling (Seibold, Cadotte, MacIvor, Thorn, & Muller, 2018; Wardle et al., 2004), understanding the role of climate control on multiple groups is of crucial importance to more realistically forecast the future of biodiversity and entire ecosystems (Cavicchioli et al., 2019; Guo, Feng, et al., 2018; Hagedorn, Gavazov, & Alexander, 2019; Schleuning et al., 2016). In addition to the taxonomic bias, assessing and comparing the role of climate and warming across groups is hindered by the varying sampling designs, study areas and spatial scales used among the existing studies (Rodríguez‐Castañeda, Hof, Jansson, & Harding, 2012). A robust comparative quantification of the role of climatic control on distributions across the various above‐ and below‐ground groups composing ecosystems was thus needed.

Here using comprehensive taxon occurrence data from the same mountainous study area in western Switzerland and a robust modelling, evaluation and prediction framework, we aim to: (a) assess to what extent the commonly used and readily available topoclimatic variables explain alone the distribution patterns of six above‐ground taxonomic groups (amphibians, reptiles, grasshoppers, butterflies, bumblebees, vascular plants) and three below‐ground ones (soil fungi, bacteria and protists); and (b) make baseline predictions and comparisons about how those taxa may respond to climate change. For this, we fitted climate suitability models with the occurrences of each taxon and six topoclimatic variables and compared models’ predictive ability and the taxa's spatial predicted probability of occurrence (PPO) under current and future climatic conditions.

MATERIALS AND METHODS

Study area

The study area is a priority area for transdisciplinary research and biodiversity conservation located in the Western Swiss Alps (46°10′–46°30′ N; 6°60′–7°10′ E). It covers an area of c. 700 km2 with elevations ranging from 370 to 3,200 m, resulting in heterogeneous climatic and topographic conditions. The annual mean temperature and precipitation sum vary from 8°C and 1,200 mm at 600 m a.s.l. to −5°C and 2,600 mm at 3,000 m a.s.l., respectively (Bouët, 1985). The distribution of energy, water and debris strongly varies with slope and aspect. Landcover is dominated by alpine grasslands, forests, glaciers and agricultural lands. More detailed descriptions of the study area are provided in Randin et al. (2006), Dubuis et al. (2011) and at www.unil.ch/centre‐montagne.

Data

The data sets used to build the species distribution models (SDMs) and to project habitat suitability under current and future climatic conditions consist of spatially explicit information on biota and topoclimatic conditions (for detailed description of the data, see Appendix S1). The data on biota comprise site‐level presence–absence information about species, genera and/or orders from nine taxonomic groups: amphibians, reptiles, grasshoppers, butterflies, bumblebees, plants, fungi, bacteria and protists (Table 2; Figures S1–S10 in Appendix S1). Above‐ground species were visually identified based on their phenotype. Amphibian and reptile data sets resulted from point observations, whereas sampling of plants, insects and soil (for amplicon sequencing of microorganisms) was targeted on non‐forested sites selected following a random stratified sampling design considering elevation, slope and aspect. Occurrences of soil bacteria, fungi and protists were obtained by amplicon sequencing of the V5 region of 16S, ITS1 and the V4 of the 18S, respectively, region of the ribosomal RNA small subunit gene. About 38,000 operational taxonomic units (OTUs) were taxonomically identified but as this number exceeded the computational power available (especially for postprocessing of spatial projections), we clustered them at genera and order levels.

TABLE 2

Data sets of the nine taxonomic groups available and used for modelling in the same study area in the Western Swiss Alps

Taxonomic group	n of taxa (total)	n of modelled taxa	n of sites	Prevalence of taxa (mean ± SD)	Niche breadth of taxa a (mean ± SD)	Data origin
Amphibians (sp.)	14	5	133	0.40 ± 0.28	1.71 ± 0.78	Schmidt and Zumbach (2019)
Reptiles (sp.)	12	12	1,144	0.09 ± 0.08	1.59 ± 0.35	Pittet (2017)
Grasshoppers (sp.)	41	21	202	0.25 ± 0.19	1.98 ± 0.31	Pradervand et al. (2013)
Butterflies (sp.)	140	78	208	0.22 ± 0.14	2.14 ± 0.29	Pellissier et al. (2012)
Bumblebees (sp.)	29	20	202	0.25 ± 0.16	2.24 ± 0.21	Pellissier et al. (2013)
Plants (sp.)	795	296	909	0.08 ± 0.08	2.02 ± 0.32	Dubuis et al. (2011)
Plants (ge)	288	160	909	0.14 ± 0.14	2.23 ± 0.36	Dubuis et al. (2011)
Fungi (ge)	190	92	103	0.34 ± 0.17	2.98 ± 0.36	Pinto‐Figueroa et al. (2019)
Fungi (or)	36	10	198	0.40 ± 0.28	3.03 ± 0.11	Pagni et al. (2013)
Bacteria (ge)	758	346	258	0.38 ± 0.29	2.92 ± 0.46	Yashiro et al. (2016)
Bacteria (or)	276	124	258	0.46 ± 0.29	3.04 ± 0.48	Yashiro et al. (2016)
Protists (ge)	496	285	220	0.38 ± 0.26	2.86 ± 0.43	Seppey et al. (2020)
Protists (or)	161	92	220	0.47 ± 0.31	2.91 ± 0.36	Seppey et al. (2020)

Abbreviations: ge, genera; or, orders; sp., species.

Unitless; calculated from variation in principal components of environmental conditions occupied; see Section 2.

The topoclimatic variables used were mean annual temperature, annual temperature range, sum of annual precipitation, topographic position index, slope and potential annual solar radiation (see Figures S11–S20 in Appendix S1). These were derived from climate data and DEM (Federal Office of Topography: swisstopo.ch) as spatial layers at a resolution of 25 m, with climatic data based on both the current conditions (the most recent climate normal: 1981–2010) and the future scenarios (IPCC's SRES A1b and A2 for the periods 2020–2050, 2045–2075 and 2070–2100; Zubler et al., 2014; meteoswiss.ch; and Figure S21 in Appendix S1). As the new Representative Concentration Pathways (RCPs; IPCC, 2014) are not yet available for the study area with the needed resolution, and as we do not aim to make realistic future predictions, but to test climate control among taxonomic groups, we consider these scenarios suitable for this study. For comparison between SRES and RCPs, see Rogelj, Meinshausen, and Knutti (2012).

Data preparation for analyses

To test the predictive power of the topoclimatic variables, we evaluated the models using a split‐sampling approach (10×) assigning 80% of the sites to model calibration and the remaining 20% to model evaluation so that each site was used eight times for calibration and twice for evaluation. The two predicted PPO values per evaluation site were then averaged to make a comparison against the observed occurrences.

Approximately half of the taxa with less than eight presences or absences within the calibration data sets and with less than one presence or absence within the evaluation data sets were removed to allow confident parameter estimation and model evaluation (Table 2). To test if the taxon's commonness or tolerance of varying environmental conditions within the study area are related to the magnitude and direction of predicted future changes, we computed prevalence (as a proportion of presence observations within sites; Figure S22 in Appendix S1) and niche breadth per taxon (Figure S23 in Appendix S1). To obtain an estimate of the niche breadth, we first applied a principal component analysis (PCA; using dudi.pca from the R‐package ade4) to the matrix of the six environmental variables over all pixels of the study area. Niche breadth was then calculated as the sum of the standard deviations of the scores of the first and second PCA axes where the taxon is present (Segurado et al., 2011).

Modelling

To allow consideration of taxa with extreme (i.e. very low or high) prevalence while still being able to incorporate all six environmental variables without over‐parametrizing the models, we used Ensemble of Small Models (ESM; Breiner, Guisan, Bergamini, & Nobis, 2015) with the corresponding functions in the R‐package ecospat (Di Cola et al., 2017). ESM first fit ‘small’ models by forming all possible bivariate combinations (i.e. two at a time) of the environmental predictors (here, ). The bivariate models are then combined to create the final ensemble model. Following the recommendations of Breiner, Nobis, Bergamini, and Guisan (2018), we fitted the bivariate models using the Generalized Boosting Model (GBM), Generalized Linear Model (GLM), Artificial Neural Network and Classification Tree Analysis. The final ensemble for each taxon was built based on the converged bivariate models weighted based on their model performance as measured by AUC (bivariate models with AUC < 0.5 were discarded from the final ensemble models). The default settings of the functions were used, except that all observations, regardless of taxon prevalence, had the same weight.

Evaluating the predictive power of topoclimatic variables across taxa

The per taxon models were evaluated with AUC (i.e. the area under the ROC curve), the maximum True Skills Statistic (maxTSS) and the maximum Cohen's KAPPA (maxKAPPA), which do not require converting the predicted PPO to presence/absence predictions (Guisan, Thuiller, & Zimmermann, 2017), as such thresholding is shown to influence the modelling outcome (Benito, Cayuela, & Albuquerque, 2013; Calabrese, Certain, Kraan, & Dormann, 2014; Liu, Berry, Dawson, & Pearson, 2005). In addition, by using multiple metrics with varying prevalence dependencies (Lawson, Hodgson, Wilson, & Richards, 2014; Lobo, Jiménez‐Valverde, & Real, 2008; Somodi, Lepesi, & Botta‐Dukat, 2017), we avoided the prevalence of taxa affecting the interpretations of the role of topoclimatic variables. The predictor contributions (sensu variable importance in biomod2—R‐package; Thuiller, Georges, Engler, & Breiner, 2016) of environmental variables per taxon and per technique were defined as the sum of weights (based on AUC) of bivariate models including the variable of interest divided by the sum of all weights. Variable contributions in the final ensemble models were then further calculated as the weighted means of these contributions with weighting based on the overall AUC of the technique.

Projections under current and future conditions

PPO of each taxon were obtained across the entire study area (i.e. for all of the 25 m × 25 m pixels defined by topoclimatic data) under current and future environmental conditions based on a full model trained with all data (i.e. outside CV used for evaluation). Some areas under current and future scenarios are not analogous to the environmental conditions present at the sampled sites (Figures S1–S10, S24 and S25 in Appendix S1). This implied that the models needed to extrapolate predictions to the environmental space outside the training and evaluation conditions of the model (Fitzpatrick & Hargrove, 2009). Thus, the changes between current and future predictions for each taxon were also quantified in the subset of the study area with analogous environmental conditions to the sites sampled under current and all future scenarios (i.e. the same mask for all taxa of a group was used to exclude non‐analogous areas from the current and all future projections; Figure S25 in Appendix S1). Furthermore, forested areas were masked from the spatial projections of open‐habitat insects, plants and soil microorganisms. The same forest mask (from Swiss Map Vector 25 BETA; swisstopo.ch) was used for the current and all future projections, which simulates a scenario where the forest cover in the study area would overall remain constant because the Swiss forest policy aims to prevent forest encroachment of pastures (following the ‘lowered agriculture production’ socio‐economic scenario of Gago‐Silva, Ray, & Lehmann, 2017). Other scenarios are, however, plausible too and we acknowledge that incorporating land use changes could affect the predictions (Guo, Lenoir, & Bonebrake, 2018). For example, forest expansion (following the ‘liberalization’ scenario of Gago‐Silva et al., 2017) could further restrict the future distributions of grassland taxa (Tasser & Tappeiner, 2002).

The examination of predictions and the changes therein was based directly on the PPO values under current and future conditions, thus avoiding the problematic use of binarization thresholds (Nenzen & Araujo, 2011). To estimate the expected magnitude of change for each taxon in the study area, we calculated the overlap of the PPO values between current and future projections. For each taxon, the overlap was defined based on the sum of differences between current and future PPO values across all pixels, that is, the Schoener's D index (here, derived using the ‘modOverlap’ function in the fuzzySim R‐package; Barbosa, 2015):where p (resp. p) denotes a taxon's PPO in pixel i proportional to the sum of taxon's PPOs in all pixels in time X (resp. Y; Barbosa, 2015; Warren, Glor, & Turelli, 2008):where q is taxon's PPO and n is number of pixels. The index varies between 1 (i.e. identical current and future predicted distributions) and 0 (i.e. completely non‐overlapping current and future predicted distributions), and thus the lower the overlap, the higher the predicted future change in taxon's distribution.

The spatial variation in the magnitude of community changes within each taxonomic group was assessed by first calculating, for each 25 m × 25 m pixel of the study area, the Euclidian distance between the current and future PPO values of all taxa, using the formula for n‐dimensional space:where p is the PPO of taxon i in the pixel based on the prediction under current conditions, q is the PPO of taxon i the pixel based on the prediction under future conditions and n is the number of taxa in the respective taxonomic group. These Euclidian distances per group per pixel per future scenario were then divided by the highest possible distance per group (; i.e. the case where for all taxa, PPOs would change from 0 to 1 or from 1 to 0, with resulting distances for all taxa being 1). The values representing the magnitude of community change in each pixel then vary between 0 and 1, where a value of 0 means that the PPO of none of the taxa of the respective taxonomic groups changed between current and future prediction, and a value of 1 means that the PPO of all taxa changed from 1 to 0 or from 0 to 1.

Finally, the direction of change (i.e. estimate of gains vs. losses in suitable habitat) was examined for each taxon and future scenario by calculating the mean and median changes in PPOs across the study area and separately at low, mid and high elevations.

Comparing ESMs to standard SDMs

The standard ensemble SDMs (with commonly used GBM, GLM, Generalized Additive Models = GAM and Random Forest = RF weighted by AUC within Biomod2‐platform; Hao, Elith, Guillera‐Arroita, & Lahoz‐Monfort, 2019; Thuiller et al., 2016; Thuiller, Lafourcade, Engler, & Araújo, 2009) were fitted for the subsets of species and taxonomic groups that included at least 50 presences and absences in each CV‐fold of the training data and at least one presence and absence in each fold of the evaluation data. A selection of model evaluations and analyses of the predicted changes were performed to compare the ESMs and standard SDMs.

RESULTS

Across all taxa and techniques, the mean ± standard deviation of the model performance is 0.747 ± 0.096 (measured by AUC; for maxTSS and maxKAPPA see Figure S26 in Appendix S2). The model performance of above‐ground species (0.823 ± 0.071) is higher than the model performance of below‐ground microorganisms (0.700 ± 0.079 for genera and 0.708 ± 0.086 for orders; Figure 1a; Figure S26 in Appendix S2). The higher model performances are statistically significant when tested with Mann–Whitney U test (U = 271,431, p < .001 when 432 above‐ground species are compared to 723 below‐ground genera and U = 81,934, p < .001 when compared to 226 below‐ground orders) and the probabilities of AUC were found to be higher for above‐ground species than for below‐ground genus or order are 0.87 and 0.83, respectively, based on the common language effect size (McGraw & Wong, 1992). For all taxonomic groups, the annual mean temperature is the most and solar radiation the least influential topoclimatic variable (Figure S27 in Appendix S2). The contributions of the other four predictors do not rank consistently among the taxonomic groups.

FIGURE 1

The role of topoclimate and climate change in defining current and future taxa distributions. (a) Performance as measured with AUC per taxon of the topoclimatic ensemble models. The model performance of above‐ground species is significantly better than the model performance of below‐ground microorganisms (p < .001; measured with Wilcoxon rank sum test). For other evaluation metrics and per technique (Generalized Linear Model, Generalized Boosting Model, Artificial Neural Network, Classification Tree Analysis), and for plant genera and soil microorganism orders, see Figure S26 in Appendix S2. (b) Spatial overlap in predicted probability of occurrence (PPO) between current and future (A2 2085 scenario) projections within the whole study area. For other scenarios, plant genera, soil microorganism orders and environmentally analogous areas, see Figure S28 in Appendix S2. Boxes in boxplots span the 25th–75th quartile, with median (black bar) and mean (orange point) in the middle. Whiskers span the lowest and highest scores, yet in maximum to 1.5*(75th–25th quartile); outlier scores are indicated by black dots

Changes between current and future predicted PPO increase linearly with the severity of the deployed climate scenarios, and the changes are mostly similar between the whole study area and the parts of the study area which exclude novel future environmental conditions (i.e. non‐analogous conditions; see Appendix S2). Forecasts show a relatively smaller overlap between the current and future PPO for reptiles, insects and plants, whereas higher overlap indicates less change especially for the soil microbial distributions (Figure 1b; Figure S28 in Appendix S2). The magnitude of overlap is positively correlated with the prevalence and niche breadth of the taxa (see Tables S1 and S2 and Figures S29 and S30 in Appendix S2).

The strongest predicted changes in taxonomic composition across the study area, as measured by the Euclidean distance in the ordination space of the PPO values between the current and future conditions, occur at the highest elevations (Figure 2). This holds for all taxonomic groups, except for amphibians, reptiles and bumblebees, for which the strongest changes are predicted at midelevations (Figure 2d,g,h; Figures S31–S42 in Appendix S2). The PCA, restricted to the non‐forested areas, divides the taxonomic groups into five categories (Figure 2j). The future variations in reptile and amphibian compositions strongly differ from the variations of other groups. The change in composition of bumblebees is opposite to that in amphibians. The predicted variation in the compositions of grasshoppers, butterflies, plants and fungi is similar, whereas soil bacteria and protists cluster separately.

FIGURE 2

Spatial variation in the magnitude of the predicted future changes in community structure. (a) Location of the study area in the Swiss Alps is marked by yellow squares. (b–i, k) Magnitude of predicted changes in community structure between current conditions and future scenario A2 for year 2085. The maximum possible magnitude is 1, meaning that predicted probability of occurrences (PPOs) of all taxa are predicted to change from 1 to 0 or from 0 to 1, and the minimum possible magnitude is 0, meaning that PPOs of none of the taxa are predicted to change between current and future prediction. Grey areas mark the forest cover masked from the predictions of taxa, as sampling targeted non‐forested sites. (j) Principal component analysis (PCA) of community structure changes of the nine taxonomic groups in the non‐forested areas shows the (dis)similarity in spatial patterns, that is, groups with lines pointing in the same direction have similar spatial patterns in the magnitude of community changes in non‐forested areas, whereas lines of varying directions indicate varying patterns. (l) Elevation of the study area. For maps and PCAs of other taxonomic ranks, scenarios and years, see Figures S31–S36 in Appendix S2, and for the relationships between the changes and elevation, and for comparison to analogous environmental space, see Figures S37–S42 in Appendix S2

There are slightly more taxa with increases than decreases in mean PPO across the study area between now and the future (Figure 3a; Figure S43 in Appendix S2), especially if only the area with analogous environments is considered (Figure 3b; Figure S43 in Appendix S2). For example, under the A2 scenario for 2070–2100, 83.3% of reptile species and 71.7% of bacteria genera showed an increase in mean PPO. In contrast, amphibians, butterflies and plants show the lowest proportional increase in mean PPO (40.0%, 38.5% and 36.5%, respectively). The mean change in PPO across the study area mainly shows a weak and non‐significant relationship to the prevalence and niche breadth of taxa (Tables S3 and S4 and Figures S44 and S45 in Appendix S2). Overall, the changes in mean PPO are stronger and often more positive at mid‐ and high‐compared to low elevations (Figure 3c–e; Figure S46 in Appendix S2). Median changes in PPOs are similar to mean changes (Figures S47–S49 in Appendix S3).

FIGURE 3

Proportions of taxa with different mean changes in predicted probability of occurrence (PPO; classified to six classes) under A2 scenario for 2085. (a) Across the total study area. (b) Across an environmentally analogous part of the study area. (c–e) At different elevational bands (low = <1,180 m a.s.l.; mid = 1,180–1,650 m a.s.l.; high > 1,650 m a.s.l.). For other scenarios and taxonomic ranks, see Figures S43 and S46 in Appendix S2. For median changes, see Figures S47–S49 in Appendix S3

For the subset of taxa modelled with both ESMs and standard SDMs (i.e. all variables incorporated simultaneously), the model performances are highly correlated (0.964–0.983; Table S5 and Figure S50 in Appendix S3), but in general, ESMs perform better for most taxonomic groups (Table S6 in Appendix S3). The mean relative variable contributions across taxa follow the same order for ESMs and standard SDMs (Figure S51 in Appendix S3). The mean changes in PPO of ESMs and standard SDMs are highly correlated, but predictions of standard SDMs show a larger variation (Figures S52–S55 and Table S7 in Appendix S3).

DISCUSSION

Our results show that topoclimatic variables alone better explain and predict the distributions of above‐ground groups (especially reptiles, grasshoppers and butterflies) than those of below‐ground microorganisms and that the forecasted future changes in soil microorganism distributions are accordingly smaller and more spatially different. The low model performance indicates that other factors, in addition to the easily obtainable topoclimatic variables, are needed to improve models and future predictions of soil microorganisms (Bahram et al., 2018; Bates et al., 2013; Fierer, Strickland, Liptzin, Bradford, & Cleveland, 2009; Lembrechts, Nijs, & Lenoir, 2019; Li et al., 2018). Importantly, Lembrechts, Lenoir, et al. (2019) showed that soil temperature better explains the distributions of close‐to‐ground growing plants than free air climate. This suggests that soil microclimate, which better expresses the temperature and moisture conditions affecting below‐ground and ground‐dwelling organisms and which can be strongly decoupled from air temperature in space and time (see also Bramer et al., 2018; Graae et al., 2012; Scherrer & Körner, 2011; Suggitt et al., 2011), likely also better explains the distributions of soil microorganisms. Furthermore, other edaphic factors (particularly soil pH) and biotic interactions are known to be major drivers for the diversity and community structure of soil fungi and bacteria (Bahram et al., 2018; Bates et al., 2013; Fierer et al., 2009; Li et al., 2018; Malard, Anwar, Jacobsen, & Pearce, 2019; Yashiro et al., 2016), and also but to a lesser extent for the diversity of protists (Seppey et al., 2020). Thus, they can also be expected to strongly contribute to explain the distribution of individual soil microorganisms.

The forecasted lower magnitudes of distributional changes of soil microorganisms, therefore, do not necessarily imply less severe changes to be expected in the future, but a smaller direct role of macroclimatic factors. As the effects of climate change may also be indirect, propagated by alterations, for example, in soil and vegetation (Alexander, Diez, Hart, & Levine, 2016; Collins, Stajich, Weber, Pombubpa, & Diez, 2018; Mod & Luoto, 2016; Yashiro et al., 2018), future scenarios are needed for soil factors in an analogous way to those developed for macroclimate (IPCC, 2014; Lembrechts et al., 2020). Furthermore, the multicollinearity should be assessed between topoclimate and missing predictors (Oliver & Morecroft, 2014; Sears, Raskin, & Angilletta, 2011), because large spatial and temporal offsets have been demonstrated for example between air and soil temperatures (Aalto, Scherrer, Lenoir, Guisan, & Luoto, 2018; Lembrechts, Lenoir, et al., 2019), and thus the used predictors can have a strong influence on the interpretations of future biodiversity changes (Lenoir et al., 2013; Lenoir, Hattab, & Pierre, 2017; Scherrer, Schmid, & Körner, 2011). Deriving and incorporating all influential and change‐prone environmental variables at ecologically relevant scales in the models and future predictions (see e.g. Lembrechts, Nijs, et al., 2019; Mod et al., 2016) would allow to assess if the here demonstrated among‐group differences in magnitudes (Figure 1b), spatial configurations (Figure 2j) and/or directions (Figure 3) of future distributional changes are realistic. Incongruent reactivity to environmental changes between biologically dependent taxa and groups (e.g. plants and bacteria; Yashiro et al., 2018) can lead to shuffled species interactions and disruptions in biotic networks (Alexander, Diez, & Levine, 2015; Araújo & Luoto, 2007; Walther, 2010). Furthermore, land use changes, such as expansion of forest cover (Gago‐Silva et al., 2017), influence both directly and indirectly (though altering the abiotic conditions) the re‐distribution of taxa (Guo, Lenoir, et al., 2018). Thus, in future studies, land use changes should be accounted for together with climate change.

While our analyses revealed important ecological differences in the expected climate change responses among the studied taxonomic groups, three commonalities also emerged. First, among the six topoclimatic variables used, annual mean temperature had the highest average relative contribution in predicting taxa distribution within our study area. While our study area covers a variety of climatic conditions due to the c. 3,000 m elevation range (equivalent to approx. >1,500 km latitudinal gradient; Montgomery, 2006), the generality of this, and our other findings, across ecoregions (Miller et al., 2018; Petitpierre et al., 2016), geographical scales (Birkhofer et al., 2012; Graf, Bollmann, Suter, & Bugmann, 2005) and other groups (e.g. below‐ground macrofauna; Decaëns, 2010) remains to be assessed. Second, the taxa with the strongest changes in suitable habitat generally have lower prevalence and narrower environmental niches. However, these changes were not systematically gains or losses of habitat, making these rare specialists sensitive yet not necessarily vulnerable to climate change (Foden et al., 2013; Moritz & Agudo, 2013). Third, we observed the strongest predicted changes in taxa distributions at mid‐ and high elevations. This pattern has been commonly observed for plants under climate change forecasts in mountain environments (Dirnböck, Dullinger, & Grabherr, 2003; Engler et al., 2011; Steinbauer et al., 2018), and here we show that it applies to other taxonomic groups as well. Nevertheless, this interpretation must be taken with caution because boundary effects are expected to affect the predictions at low elevations. Boundary effects may arise when omitting taxa occurrences from warmer conditions outside the study area (i.e. resulting in erroneous predictions based on truncated response curves; Hannemann, Willis, & Macias‐Fauria, 2016) and restricting the taxa to the ones actually found in the study area (i.e. omitting the taxa migrating and invading from lower altitudes as a result of climate tracking; Lenoir, Gégout, Marquet, de Ruffray, & Brisse, 2008; Menéndez, González‐Megías, Jay‐Robert, & Marquéz‐Ferrando, 2014; Petitpierre et al., 2016).

Our study also demonstrated the importance of two methodological developments, which have so far rarely been considered in climate‐biodiversity impact studies involving a large number of taxa. Notably, the ESM approach allows more taxa to be included, particularly less frequent and rare taxa, and thus allows more complete modelling of local communities (Breiner et al., 2015, 2018; Lomba et al., 2010), and tends to outperform standard SDMs even for more common species with large sample sizes (Breiner et al., 2015). Interestingly, the predicted changes in PPO were generally lower than those predicted by standard SDMs, a finding never reported before. The likely cause for this is that the bivariate models in ESMs are less prone to multicollinearity that can bias the estimated effects of variables within standard SDMs, where all predictors are included simultaneously (Graham, 2003). Second, we avoided here the potential bias and loss of information resulting from the traditionally performed thresholding of predictions, because we applied metrics that do not require binarizing probabilities into presences–absences for evaluating the models and assessing the future changes (Jiménez‐Valverde & Lobo, 2007; Lawson et al., 2014; Nenzen & Araujo, 2011; Scherrer, Mod, & Guisan, 2020; Vaughan & Ormerod, 2005).

Finally, we would like to point out four insights related to multigroup biodiversity assessments. First, global data sets would be ideal to test the general effects of climate on multiple taxa and groups. However, existing global data sets are strongly biased in terms of taxonomy (Lenoir et al., 2020; Lenoir & Svenning, 2015; Troudet, Grandcolas, Blin, Vignes‐Lebbe, & Legendre, 2017) and spatial coverage (Boakes et al., 2010) affecting SDMs and hindering inter‐taxa/‐group comparisons (Beck, Böller, Erhardt, & Schwanghart, 2014). Here we reduced considerably these potential biases by sampling many taxa through a similar random stratified sampling design (Hirzel & Guisan, 2002) in a same region, where variation in elevation across the study area creates wide environmental gradients. Second, we acknowledge that our models used input data that were taxonomically attributed differently among the studied groups, as species‐level information is currently not available and difficult to determine for soil microorganisms unlike for above‐ground groups. One notable reason for this is that macroorganisms can be taxonomically catalogued phenotypically, whereas microorganisms were determined from DNA sequencing data only. As a proxy for microorganisms at the species level, we could have used the frequently deployed OTUs, but their total number (>38,000) exceeded the available storage space and computational capacity regarding the post‐processing of data (e.g. calculating Euclidian distances in n‐dimensional space). Hence, soil microorganisms were aggregated to the putative level of genera. While this might have added noise in the data for model calibration, we tested this by artificially aggregating the plant species to genera and the microbial genera to orders. This did not systematically affect model performance, and the differences in model performance between taxonomic resolutions (on average 0.06) remained minor compared to the differences in performance between above‐ and below‐ground taxa (on average 0.13). Third, we focused on identifiable groups to present an ecological narrative for key aspects of mountain ecosystems (e.g. amphibians/reptiles—predation, grasshoppers—herbivory, bumblebees—pollination, plants—primary production, fungi—decomposition, bacteria/protists—nutrient cycling), but we acknowledge that interpretations of climate influence might have been different if the taxa had been organized by another grouping (Dormann & Woodin, 2002; Harrison et al., 2010; Peay, Kennedy, & Talbot, 2016; Soudzilovskaia et al., 2013). Fourth, despite using state‐of‐the‐art methodologies, there might be biases in data, models and projections related to taxa detection (Benoit, Jackson, & Ridgway, 2018), sampling methods (e.g. point observations of amphibians and reptiles vs. presence–absence observations of other groups; Aarts, Fieberg, & Matthiopoulos, 2012; field observations vs. eDNA; Porter & Hajibabaei, 2018; Yoccoz, 2012; including observations also from forested sites for amphibians and reptiles while restricting the analyses to grasslands for other groups) and/or amplification of eDNA (Fernandes, Scherrer, & Guisan, 2019; Pinol, Mir, Gomez‐Polo, & Agusti, 2015). However, these issues tend to affect all intergroup biodiversity studies and should not alter specifically our findings compared to all others.

The need to understand the role of climate on biodiversity is particularly invigorated by the current threat of climate change on ecosystems (IPBES, 2019), and realistic future forecasts require incorporating multiple taxonomic and functional groups. Our models indicate both commonalities and differences in climate control among the studied groups. For all taxonomic groups, mean annual temperature was the most important topoclimatic predictor, and the strongest changes were predicted for rare taxa and at high altitudes. However, our results clearly show that the control of topoclimate is stronger on the distribution of above‐ground macroorganisms than of below‐ground microorganisms, for which soil microclimate likely better expresses temperature and moisture‐related constraints. This finding warns against potential disruptions in biotic networks but also implies that we still need to delineate the most ecologically meaningful environmental variables to comprehensively describe the distributions of soil microbiota and how the latter will be affected by climate and soil changes. Due to the important and intricate biological relationships among the different taxonomic groups, realistic forecasts for all members of the ecosystems will be crucial for the optimal assessment of future patterns of biodiversity. For this, future scenarios of other variables than macroclimatic ones, such as soil microclimatic and geochemical conditions, are also needed.

Appendix S1‐S3

Click here for additional data file.