# Multi-trophic β-diversity mediates the effect of environmental gradients on the turnover of multiple ecosystem functions

Camille Martinez-almoyna, Wilfried Thuiller, Loïc Chalmandrier, Marc

Ohlmann, Arnaud Foulquier, Jean-Christophe Clement, Lucie Zinger, Tamara

Münkemüller

► To cite this version:

Camille Martinez-almoyna, Wilfried Thuiller, Loïc Chalmandrier, Marc Ohlmann, Arnaud Foulquier, et al. Multi-trophic β-diversity mediates the effect of environmental gradients on the turnover of multiple ecosystem functions. Functional Ecology, 2019, 33 (10), pp.2053-2064. 10.1111/1365-2435.13393. hal-02276279

HAL Id: hal-02276279

https://hal.science/hal-02276279

Submitted on 20 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Multi-trophic β-diversity mediates the effect of environmental

2 gradients on the turnover of multiple ecosystem functions

3

4 Camille Martinez-Almoyna¹*, Wilfried Thuiller¹, Loïc Chalmandrier², Marc Ohlmann¹, Arnaud

5 Foulquier¹, Jean-Christophe Clément³, Lucie Zinger⁴ and Tamara Münkemüller¹,

6

7 ¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Laboratoire d'Ecologie Alpine, F-38000,

8 Grenoble, France

9 ²Department of Botany, University of Wyoming, Laramie, Wyoming, USA

10 ³CARRTEL, UMR 0042 Univ. Savoie Mont Blanc – INRA, FR-73376 Le Bourget-du-Lac, France

11 ⁴Ecole Normale Superieure, PSL Research University, CNRS, Inserm, Institut de Biologie de l'Ecole

12 Normale Superieure (IBENS), F-75005, Paris France

13 *Correspondence: Email: cmartinezalm@gmail.com

14

15 Keywords: biodiversity-ecosystem functioning, path analyses, multi-trophic, beta-diversity, Alpine

16 systems, elevation gradient.

17
18 Abstract

1. Much effort has been devoted to better understanding the effects of environment and biodiversity on ecosystem functioning. However, few studies have moved beyond measuring biodiversity as species richness of a single group and/or focusing on a single ecosystem function. While there is a growing recognition that along environmental gradients, the compositional turnover of multiple trophic groups influences not only productivity but multiple ecosystem functions, we do not know yet which components of multi-trophic β-diversity influence which ecosystem functions.

2. Here, we captured the biodiversity found in soils using environmental DNA to study total soil multi-trophic β-diversity (between all taxa regardless of their trophic group association), horizontal β-diversities (β-diversities within trophic groups) and vertical β-diversity (β-diversity across trophic groups) along a 1000-m elevational gradient in the French Alps. Using path analyses, we quantified how these β-diversity components mediate the effects of environmental turnover on the turnover of multiple ecosystem functions (i.e. productivity, N-cycling, N-leaching) and overall multifunctionality.

3. While we found a strong direct effect of soil properties on the turnover of multiple ecosystem functions, we also found an indirect effect of climate and soil properties through multi-trophic β-diversity. More specifically, only total multi-trophic β-diversity and the horizontal β-diversity of saprophytic fungi were strongly related to the turnover of multifunctionality and, to a lower extent, the turnover of productivity and N-cycling. Our results suggest that decomposition processes and resulting nutrient availability are key to understand how ecosystem functions change along edaphic and climatic gradients in alpine ecosystems.

4. Beyond alpine systems, our study stresses the paramount importance of considering the mediating role of multi-trophic diversity in environmental change impacts on different ecosystem functions and their multifunctionality.

41
42

# INTRODUCTION

43

44

45

46

47

48

49

50

51

52

53

The idea that biodiversity enhances ecosystem functioning and acts as a buffer against environmental changes is now a consensus in ecology (Cardinale et al., 2012; D. U. Hooper et al., 2005). Indeed, diverse communities exploit a wide range of available resources (complementarity effect) and often host the best performing species (selection effect, Hines et al., 2015; Loreau et al., 2001), which ultimately leads to enhanced ecosystem functioning (Bradford et al., 2014; Hector & Bagchi, 2007). This idea mostly comes from studies focusing on diversity within single trophic levels, single ecosystem functions (mainly productivity), and often carried out under controlled conditions (Cardinale et al., 2012; Loreau, 2010). However, understanding how environmental changes across space and time will affect natural ecosystems requires to consider the multiple trophic levels in natural communities and the various ecosystem functions carried out by multi-trophic diversity (Brose & Hillebrand, 2016; Mayor et al., 2017; Soliveres et al., 2016).

54

55

56

57

58

59

60

61

62

63

64

65

Environmental conditions such as climate, available carbon or nitrogen vary by nature over space and time, for example along elevational gradients (Körner, 2003). They directly influence ecosystem functions such as litter decomposition, N-cycling, and primary productivity by regulating the rates of bio-chemical processes (Sveinbjörnsson, Abadie, & Butler, 1995). These environmental conditions also select specific plant and soil assemblages and there is increasing recognition that the taxonomic and functional composition of these species assemblages will then influence the ecosystem functions as well (Hautier et al., 2018; Mori, Isbell, & Seidl, 2018; van der Plas et al., 2016). It has been recently highlighted that a better process-level understanding of the roles and functions of biodiversity may be gained through the study of the turnover of ecological communities, i.e. their β-diversity, a so far largely underexplored facet of biodiversity (Mori et al., 2018). Therefore, here we aim at deciphering direct and indirect effects of the environment by jointly considering the turnover of environmental conditions, ecosystem functions and β-diversity (van der Plas et al., 2016).

66

67

Moreover, recent studies have highlighted the importance of accounting for multi-trophic communities to fully understand the role of β-diversity in ecosystem functions (Hines et al., 2015;
Soliveres et al., 2016). Even if experimental work has identified effects at single trophic levels and for single ecosystem functions (Hines et al., 2015), these effects are likely to be more complex in natural systems. Indeed, the effects of different groups of species (same trophic level and/or similar functional characteristics, called trophic groups hereafter) can be simply additive, show synergies or go in different directions (Hines et al., 2015). For example, symbiotic fungi and herbivorous nematodes can have opposite effects on the biomass production of individual plant species (Brussaard, 1998). Therefore, teasing apart the effects of plant, symbiotic fungi and herbivorous nematodes diversity on ecosystem functioning requires a multi-trophic approach. Recently, Soliveres et al. (2016), showed that the richness variation of three trophic groups needs to be known to properly explain a single ecosystem function. In addition, depending on the amount of functional redundancy between trophic groups, their relative importance may change along environmental gradients (Setälä, Berg, & Jones, 2005). In other words, expanding the traditional focus on horizontal β-diversity (i.e. within a single trophic group) to a vertical perspective of multi-trophic β-diversity (i.e. variation across trophic groups) and a total multi-trophic β-diversity (i.e. of all taxa regardless of their trophic group) should critically improve our understanding of the role of multi-trophic β-diversity on ecosystem functioning (Figure 1).

Finally, natural ecosystems are characterized by many interconnected ecosystem functions. It is likely that distinct ecosystem functions respond differently to various components of multi-trophic β-diversity. A changing environmental context may also alter the links between ecosystem functions and multi-trophic β-diversity (Bradford et al., 2014; Setälä et al., 2005). The influence of horizontal and vertical β-diversities on ecosystem functioning could thus be mis-quantified when focusing on a single function, e.g. primary productivity, as a proxy for overall ecosystem functioning (Hector & Bagchi, 2007). Indeed, whereas turnover of productivity may be driven by the horizontal β-diversities of plants and fungi, turnover of nutrients may depend more on the horizontal β-diversity of bacteria (Moore et al., 2004; Wall, Bardgett, Behan-Pelletier, Herrick, Jones, Ritz, ... & van der Putten, 2012). The turnover of N-leaching, usually more related to abiotic characteristics of the ecosystem, might
instead be independent of multi-trophic β-diversity. Thus, studying horizontal and vertical β-diversity effects on different ecosystem functions should clarify how multifunctionality along environmental gradients is driven by different trophic groups (Delgado-Baquerizo et al., 2016).

Here, we study the direct and indirect effects of environmental turnover on the turnover of multiple ecosystem functions in mountain grasslands by explicitly considering the role of soil horizontal and vertical β-diversities (Figure 1). Along a 1000 m elevation gradient, we surveyed climate and soil properties, multi-trophic biodiversity of the most important trophic groups of the plant-soil compartment (as inferred from environmental DNA, Kress, García-Robledo, Uriarte, & Erickson, 2015; Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012), and several *in-situ* ecosystem processes (i.e. primary productivity, N-leaching, N-cycling, and multifunctionality). We first tested how total multi-trophic β-diversity mediates the effects of environmental turnover on turnover of multiple ecosystem functions and whether these effects varied between ecosystem functions. Second, we tested the independent effects of vertical multi-trophic β-diversity and horizontal β-diversities on the spatial turnover of multiple ecosystem functions. Finally, we jointly analyzed the relative contributions of horizontal β-diversities of different trophic groups to the spatial turnover of multiple ecosystem functions.

110
111 MATERIEL AND METHODS

112 Study site

113 The study was carried out in the French Alps, along a continuous elevation gradient (1750-
114 2725m a.s.l., Arves Massif, 45.12°N, 6,40°E) in a cow-grazed pasture. Subalpine grasslands dominate
115 the bottom of the gradient while sparsely vegetated alpine meadows characterize higher elevations.
116 Ten sites were sampled in summer 2012 along the same south-facing slope and they were separated
117 by an 100m elevation difference and an average distance of 340m. Each site consisted of two
118 10mx10m plots with homogeneous vegetation (see Chalmandrier et al., 2017, for more details).

119

120 Environmental variables

121 We considered four bioclimatic variables known to be strong drivers of community assembly in
122 mountain systems (Körner, 2003): soil temperature, solar radiation, growing season length (GSL) and
123 number of frost days. Growing season length and annual number of frost days were derived from
124 Landsat 7 and 8 imageries (Carlson, Choler, Renaud, Dedieu, & Thuiller, 2015). Solar radiations were
125 calculated with the area solar radiation tool in ArcGIS (version 10.2, 2013; Redlands, CA, USA) using a
126 2-m LiDAR digital elevation model and the clear sky model set to a sky size of 2800 pixels (Carlson et
127 al., 2015). Soil temperature was measured in the field.

128 To capture small scale variation, we added local topographic variation using elevation, slope, the
129 topographic wetness index (TWI) and the topographic position index (TPI) with a 50cm resolution
130 Digital Elevation Model derived from the airborne LiDAR data acquired the year of sampling.

131 Finally, we measured soil properties from three soil cores (10 cm depth, 5 cm diameter) taken in
132 a 1.25m x 1.25m subplot at the center of each plot. Two soil cores were weighed, 5mm-sieved and
133 pooled for quantifying soil moisture, soil organic matter content, pH, soil texture using standard
134 protocols (Robertson, Coleman, Sollins, & Bledsoe, 1999), and total soil N and C content using an
135 elemental analyzer (FlashEA1112, ThermoFisher Scientific, MA, USA). We added 100 mL of distilled
136 water to saturate the third soil core and calculate the water-filled pore space (WFPS). The soil core
137 was then air-dried, sieved and grounded to calculate gravel weight, total porosity and apparent
density (Legay et al., 2014).
138

139
140 We found that the first two axes of a principal component analyses run for all produced (and
141 normalized) environmental variables (figure S1) captured respectively 48.1% and 19.4% of the
142 environmental variability. We then split the environmental variables in two groups corresponding to
143 the major environmental gradients of the study area: a climatic gradient and a gradient of soil
144 properties.
145

Climate
---

146 To measure the turnover of climate (and topography), we focused on the set of environmental
147 variables that load strongest on the first axis of the PCA: solar radiation, soil temperature, mean and
148 variance of GSL, number of frost days, percentage of WFPS, TPI, TWI, altitude and slope. We then
149 computed Euclidean distances between each pair of plots for these environmental variables, and
150 defined these pairwise distances as a measure of *climatic turnover*.

Soil properties
---

152 To measure the turnover of soil properties, we focused on the set of environmental variables
153 that load strongest on the second axis of the PCA: gravel mass, apparent density, soil porosity, soil
154 pH, percentage of organic matter and of total N. We then computed Euclidean distances between
155 each pair of plots for these environmental variables, and defined these pairwise distances as a
156 measure of *soil properties turnover*.

157
158
**Multi-trophic data**
159

Sampling
---

160 To get an estimate of the diversity found in each of the sampled plot, we analyzed soils samples
161 through DNA metabarcoding. In each plot, we collected 21 soil samples (10 cm depth, 5 cm
diameter) to account for small-scale heterogeneity (ten samples on each diagonal and a
supplementary central one). Soil biodiversity was measured with four different DNA markers. Two
163
universal markers, one amplifying all Eukaryotes (v6-v7 region of the 18S rRNA gene) and one
164
amplifying all Bacteria (v5-v6 region of the 16S rRNA gene), were used to obtain a general overview
165
of the soil multi-trophic composition of the plots. We also used two additional markers focused on
166
fungi (Internal Transcribed Spacer 1) and vascular plants (Chloroplast trnL-P6 loop). Molecular
167
analyses and data curation are presented in Ohlmann et al. (2018). We pooled the 21 samples per
168
plot to obtain a single community of molecular taxonomic units (MOTUs) per plot and converted the
169
data into presence-absences. The curated sequencing data as well as associated metadata are
170
available on the Dryad Digital Repository under accession https://doi.org/10.5061/dryad.5b58400.
171

172

173 Definition of the trophic groups

Following Ohlmann et al. (2018), we grouped the MOTUs based on their taxonomic affiliation,
174
shared trophic resources and main functions, such as organic matter decomposition. Hereafter, we
175
will refer to them as *trophic groups*, while other definitions like trophic guilds or tropho-functional
176
groups could have been given.
177

We first assigned all vascular plant MOTUs to a single trophic group of primary producers, which
178
directly contribute to ecosystem productivity and to the quality of litter inputs to soils with potential
179
impacts on nitrogen cycling and loss (Hooper & Vitousek, 1998). Second, we considered microbial
180
MOTUs (i.e. bacteria and fungi) as decomposers. We distinguished bacteria from fungi since their
181
metabolic functions do not fully overlap (De Boer, Folman, Summerbell, & Boddy, 2005). We further
182
used the FUNguid database (Nguyen et al., 2016) distinguish symbiotic (i.e. mycorrhizal), saprophytic
183
and pathogenic fungi. Finally, the diversity of soil micro- and mesofauna regulates the activity and
184
abundance of soil microbes and associated nutrient recycling by shredding/grazing leaf litter or
185
feeding directly on microbes (Crowther et al., 2015; Pulleman et al., 2012). They should thus shape
186
the multi-trophic β-diversity and have a determinant impact on multiple ecosystem functions. We
187
thus extracted from the Eukaryotes-marker, all MOTUs assigned to nematodes, springtails, and
188
189 oribatid mites which are the most abundant belowground multicellular organisms (Wall et al. 2012).
190 We further distinguished bacterivore from herbivore/fungivore nematodes using the NEMguard
191 database (Nguyen et al., 2016), where nematode species trophic modes are defined on the basis of
192 the morphology of their mouth.
193

194 Different components of multi-trophic β-diversity

195 First, we computed the true turnover (Baselga, 2010) of the whole soil community (i.e. ignoring
196 trophic groups) as a measure of the *total multi-trophic β-diversity* using the Simpson dissimilarity
197 index between each pair of plots (Eq. 1).
198

$$
\beta_{kl} = \frac{\min(\|P_{k.}\| - \|P_{k.P_l.}\|, \|P_l.\| - \|P_{k.P_l.}\|)}{\|P_{k.P_l.}\| + \min(\|P_{k.}\| - \|P_{k.P_l.}\|, \|P_l.\| - \|P_{k.P_l.}\|)} \quad (1)
$$

199 Where *k* and *l* are the focal plots (such as 1 ≤ *k* ≤ *K* ∧ 1 ≤ *l* ≤ *K*, *K* the total number of plots),
200 *P* = (*P*<sub>*ki*</sub>)<sub>1≤*k*≤*K*</sub> is the community matrix of the whole soil community (*n* is the total number of
201 MOTUs and *P*<sub>*ki*</sub> ∈ {0,1} denotes the presence of the MOTU *i* in plot *k*). *P*<sub>*k*</sub> is the community vector of
202 the plot *k* and ||*P*<sub>*k*</sub>|| = Σ<sub>*i*=1</sub><sup>*n*</sup> *P*<sub>*ki*</sub> denotes the MOTUs richness of the plot *k*. ||*P*<sub>*k*</sub>. *P*<sub>*l*</sub>|| is the number of
203 shared MOTUs between plot *k* and *l*.
204

204 Second, we calculated the true turnover (Baselga, 2010) of each single trophic group as
205 measures of their horizontal β-diversities using the Simpson dissimilarity index between each pair of
206 plots (Eq. 2).
207

$$
\beta_{kl}^{g} = \frac{\min(\|P_{k.}^{g}\| - \|P_{k.P_{l.}^{g}}\|, \|P_{l.}^{g}\| - \|P_{k.P_{l.}^{g}}\|)}{\|P_{k.P_{l.}^{g}}\| + \min(\|P_{k.}^{g}\| - \|P_{k.P_{l.}^{g}}\|, \|P_{l.}^{g}\| - \|P_{k.P_{l.}^{g}}\|)} \quad (2)
$$

208 Where *k* and *l* are the focal plots (such as 1 ≤ *k* ≤ *K* ∧ 1 ≤ *l* ≤ *K*, *K* the total number of plots), *g* is
209 the focal trophic group, *P*<sup>*g*</sup> = (*P*<sub>*ki*</sub><sup>*g*</sup>)<sub>1≤*k*≤*K*</sub> is the community matrix of the group *g* (*n*<sub>*g*</sub> is the total
210 number of MOTUs and *P*<sub>*ki*</sub><sup>*g*</sup> ∈ {0,1} denotes the presence of the MOTU *i* in plot *k*). *P*<sub>*k*</sub><sup>*g*</sup> is the
211 community vector (of the group *g*) of the plot *k* and ||*P*<sub>*k*</sub><sup>*g*</sup>|| = Σ<sub>*i*=1</sub><sup>*n*</sup> *P*<sub>*ki*</sub><sup>*g*</sup> denotes the richness of the
212 plot *k* for the group *g*.
213
213

214

215

216

And finally, we calculated the richness variation across trophic groups as a measure of the vertical b-diversity using the Bray-Curtis (BC) dissimilarity index between each pair of plots (Legendre & Legendre, 1998):

217

$$BC\beta_{kl} = \frac{2\sum_{g=1}^{G} \min(\|P_{k,g}\|, \|P_{l,g}\|)}{\sum_{g=1}^{G} \|P_{k,g}\| + \|P_{l,g}\|} \quad (3)$$

218

Where G is the total number of trophic groups.

219

220

## Multifunctionality data

221

### Productivity

222

223

224

225

226

227

228

229

230

231

232

233

234

### N-Cycling

235

236

237

We measured turnover of primary productivity using two types of data. First, we measured the total green biomass of the year in situ, which constitutes a good proxy of the annual biomass increment. We harvested all the vegetation at the ground level on a 1mx1m plot, sorted out the current year live and recently senescent material, dried and weighed it. We used the spatial mean and variability of the Normalized Difference Vegetation Index (NDVI), a proxy of the photosynthetic activity (Yoder & Waring, 1994), as a complementary measure. For each plot, we computed spatial mean and variance of NDVI using an air-borne 15cm resolution infrared picture taken at the productivity peak in August 2012 (Yoder & Waring, 1994) for the whole elevation gradient. We then calculated their Euclidean distance for each pair of plots. Productivity turnover was then defined as the Euclidean distances between each pair of plots for the three measures: total green biomass, spatial mean and variance of NDVI.

To estimate turnover of N-cycling, we used the same two soil cores as for soil physico-chemical measurements (see environmental data). We estimated the decomposability of organic matter by calculating the C:N ratio. Soil nutrients (nitrate ($NO_3^-$), ammonium ($NH_4^+$), total dissolved nitrogen
238 (TDN) and dissolved organic nitrogen (DON) were determined from soil extracts with 0.5 M K₂SO₄
239 (Jones & Willett, 2006) using a photometric analyzer (Gallery Plus, Thermo Fisher Scientific, Waltham,
240 Massachusetts, USA). We estimated the relative availability of different forms of N by calculating the
241 ratio of total dissolved nitrogen on total nitrogen (N-TDN/N) and the ratio of dissolved organic
242 nitrogen on mineral nitrogen (N-DON/(N-NO₃+N-NH₄)) because dissolved organic N, particularly
243 amino acids, constitutes a large portion of the N budget of plants in these high latitude ecosystems
244 (Loomis, Ruess, Sveinbjörnsson, & Kielland, 2006). Finally, potential nitrogen mineralization rates
245 were used to estimate the ability of the ecosystem to mobilize soil organic matter and supply mineral
246 N (Legay et al., 2016), and measured based on anaerobic incubations of fresh soil subsamples (dark, 7
247 days, 40 °C), during which organic N was mineralized and accumulated as NH₄⁺ (Waring & Bremner,
248 1964; Wienhold, 2007). The difference between NH₄⁺ contents in a given sample before (t1) and after
249 the anaerobic incubation (t2) gave PNM = [(NH₄⁺-N)t2(NH₄⁺-N)t1]/dw/7 days. *N-cycling turnover* was
250 then defined as the Euclidean distance between each pair of plots for C/N, N-TDN/N, potential
251 mineralization and N-DON/(N-NO₃+N-NH₄).

252

253 N-leaching

254 As a proxy for nitrogen losses through leaching we measured potential ammonium and nitrate
255 leaching from the percolate. Percolated water was filtered through a Whatmann filter paper n°42
256 (2.5µm pore size) and analyzed using a photometric analyzer to calculate potentially leached NH₄⁺-N
257 and NO₃⁻-N (Legay et al., 2016). *N-leaching turnover* was then defined as the Euclidean distance
258 between each pair of plots for potential ammonium and nitrate leaching.

259

260 Multifunctionality

261 We calculated *multifunctionality turnover* as the Euclidean distance between each pair of plots
262 for all measured ecosystem processes. Correlations of each ecosystem function to the
multifunctionality turnover are presented in Appendix (Figure S2). We also run all the analyses for each ecosystem process separately.

265

## 266 Statistical Analyses

To address our first two objectives, we performed a set of path analyses using the lavaan package in R (Figure 2, R Core Team, 2018; Rosseel, 2012). All variables were first centered and normalized. First, we built four path analyses independently for each of the three single ecosystem functions and multifunctionality to determine the effects of climatic and soil properties turnover and of total multi-trophic β-diversity. Second, we built four path analyses independently for each of the three single ecosystem functions and multifunctionality to determine the effects of environmental turnover (turnover of climate and soil properties) and of vertical β-diversity. Finally, we built 36 independent path analyses independently to assess the effect of each of the 9 trophic groups on each of the three single ecosystem functions and multifunctionality while including in each model the effects of the two sets of variables (climate and soil properties) for environmental turnover.

To address our final objective, we built an integrated path model for each ecosystem function and multifunctionality. We selected variables in a semi-explorative way. We always included climatic and soil properties turnover to account for their direct effects on both biodiversity and ecosystem functions. In respect to multi-trophic β-diversity, we only included the trophic groups that had a significant horizontal effect and estimated their direct effects on turnover of multiple ecosystem functions. Since we did not have any prior expectation on their causal relationships, we did not specify any directional link between them. However, we let free covariances between each pair of trophic groups to account for the effect of their probable interactions. In turn, all our path analyses were saturated and we thus could not test the structure of the models.

For all above-mentioned path analyses, we tested for the significance of the effects using a parametric bootstrap procedure based on 10,000 sampling with simultaneous replacement of rows and columns to quantify confidence intervals for the parameters (Fournet et al., 2018). This
289 approach is suitable when modelling pairwise distances where the elements of one row/column are

290 not independent of each other (Fourtune et al., 2018). Parameters are significant when the

291 confidence intervals do not include zero.

292
293 **RESULTS**

294
Turnover in soil properties was a consistent and strong direct predictor of multifunctionality and
295
of the three ecosystem functions, while climatic turnover was not (Figure 3 & 4). Total multi-trophic
296
β-diversity was a significant predictor of productivity and multifunctionality turnover and almost as
297
important as soil properties turnover. However, it did not explain N-cycling and N-leaching turnover
298
(Figure 3).

299
In contrast, when using vertical β-diversity instead of total multi-trophic β-diversity, the diversity
300
effect vanished and only soil properties turnover had a significant effect on turnover of ecosystem
301
functions, including multifunctionality (Figure 4). In other words, the effect of total multi-trophic β-
302
diversity on turnover of multiple ecosystem functions was not driven by vertical β-diversity.

303
Finally, focusing on the independent effects of horizontal β-diversities of the different trophic
304
groups, the most important predictors of multifunctionality and productivity turnover were soil
305
properties turnover and horizontal β-diversities of saprophytic fungi, followed by symbiotic fungi,
306
herbi-fungivore nematodes, pathogenic fungi and bacteria (Figure 5, S3). For N-cycling turnover, the
307
most important predictors were soil properties turnover and horizontal β-diversities of saprophytic
308
fungi, followed by symbiotic fungi and pathogenic fungi (Figure 5, S3). No component of multi-
309
trophic β-diversity was linked to N-leaching turnover (Figure 5), for which soil properties turnover
310
was the only significant predictor (Figure S7). For the non-aggregated ecosystem processes, different
311
subsets of trophic groups were the best predictors, but the general trends remained (Figure S7).

312

313
The integrated path models that considered several trophic groups together, but were built
314
independently for each ecosystem function and multifunctionality, suggested how the pieces of the
315
puzzle can be put together (Figure 6). Consistent with the results above, soil properties were of
316
primary importance for all considered ecosystem functions and multifunctionality directly, whereas
317
climate was not (Figure 6, Figure S4 & S5). Similar figures emerged when considering each ecosystem
318
process independently (Figure S7).
319 Climatic turnover and soil properties turnover influenced horizontal β-diversities of most of the trophic groups (Figure S6). Pathogenic fungi (Figure 6, Figure S4 & S5, Table S1) and, to a lower extent, saprophytic fungi (Table S1) were directly affected by climate. Bacteria and, to a lesser extent, herbi-fungivorous nematodes and saprophytic fungi were directly affected by soil properties (Figure 6, Figure S4 & S5, Table S1).

324 Importantly, all trophic groups were interlinked and were thus involved in the indirect effects of climate and soil properties on multiple ecosystem functions. Pathogenic fungi were strongly linked to symbiotic and saprophytic fungi (Figure 6, Figure S4 & S5). Bacteria were linked especially to saprophytic fungi (Figure 6, Figure S4) and to a lesser extent to symbiotic fungi (Table S1).

328 Finally, the only significant biotic predictor for the turnover of multifunctionality was the horizontal β-diversity of saprophytic fungi (Figure 6). Additionally, horizontal β-diversity of herbi-fungivorous nematodes showed a strong effect even if not significant (Table S1). Results for productivity turnover showed comparable trends but with weaker effects (Figure S4). The combined horizontal β-diversities did not explain N-cycling turnover (Figure S5). N-leaching turnover was not considered here since it was not positively linked to any of the multi-trophic β-diversity components in the previous path analyses. When we considered the non-aggregated ecosystem processes in the same way, different trophic groups, also mostly fungi, had a significant effect on different ecosystem processes: saprotroph fungi on dissolved N:total N, symbiotic fungi on soil C:N, pathogenic on dissolved organic N:dissolved inorganic N and fungi-herbivorous nematodes on total green biomass (Figure S7.C, S7.D).

339

340
341

DISCUSSION

342

343

344

345

346

Brose & Hillebrand (2016) recently called for a shift from traditional biodiversity and ecosystem functioning experiments under controlled conditions to natural systems with (i) multi-trophic communities, (ii) varying environmental conditions, and (iii) larger spatial scales. Here, we tackled this challenge by using multi-trophic β-diversity measures to investigate how they mediate the influence of environmental turnover on the turnover of multiple ecosystem functions.

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

In our study along an elevation gradient in the French Alps, we found that environmental turnover influenced ecosystem functioning turnover both directly and indirectly via biodiversity turnover. Surprisingly, only the turnover of soil but not climatic properties contributed significantly to the direct pathway. Due to varying limitations in soil organic matter and water availability, decomposition rates and nutrient retention times changed across the gradient, with direct consequences on N-cycling, productivity and N-leaching (Gavazov, 2010; Louca et al., 2018). In contrast to our findings, a large body of literature on size reduction of alpine plants with increasing elevation led us expect that differences in temperature-related variables (e.g. Growing Season Length, number of frost days and solar radiations) captured in our climatic turnover variable (PCA axis 1, Figure S1) would drive at least productivity turnover (Choler, 2005; Körner, 2003; Schöb, Kammer, Choler, & Veit, 2009; Sundqvist, Sanders, & Wardle, 2013). However, a few recent studies that also included the indirect effect of climate through biodiversity support our results. Together this suggests that under natural conditions the direct effect of climate on ecosystem functioning can be overridden and/or compensated by complex cascading effects from microorganisms to plants. For example, Crowther et al. (2015) showed that turnover of ecosystem functions can be bottom-up regulated by the climatic conditions in harsh environments but top-down regulated through consumption of the microbes by the soil fauna, when environmental constraints are less extreme. Also, at the global scale, Delgado-Baquerizo et al. (2016) reported only a negligible effect of temperature on multifunctionality, while the effect of microbial community diversity was always strong and positive. Our results also demonstrated that on top of its direct response to soil
properties, turnover of multiple ecosystem functions depends on soil and climatic turnovers through
their effects on the β-diversities of several trophic groups (Figure S6).

In order to untangle these cascading effects, we contrasted the effects of several multi-trophic α-
and β-diversity aspects on multiple ecosystem functions. Unlike total β-diversity (Figure 3), total α-
diversity (Figure S8) appeared to be a weak predictor of ecosystem functioning, suggesting that a
number of species perform similar functions in soil communities, i.e. they are functionally redundant.
Also, vertical β-diversity did not cause turnover of ecosystem functions (Figure 4), indicating that
whatever their richness, trophic groups cluster particularly functionally redundant species (Louca et
al., 2018; Setälä et al., 2005). Nevertheless, intra-group functional redundancy level may depend on
the studied trophic group. For example, Mori et al. (2015) argue that changes in the composition of
fungal communities can modify ecosystem functioning by altering which ecosystem functions are
preferentially provided, while Louca et al. (2018) showed that bacterial taxa have similar broad
metabolic potential, leading most ecosystem functions to be barely dependent on the taxonomic
composition of bacterial communities. Likewise, in our study, the relationship between α-diversity
and ecosystem functions was linear and positive only for a small selection of trophic groups and
ecosystem functions (Figure S8). Together, these results suggest that a few of our trophic groups
contain functionally complementary species, pivotal for ecosystem functioning.

To avoid an overestimation of the direct effect of these key trophic groups on ecosystem functioning
(Ohlmann et al., 2018), the horizontal β-diversities of their interaction partners partaking in the
studied ecosystem functions should not be omitted (Brose & Hillebrand, 2016; Hines et al., 2015;
Thompson et al., 2012). As such, we found that most trophic groups influenced ecosystem functions
when considered alone (figure 5, figure S6). However, when the horizontal β-diversities of their
interaction partners were added, only the saprophytic fungi β-diversity was significantly linked to the
turnover of multifunctionality, and to a lesser extent to productivity and N-cycling (figure 6, figure S4,
figure S5, Table S1). Saprophytic fungi are critical in the decomposition process because they produce
exo-enzymes processing complex organic compounds (Baptist et al., 2008; Moore et al., 2004; Wall
et al. 2012). The resulting simpler compounds, scarce in alpine environments, become available for plant and microbial absorption (Moore et al., 2004; Wall et al., 2012). Changes in saprophytic fungi communities should thus alter nutrient availability, which impacts productivity, N-cycling and multifunctionality (Mori et al., 2016; Valencia et al., 2018). Nevertheless, also the other trophic groups β-diversities indirectly drove ecosystem functions turnover through their strong links to the saprophytic fungi β-diversity. In line with a growing number of studies, this suggests that species links in soil food webs could be pivotal for ecosystem functioning and that they need to be accounted for in addition to just the β-diversity of species within (horizontal) and across (vertical) trophic groups (Ohlmann et al., 2019; Setälä et al., 2005; Thompson et al., 2012). In particular, we found that specific regulators of fungi (e.g. herbi-fungivorous nematodes, figure 6, S4) and competitors (e.g. bacteria and other fungal trophic groups, De Boer et al., 2005; Nazir, Warmink, Boersma, & Van Elsas, 2010) tend to have an indirect influence on multiple ecosystem functions through their influence on the compositional turnover of saprophytic fungi. In turn, climatic variables and ecosystem functions are not linked by a single indirect pathway but by multiple indirect pathways with saprophytic fungi as a corner stone. This result is in line with several other studies (Mori et al., 2016; Setälä et al., 2005; Valencia et al., 2018).

Our results on both an integrated measure of multifunctionality and more specific measures of three important functions of alpine ecosystems (i.e. productivity, N-cycling, N-leaching, Sundqvist et al., 2013) are also consistent with several studies arguing that major biogeochemical pathways are regulated by a reduced number of energy pathways, broadly supported by some key trophic groups and sometimes strongly driven by some limiting environmental conditions (Louca et al., 2018). At the highest scale of aggregation of multiple ecosystem functions, saprophytic fungi being keystones of multifunctionality suggests that, due to the strong nitrogen limitations in alpine ecosystems (Loomis et al., 2006; Robson, Lavorel, Clement, & Roux, 2007; Sundqvist et al., 2013; Sveinbjörnsson et al., 1995), decomposition is acting as a bottle neck for other ecosystem functions (Setälä et al., 2005; Wardle, 2002). Indeed, ammonium and nitrate produced by decomposers can be leached out,
419 adsorbed on soil particles, denitrified, volatilized, or assimilated by plants and bacteria. The balance
420 between these processes mainly depends on abiotic soil conditions which control water and solute
421 transfers, and on the nitrogen demand of plant and microorganisms (Robson et al., 2007) especially
422 during the productivity peak when plant N-uptake is the highest (Legay et al., 2016). Therefore,
423 weaker effects of each trophic group horizontal β-diversity on the turnover of productivity and N-
424 cycling could be explained by a functional complementary effect between trophic groups of fungi and
425 their specific biotic regulators, and feedbacks between trophic groups and soil. Adding other, and so
426 far, a bit less studied pathways (e.g. P-cycling, Leff et al., 2015), may have uncovered additional direct
427 and indirect links of environmental turnover and functioning turnover via β-diversities.
428

428

429 To conclude, looking at multiple ecosystem functions showed an important mediating role of multi-
430 trophic β-diversity. Our finding particularly supports the paramount importance of decomposers,
431 especially saprophytic fungi, and of their interaction partners, in nitrogen limited alpine systems. In
432 this sense, our study supports recent calls to explicitly integrate known interactions between trophic
433 groups (Thompson et al., 2012) and for temporal monitoring of the ecosystems (Brose & Hillebrand,
434 2016), two missing steps towards a better understanding of functional complementarity and
435 selection effects of diversity and of their role in multiple ecosystem functions (Hines et al., 2015) and
436 their trade-offs in natural systems.
437

436
437 Author Contributions statement:

438 CMA, TM and WT conceptualized the overall study and methodology. LC, JCC, TM, WT and LZ
439 collected the data. LZ run the bioinformatic pipelines. MO built the trophic groups and CMA run all
440 statistical analyses. CMA led the interpretation of the results and the writing with contributions of all
441 authors.

442

443 Acknowledgements

444 We thank the large team of researchers and students that helped collect the data. The research
445 received funding from the French Agence Nationale de la Recherche (ANR) through the GlobNets
446 (ANR-16-CE02-0009) project, and from "Investissement d'Avenir" grants managed by the ANR
447 (Trajectories: ANR-15-IDEX-02; Montane: OSUG@2020: ANR-10-LAB-56).

448

449 Data accessibility statement:

450 Data will be provided after acceptance.

451

452 References:

453 Baptist, F., Zinger, L., Clement, J. C., Gallet, C., Guillemin, R., Martins, J. M. F., ... Geremia, R. (2008).

454 Tannin impacts on microbial diversity and the functioning of alpine soils: A multidisciplinary approach. *Environmental Microbiology*, 10(3), 799-809. doi:10.1111/j.1462-2920.2007.01504.x

456 Baselga, A. (2010). Partitioning the turnover and nestedness components of beta diversity. *Global Ecology and Biogeography*, 19(1), 134–143. doi:10.1111/j.1466-8238.2009.00490.x

458 Bradford, M. A., Wood, S. A., Bardgett, R. D., Black, H. I. J., Bonkowski, M., Eggers, T., ... Jones, T. H. (2014). Discontinuity in the responses of ecosystem processes and multifunctionality to altered soil community composition. *Proceedings of the National Academy of Sciences*, 111(40), 14478–14483. doi:10.1073/pnas.1413707111

462 Brose, U., & Hillebrand, H. (2016). Biodiversity and ecosystem functioning in dynamic landscapes.
463 *Philosophical Transactions of the Royal Society B: Biological Sciences*, 371(1694), 20150267.

464 doi:10.1098/rstb.2015.0267

465 Brussaard, L. (1998). Soil fauna, guilds, functional groups and ecosystem processes. *Applied Soil*

466 *Ecology*, 9(1–3), 123–135. doi:10.1016/S0929-1393(98)00066-3

467 Cardinale, B. J., Duffy, J. E., Gonzalez, A., Hooper, D. U., Perrings, C., Venail, P., ... Grace, J. B. (2012).

468 Biodiversity loss and its impact on humanity. *Nature*, 486(7401), 59–67.

469 doi:10.1038/nature11148

470 Carlson, B. Z., Choler, P., Renaud, J., Dedieu, J. P., & Thuiller, W. (2015). Modelling snow cover

471 duration improves predictions of functional and taxonomic diversity for alpine plant

472 communities. *Annals of Botany*, 116(6), 1023–1034. doi:10.1093/aob/mcv041

473 Chalmandrier, L., Münkemüller, T., Colace, M.-P., Renaud, J., Aubert, S., Carlson, B. Z., ... Thuiller, W.

474 (2017). Spatial scale and intraspecific trait variability mediate assembly rules in alpine

475 grasslands. *Journal of Ecology*, 105(1), 277–287. doi:10.1111/1365-2745.12658

476 Choler, P. (2005). Consistent Shifts in Alpine Plant Traits along a Mesotopographical Gradient. *Arctic*,

477 *Antarctic, and Alpine Research*, 37(4), 444–453. doi:10.1657/1523-

478 0430(2005)037[0444:CSIAPT]2.0.CO;2

479 Crowther, T. W., Thomas, S. M., Maynard, D. S., Baldrian, P., Covey, K., Frey, S. D., ... Bradford, M. A.

480 (2015). Biotic interactions mediate soil microbial feedbacks to climate change. *Proceedings of*

481 *the National Academy of Sciences*, 112(22), 7033–7038. doi:10.1073/pnas.1502956112

482 De Boer, W., Folman, L. B., Summerbell, R. C., & Boddy, L. (2005). Living in a fungal world: Impact of

483 fungi on soil bacterial niche development. *FEMS Microbiology Reviews*, 29(4), 795–811.

484 doi:10.1016/j.femsre.2004.11.005

485 Delgado-Baquerizo, M., Maestre, F. T., Reich, P. B., Jeffries, T. C., Gaitan, J. J., Encinar, D., ... Singh, B.

486 K. (2016). Microbial diversity drives multifunctionality in terrestrial ecosystems. *Nature*

487 *Communications*, 7, 10541. doi:10.1038/ncomms10541

488 Fourtune, L., Prunier, J. G., Paz-Vinas, I., Loot, G., Veyssière, C., & Blanchet, S. (2018). Inferring
489 Causalities in Landscape Genetics: An Extension of Wright's Causal Modeling to Distance

490 Matrices. *The American Naturalist*, **191**(4), 491–508. doi:10.1086/696233

491 Gavazov, K. S. (2010). Dynamics of alpine plant litter decomposition in a changing climate. *Plant and

492 Soil*, **337**(1–2), 19–32. doi:10.1007/s11104-010-0477-0

493 Hautier, Y., Isbell, F., Borer, E. T., Seabloom, E. W., Harpole, W. S., Lind, E. M., ... Hector, A. (2018).

494 Local loss and spatial homogenization of plant diversity reduce ecosystem multifunctionality.

495 *Nature Ecology & Evolution*, **2**(1), 50–56. doi:10.1038/s41559-017-0395-0

496 Hector, A., & Bagchi, R. (2007). Biodiversity and ecosystem multifunctionality. *Nature*, **448**(7150),

497 188–190. doi:10.1038/nature05947

498 Hines, J., van der Putten, W. H., De Deyn, G. B., Wagg, C., Voigt, W., Mulder, C., ... Eisenhauer, N.

499 (2015). Towards an Integration of Biodiversity–Ecosystem Functioning and Food Web Theory to

500 Evaluate Relationships between Multiple Ecosystem Services. In Guy Woodward & David A.

501 Bohan (Eds.), *Advances in Ecological Research* (Vol. 53, pp. 161–199).

502 doi:10.1016/bs.aecr.2015.09.001

503 Hooper, D. U., Chapin, F. S., Ewel, J. J., Hector, A., Inchausti, P., Lavorel, S., ... Wardle, D. A. (2005).

504 Effects of Biodiversity on Ecosystem Functioning: A Consensus of Current Knowledge. *Ecological*

505 *Monographs*, **75**(1), 3–35. doi:10.1890/04-0922

506 Hooper, David U., & Vitousek, P. M. (1998). Effects of Plant Composition and Diversity on Nutrient

507 Cycling. *Ecological Monographs*, **68**(1), 121. doi:10.2307/2657146

508 Jones, D. L., & Willett, V. B. (2006). Experimental evaluation of methods to quantify dissolved organic

509 nitrogen (DON) and dissolved organic carbon (DOC) in soil. *Soil Biology and Biochemistry*, **38**(5),

510 991–999. doi:10.1016/j.soilbio.2005.08.012

511 Körner, C. (2003). *Alpine Plant Life*. Berlin, Heidelberg: Springer Berlin Heidelberg. doi:10.1007/978-3-

512 642-18970-8

513 Kress, W. J., García-Robledo, C., Uriarte, M., & Erickson, D. L. (2015). DNA barcodes for ecology,

514 evolution, and conservation. *Trends in Ecology & Evolution*, **30**(1), 25–35.
doi:10.1016/j.tree.2014.10.008

Leff, J. W., Jones, S. E., Prober, S. M., Barberán, A., Borer, E. T., Firn, J. L., ... Fierer, N. (2015).

Consistent responses of soil microbial communities to elevated nutrient inputs in grasslands

across the globe. *Proceedings of the National Academy of Sciences*, 112(35), 10967–10972.

doi:10.1073/pnas.1508382112

Legay, N., Baxendale, C., Grigulis, K., Krainer, U., Kastl, E., Schloter, M., ... Lavorel, S. (2014).

Contribution of above- and below-ground plant traits to the structure and function of grassland

soil microbial communities. *Annals of Botany*, 114(5), 1011–1021. doi:10.1093/aob/mcu169

Legay, N., Lavorel, S., Baxendale, C., Krainer, U., Bahn, M., Binet, M.-N., ... Bardgett, R. D. (2016).

Influence of plant traits, soil microbial properties, and abiotic parameters on nitrogen turnover

of grassland ecosystems. *Ecosphere*, 7(11), e01448. doi:10.1002/ecs2.1448

Legendre, P., & Legendre, L. (1998). *Numerical ecology*. Elsevier.

Loomis, P. F., Ruess, R. W., Sveinbjörnsson, B., & Kielland, K. (2006). Nitrogen cycling at treeline:

Latitudinal and elevational patterns across a boreal landscape. *Ecoscience*, 13(4), 544–556.

doi:10.2980/1195-6860(2006)13[544:NCATLA]2.0.CO;2

Loreau, M. (2010). Linking biodiversity and ecosystems: towards a unifying ecological theory.

*Philosophical Transactions of the Royal Society B: Biological Sciences*, 365(1537), 49–60.

doi:10.1098/rstb.2009.0155

Loreau, M., Naeem, S., Raffaelli, D., Schmid, B., Bengtsson, J., Grime, J. P., ... Huston, M. A. (2001).

Biodiversity and Ecosystem Functioning: Current Knowledge and Future Challenges. *Science*,

294(5543), 804–808. doi:10.1126/science.1064088

Louca, S., Martiny, A. C., Polz, M. F., Doolittle, W. F., Nelson, M. B., Hahn, A. S., ... Hallam, S. J. (2018).

Function and functional redundancy in microbial communities. *Nature Ecology and Evolution*, 2,

936–943.

Mayor, J. R., Sanders, N. J., Classen, A. T., Bardgett, R. D., Clément, J.-C., Fajardo, A., ... Wardle, D. A.

(2017). Elevation alters ecosystem properties across temperate treelines globally. *Nature*,
541 542(7639), 91-95. doi:10.1038/nature21027
542 Moore, J. C., Berlow, E. L., Coleman, D. C., De Suiter, P. C., Dong, Q., Hastings, A., ... Wall, D. H. (2004).
543 Detritus, trophic dynamics and biodiversity. *Ecology Letters*, 7(7), 584-600. doi:10.1111/j.1461-0248.2004.00606.x
544
545 Mori, A. S., Isbell, F., Fujii, S., Makoto, K., Matsuoka, S., & Osono, T. (2016). Low multifunctional redundancy of soil fungal diversity at multiple scales. *Ecology Letters*, 19(3), 249-259.
546
547 doi:10.1111/ele.12560
548 Mori, A. S., Isbell, F., & Seidl, R. (2018). β-Diversity, Community Assembly, and Ecosystem Functioning. *Trends in Ecology & Evolution*, 33(7), 549-564. doi:10.1016/j.tree.2018.04.012
549
550 Nazir, R., Warmink, J. A., Boersma, H., & Van Elsas, J. D. (2010). Mechanisms that promote bacterial fitness in fungal-affected soil microhabitats. *FEMS Microbiology Ecology*, 71(2), 169-185.
551
552 doi:10.1111/j.1574-6941.2009.00807.x
553 Nguyen, N. H., Song, Z., Bates, S. T., Branco, S., Tedersoo, L., Menke, J., ... Kennedy, P. G. (2016).
554 FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild.
555 *Fungal Ecology*, 20, 241-248. doi:10.1016/j.funeco.2015.06.006
556 Ohlmann, M., Mazel, F., Chalmandrier, L., Bec, S., Coissac, E., Gielly, L., ... Thuiller, W. (2018).
557 Mapping the imprint of biotic interactions on β-diversity. *Ecology Letters*, 21(11), 1660-1669.
558 doi:10.1111/ele.13143
559 Ohlmann, M., Miele, V., Dray, S., Chalmandrier, L., O'Connor, L., & Thuiller, W. (2019). Diversity indices for ecological networks: a unifying framework using Hill numbers. *Ecology Letters*, 22(4), 737-747. doi:10.1111/ele.13221
560
561
562 Pulleman, M., Creamer, R., Hamer, U., Helder, J., Pelosi, C., Pérès, G., & Rutgers, M. (2012). Soil biodiversity, biological indicators and soil ecosystem services-an overview of European approaches. *Current Opinion in Environmental Sustainability*, 4(5), 529-538.
563
564
565 doi:10.1016/j.cosust.2012.10.009
566 R Core Team. (2018). R: A Language and Environment for Statistical Computing. Vienna, Austria.
567 Robertson, G. P., Coleman, D. C., Sollins, P., & Bledsoe, C. S. (1999). *Standard soil methods for long-term ecological research* (Vol. 2). Oxford University Press
568
569 Robson, T. M., Lavorel, S., Clement, J. C., & Roux, X. Le. (2007). Neglect of mowing and manuring leads to slower nitrogen cycling in subalpine grasslands. *Soil Biology and Biochemistry, 39*(4), 930-941. doi:10.1016/j.soilbio.2006.11.004
570
571
572 Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. *Journal of Statistical Software, 48*(2), 1-93. doi:10.18637/jss.v048.i02
573
574 Schöb, C., Kammer, P. M., Choler, P., & Veit, H. (2009). Small-scale plant species distribution in snowbeds and its sensitivity to climate change. *Plant Ecology, 200*(1), 91-104.
575
576 doi:10.1007/s11258-008-9435-9
577 Setälä, H., Berg, M. P., & Jones, T. H. (2005). Trophic structure and functional redundancy in soil communities. In R. Bardgett, M. Usher, & D. Hopkins (Eds.), *Biological Diversity and Function in Soils* (pp. 236-249). Cambridge: Cambridge University Press.
578
579
580 doi:10.1017/CBO9780511541926.014
581 Soliveres, S., van der Plas, F., Manning, P., Prati, D., Gossner, M. M., Renner, S. C., ... Allan, E. (2016). Biodiversity at multiple trophic levels is needed for ecosystem multifunctionality. *Nature*, 536(7617), 456-459. doi:10.1038/nature19092
582
583
584 Sundqvist, M. K., Sanders, N. J., & Wardle, D. A. (2013). Community and Ecosystem Responses to Elevational Gradients: Processes, Mechanisms, and Insights for Global Change. *Annual Review of Ecology, Evolution, and Systematics, 44*(1), 261-280. doi:10.1146/annurev-ecolsys-110512-135750
585
586
587
588 Sveinbjörnsson, B., Abadie, W., & Butler, A. (1995). Soil carbon and nitrogen mineralization at different elevations in the Chugach Mountains of South-Central Alaska, U.S.A. *Arctic and Alpine Research, 27*, 29-37. doi:10.2307/1552065
589
590
591 Taberlet, P., Coissac, E., Hajibabaei, M., & Rieseberg, L. H. (2012, April 1). Environmental DNA. *Molecular Ecology*. Wiley/Blackwell (10.1111/j.1365-294X.2012.05542.x)
592
593 Thompson, R. M., Brose, U., Dunne, J. A., Hall, R. O., Hladyz, S., Kitching, R. L., ... Tylianakis, J. M.

594 (2012). Food webs: Reconciling the structure and function of biodiversity. *Trends in Ecology and Evolution*. Elsevier Ltd. doi:10.1016/j.tree.2012.08.005

595

596 Valencia, E., Gross, N., Quero, J. L., Carmona, C. P., Ochoa, V., Gozalo, B., ... Maestre, F. T. (2018).

597 Cascading effects from plants to soil microorganisms explain how plant species richness and

598 simulated climate change affect soil multifunctionality, (June), 5642–5654.

599 doi:10.1111/gcb.14440

600 van der Plas, F., Manning, P., Soliveres, S., Allan, E., Scherer-Lorenzen, M., Verheyen, K., ... Fischer, M.

601 (2016). Biotic homogenization can decrease landscape-scale forest multifunctionality.

602 *Proceedings of the National Academy of Sciences*, 113(13), 3557–3562.

603 doi:10.1073/pnas.1517903113

604 Wall, D.H., Bardgett, R.D., Behan-Pelletier, V., Herrick, J.E., Jones, T.H., Ritz, K., Six, J., Strong, D.R. &

605 van der Putten, W. H. (2012). *Soil Ecology and Ecosystem Services*. (D. H. Wall, R. D. Bardgett, V.

606 Behan-Pelletier, J. E. Herrick, T. H. Jones, K. Ritz, ... W. H. van der Putten, Eds.), *Soil Ecology and Ecosystem Services*. Oxford University Press.

607

608 Wardle, D. A. (2002). *Communities and ecosystems : linking the aboveground and belowground components*. Princeton University Press.

609

610 Waring, S. A., & Bremner, J. M. (1964). Ammonium Production in Soil under Waterlogged Conditions

611 as an Index of Nitrogen Availability. *Nature*, 201(4922), 951–952. doi:10.1038/201951a0

612 Wienhold, B. J. (2007). Comparison of Laboratory Methods and an In Situ Method for Estimating

613 Nitrogen Mineralization in an Irrigated Silt-Loam Soil. *Communications in Soil Science and Plant Analysis*, 38(13–14), 1721–1732. doi:10.1080/00103620701435498

614

615 Yoder, B. J., & Waring, R. H. (1994). The normalized difference vegetation index of small Douglas-fir

616 canopies with varying chlorophyll concentrations. *Remote Sensing of Environment*, 49(1), 81–

617 91. doi:10.1016/0034-4257(94)90061-2

618
## Figure Captions

**Figure 1:** Environmental turnover (box 1) affects turnover of multiple ecosystem functions (box 3) directly and indirectly through multi-trophic β-diversity (box 2) along a 1000m elevation gradient (box 4). In our example (turnover between plots A, B and C, box 2 and 4) total multi-trophic β-diversity is high between A and B but low between A and C; horizontal β-diversity is high between plots A and B for mites and springtails and between A and C for plants but low between B and C for springtails (community C is nested in B); vertical β-diversity (variation of MOTUs richness across the groups) is high between B and C but low between A and C.

**Figure 2:** Structure of the tested path models.

**Figure 3:** Total multi-trophic β-diversity mediates the effect of environmental turnover on the turnover of multiple ecosystem functions. Only direct standardized effects are represented. Direct standardized effect of total multi-trophic β-diversity is highlighted with a grey area. Significant and non-significant standardized effects are represented with black and grey points and their associated confidence intervals, respectively. Direct standardized effects of climatic and soil turnover on the vertical β-diversity are represented in appendix (Figure S6).

**Figure 4:** Vertical β-diversity does not mediate the environmental turnover effect on the turnover of multiple ecosystem functions. Only direct standardized effects are represented. Direct standardized effect of vertical β-diversity is highlighted with a grey area. Significant and non-significant standardized effects are represented with black and grey points and their associated confidence intervals, respectively. Direct standardized effects of climatic and soil turnover on the vertical β-diversity are represented in figure S6.

**Figure 5:** Horizontal β-diversities of several trophic groups mediate the environmental turnover effect on the turnover of multiple ecosystem functions. Only direct standardized effects are represented.
represented. Significant and non-significant standardized effects are represented with black and grey points and intervals, respectively. Direct standardized effects of climatic and soil turnover on the horizontal β-diversities and on the multiple ecosystem functions are in figure S5 & S6.

**Figure 6:** Integrated path model highlighting the direct and indirect effects of climate and soil turnover on turnover of ecosystem multifunctionality. The size of the arrows is proportional to the size of the associated standardized path coefficients (only for significant paths). Dotted grey lines represent non-significant paths. Paths with double arrows represent correlations.
<figure><img src="image_2.png" /><figcaption>Figure 1</figcaption></figure>

652

653

654
655

<figure><img src="image_3.png" /><figcaption>Figure 2</figcaption></figure>

656

657

658

659

660
<figure><img src="image_4.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6</figcaption></figure>