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

Handling Editor: Charles Fox 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 1,000 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 soil properties and climatic gradients in alpine ecosystems. 4. By demonstrating how saprophytic fungi and their associated trophic groups can offset the direct responses of multiple ecosystem functions to environmental change, our study highlights the paramount importance of multi‐trophic diversity for better understanding ecosystem multifunctionality in a changing world.


| INTRODUC TI ON
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;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).
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, that is 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 .
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;. 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, for example 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 et al., 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 1,000 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 analysed the relative contributions of horizontal β-diversities of different trophic groups to the spatial turnover of multiple ecosystem functions.

| Study site
The study was carried out in the French Alps, along a continuous elevation gradient (1,750-2,725 m, Arves Massif, 45.12°N, 6.40°E) in a cow-grazed pasture. Subalpine grasslands dominate the bottom of the gradient while sparsely vegetated alpine meadows characterize higher elevations. Ten sites were sampled in summer 2012 along the same south-facing slope, and they were separated by an 100 m elevation difference and an average distance of 340 m. Each site consisted of two 10 m × 10 m plots with homogeneous vegetation (see Chalmandrier et al., 2017, for more details).

| Environmental variables
We considered four bioclimatic variables known to be strong drivers of community assembly in mountain systems (Körner, 2003): soil temperature, solar radiation, growing season length (GSL) and number of frost days. Growing season length and annual number of frost days were derived from Landsat 7 and 8 imageries (Carlson, Choler, Renaud, Dedieu, & Thuiller, 2015). Solar radiations were calculated with the area solar radiation tool in ArcGIS (version 10.2, 2013; Redlands, CA, USA) using a 2-m LiDAR digital elevation model and the clear sky model set to a sky size of 2,800 pixels (Carlson et al., 2015). Soil temperature was measured in the field.
To capture small-scale variation, we added local topographic variation using elevation, slope, the topographic wetness index (TWI) and the topographic position index (TPI) with a 50 cm resolution digital elevation model derived from the airborne LiDAR data acquired the year of sampling.
Finally, we measured soil properties from three soil cores (10 cm depth, 5 cm diameter) taken in a 1.25 m × 1.25 m subplot at the centre of each plot. Two soil cores were weighed, 5 mm sieved and pooled for quantifying soil moisture, soil organic matter content, pH, soil texture using standard protocols (Robertson, Coleman, Sollins, & Bledsoe, 1999), and total soil N and C content using an elemental analyzer (FlashEA1112, Thermo Fisher Scientific). We added 100 ml of distilled water to saturate the third soil core and calculate the water-filled pore space (WFPS).
The soil core was then air-dried, sieved and grounded to calculate gravel weight, total porosity and apparent density (Legay et al., 2014).
We found that the first two axes of a principal component analyses run for all produced (and normalized) environmental variables F I G U R E 1 Environmental turnover (box 1) affects turnover of multiple ecosystem functions (box 3) directly and indirectly through multi-trophic βdiversity (box 2) along a 1,000 m elevation gradient (box 4). In our example (turnover between plots A, B and C, boxes 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 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 S1) captured, respectively, 48.1% and 19.4% of the environmental variability. We then split the environmental variables in two groups corresponding to the major environmental gradients of the study area: a climatic gradient and a gradient of soil properties.

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

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

| Sampling
To get an estimate of the diversity found in each of the sampled plot, we analysed soils samples 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 universal markers, one amplifying all eukaryotes (v6-v7 region of the 18S rRNA gene) and one amplifying all bacteria (v5-v6 region of the 16S rRNA gene), were used to obtain a general overview of the soil multi-trophic composition of the plots. We also used two additional markers focused on fungi (internal transcribed spacer 1) and vascular plants (chloroplast trnL-P6 loop). Molecular analyses and data curation are presented in Ohlmann et al. (2018). We pooled the 21 samples per plot to obtain a single community of molecular taxonomic units (MOTUs) per plot and converted the data into presence-absences. The curated sequencing data as well as associated metadata are available on the Dryad Digital Repository under accession https ://doi.org/10.5061/ dryad.5b58400. Following Ohlmann et al. (2018), we grouped the MOTUs based on their taxonomic affiliation, shared trophic resources and main functions, such as organic matter decomposition. Hereafter, we will refer to them as trophic groups, while other definitions like trophic guilds or tropho-functional groups could have been given.

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

| Different components of multi-trophic βdiversity
First, we computed the true turnover (Baselga, 2010) of the whole soil community (i.e. ignoring trophic groups) as a measure of the total multi-trophic β-diversity using the Simpson dissimilarity index between each pair of plots (Equation 1).
where k and l are the focal plots (such as 1 ≤ k ≤ K ∧ 1 ≤ l ≤ K, K the total number of plots), P = P ki 1≤k≤K 1≤i≤n is the community matrix of the whole soil community (n is the total number of MOTUs and P ki ∈ 0,1 denotes the presence of the MOTU i in plot k). P k. is the community vector of the plot k and � � � is the number of shared MOTUs between plot k and l.
Second, we calculated the true turnover (Baselga, 2010)  ( where k and l are the focal plots (such as 1 ≤ k ≤ K ∧ 1 ≤ l ≤ K, K the total number of plots), g is the focal trophic group, P g = P g ki 1≤k≤K 1≤i≤ng is the community matrix of the group g (n g is the total number of MOTUs, and P g ki ∈ 0,1 denotes the presence of the MOTU i in plot k). P g k. is the community vector (of the group g) of the plot k, and � � � P g ki denotes the richness of the plot k for the group g.
And finally, we calculated the richness variation across trophic groups as a measure of the vertical β-diversity using the Bray-Curtis (BC) dissimilarity index between each pair of plots (Legendre & Legendre, 1998): where G is the total number of trophic groups.

| Productivity
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 1 m × 1 m 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 airborne 15 cm 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.

| N-cycling
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  (Loomis, Ruess, Sveinbjörnsson, & Kielland, 2006). Finally, potential nitrogen mineralization rates were used to estimate the ability of the ecosystem to mobilize soil organic matter and supply mineral N (Legay et al., 2016), and measured based on anaerobic incubations of fresh soil subsamples (dark, 7 days, 40°C), during which organic N was mineralized and accumulated as NH 4 + (Waring & Bremner, 1964;Wienhold, 2007). The difference between NH 4 contents in a given sample before (t1) and after the anaerobic incubation (t2) gave PNM = [(NH 4 + -N)t2(NH 4 + -N)t1]/dw/7 days. N-cycling turnover was then defined as the Euclidean distance between each pair of plots for C/N, N-TDN/N, potential mineralization and N-DON/ (N-NO 3 + N-NH 4 ).

| N-leaching
As a proxy for nitrogen losses through leaching, we measured potential ammonium and nitrate leaching from the percolate.
Percolated water was filtered through a Whatman filter paper n°42 (2.5 µm pore size) and analysed using a photometric analyzer to calculate potentially leached NH 4 + -N and NO 3 − -N (Legay et al., 2016). N-leaching turnover was then defined as the Euclidean distance between each pair of plots for potential ammonium and nitrate leaching.

| Multifunctionality
We calculated multifunctionality turnover as the Euclidean distance between each pair of plots 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.

| Statistical analyses
To address our first two objectives, we performed a set of path analyses using the lavaan package in R ( 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 (Fourtune et al., 2018). This approach is suitable when modelling pairwise distances where the elements of one row/column are not independent of each other (Fourtune et al., 2018). Parameters are significant when the confidence intervals do not include zero.

| RE SULTS
Turnover in soil properties was a consistent and strong direct predictor of multifunctionality and of the three ecosystem functions, while climatic turnover was not (Figures 3 and 4). Total multi-trophic β-diversity was a significant predictor of productivity and multifunctionality turnover and almost as important as soil properties turnover. However, it did not explain N-cycling and N-leaching turnover ( Figure 3).
In contrast, when using vertical β-diversity instead of total multi-trophic β-diversity, the diversity effect vanished and only soil properties turnover had a significant effect on turnover of ecosystem functions, including multifunctionality ( Figure 4). In other words, the effect of total multi-trophic β-diversity on turnover of multiple ecosystem functions was not driven by vertical β-diversity.
Finally, focusing on the independent effects of horizontal βdiversities of the different trophic groups, the most important F I G U R E 2 Structure of the tested path models F I G U R E 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 multitrophic β-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 total β-diversity are represented in appendix ( Figure S6) predictors of multifunctionality and productivity turnover were soil properties turnover and horizontal β-diversities of saprophytic fungi, followed by symbiotic fungi, herbi-fungivore nematodes, pathogenic fungi and bacteria ( Figure 5; Figure S3). For N-cycling turnover, the most important predictors were soil properties turnover and horizontal β-diversities of saprophytic fungi, followed by symbiotic fungi and pathogenic fungi ( Figure 5; Figure S3). No component of multitrophic β-diversity was linked to N-leaching turnover ( Figure 5), for which soil properties turnover was the only significant predictor ( Figure S3). For the non-aggregated ecosystem processes, different subsets of trophic groups were the best predictors, but the general trends remained ( Figure S7).
The integrated path models that considered several trophic groups together, but were built independently for each ecosystem function and multifunctionality, suggested how the pieces of the puzzle can be put together ( Figure 6). Consistent with the results above, soil properties were of primary importance for all considered ecosystem functions and multifunctionality directly, whereas climate was not (Figure 6; Figures S4 and S5). Similar figures emerged when considering each ecosystem process independently ( Figure S7).
Climatic turnover and soil properties turnover influenced horizontal β-diversities of most of the trophic groups ( Figure S6).
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; Figures S4 and S5). Bacteria were linked especially to saprophytic fungi ( Figure 6; Figure S4) and to a lesser extent to symbiotic fungi (Table S1).
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 ( Figures S7c,d).
F I G U R E 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 Brose and Hillebrand (2016)  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.

| D ISCUSS I ON
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. GSL, 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 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, that is 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. (2016) 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 F I G U R E 6 Integrated path model highlighting the direct and indirect effects of climate and soil properties 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 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; Figures S4 and 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 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; Figure  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, adsorbed on soil particles, denitrified, volatilized or assimilated by plants and bacteria.
The balance between these processes mainly depends on abiotic soil conditions which control water and solute transfers, and on the nitrogen demand of plant and micro-organisms (Robson et al., 2007) especially during the productivity peak when plant N-uptake is the highest (Legay et al., 2016). Therefore, weaker effects of each trophic group horizontal β-diversity on the turnover of productivity and N-cycling could be explained by a functional complementary effect between trophic groups of fungi and their specific biotic regulators, and feedbacks between trophic groups and soil. Adding other, and so far, a bit less studied pathways (e.g. P-cycling, Leff et al., 2015), may have uncovered additional direct and indirect links of environmental turnover and functioning turnover via β-diversities.
To conclude, looking at multiple ecosystem functions showed an important mediating role of multi-trophic β-diversity. Our finding particularly supports the paramount importance of decomposers, especially saprophytic fungi, and of their interaction partners, in nitrogen limited alpine systems. In this sense, our study supports recent calls to explicitly integrate known interactions between trophic groups (Thompson et al., 2012) and for temporal monitoring of the ecosystems (Brose & Hillebrand, 2016), two missing steps towards a better understanding of functional complementarity and selection effects of diversity and of their role in multiple ecosystem functions (Hines et al., 2015) and their trade-offs in natural systems.

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