Energy and physiological tolerance explain multi-trophic soil diversity in temperate mountains

Aim: Although soil biodiversity is extremely rich and spatially variable, both in terms of species and trophic groups, we still know little about its main drivers. Here, we contrast four long- standing hypotheses to explain the spatial variation of soil multitrophic diversity: energy, physiological tolerance, habitat heterogeneity and resource heterogeneity. Methods: We built on a large- scale observatory across the French Alps (Orchamp) made of seventeen elevational gradients ( ~ 90 plots) ranging from low to very high altitude (280– 3,160 m), and encompassing large variations in climate, vegetation and pedological conditions. Biodiversity measurements of 36 soil trophic groups were obtained through environmental DNA metabarcoding. Using a machine learning approach, we assessed (1) the relative importance of predictors linked to different ecological hypotheses in explaining overall multi- trophic soil biodiversity and (2) the consistency of the response curves across trophic groups. Results: We showed that predictors associated with the four hypotheses had a sta-tistically significant influence on soil multi- trophic diversity, with the strongest support for the energy and physiological tolerance hypotheses. Physiological tolerance explained spatial variation in soil diversity consistently across trophic groups, and was an especially strong predictor for bacteria, protists and microfauna. The effect of energy was more group- specific, with energy input through soil organic matter strongly affecting groups related to the detritus channel. Habitat and resource heterogeneity had overall weaker and more specific impacts on biodiversity with habitat heterogeneity affecting mostly autotrophs, and resource heterogeneity affecting bacterivores, phytophagous insects, enchytraeids and saprotrophic fungi. Main Conclusions: Despite the variability of responses to the environmental drivers found across soil trophic groups, major commonalities on the ecological processes


| INTRODUC TI ON
With the ever-increasing availability of biodiversity information, a global synthesis on the major ecological determinants of broadscale biodiversity patterns is starting to emerge (Belmaker & Jetz, 2015;Braga et al., 2019;Pontarp et al., 2019;Thuiller et al., 2020).
This general understanding is pivotal to forecast how biodiversity responds to natural and anthropogenic changes (McGill et al., 2015;Urban et al., 2016). Yet, most of the empirical support is grounded on specific aboveground macroorganisms, in particular vertebrates and plants. Comparatively, soil biodiversity has been largely less studied (Guerra et al., 2020), although it represents one quarter of global diversity and is essential for decomposition, nutrient cycling or carbon sequestrations Wagg et al., 2014). Therefore, it remains unclear whether the ecological hypotheses that hold true for aboveground systems, such as the energy or the habitat heterogeneity hypotheses, also apply to the massive bulk of belowground biodiversity (Bardgett et al., 2005;Decaëns, 2010).
Historically, the complexity of studying the soil compartment, for example, complex physical structure (Young & Crawford, 2004), taxonomic impediment (Decaëns, 2010), scale of approach (Bardgett et al., 2005;Ettema & Wardle, 2002;Thakur et al., 2020), has hindered the integration of soil biodiversity into a broader ecological hypothesis testing framework. Yet, our ability to study soil biodiversity at large spatial scales is constantly improving with joint taxonomic efforts, the development of new sampling technologies (e.g. Oliverio et al., 2020 for protists), unveiling their environmental drivers. Yet, whether soil biodiversity at all its taxonomic and trophic levels responds to the same ecological drivers as aboveground diversity and follows similar trends remains to be tested. For such tests, the integration of spatial scales and the scale at which organisms are analysed together is pivotal (Thakur et al., 2020;White et al., 2020).
Among the hypotheses formulated to explain the spatial variation of biodiversity, theory and support from empirical studies on plants and other aboveground organisms have led to four major ecological hypotheses: the "energy hypothesis", the "physiological tolerance hypothesis", the "habitat heterogeneity hypothesis" and the "resource heterogeneity hypothesis" (Figure 1). Yet, these hypotheses have been seldom tested in a single framework for soil organisms (Decaëns, 2010;Thakur et al., 2020), and even less at the scale of the whole soil biota. Observing diversity patterns of soil organisms in nature, that is, the relationship between various relevant predictors and soil diversity, is a first step to test whether these ecological hypotheses apply to the wide range of soil organisms (Shade et al., 2018).
The "energy hypothesis" predicts a positive relationship between diversity and energy. An increasing amount of energy (i.e. thermic, solar or chemical) promotes diversity across trophic levels by increasing speciation rates and/or the number of species populations, and thereby reducing local extinction (Evans et al., 2005;Wright, 1983).
An extension of the hypothesis predicts a hump-shaped relationship with a decrease in diversity at high energy levels due to exclusive competition (Mittelbach et al., 2001). Plant productivity is traditionally used as a primary energy measure, because it accounts for water limitations in the transformation of solar energy into available resources, and because plants are the main basal resource (primary producers) for aboveground organisms (Currie et al., 2004;Evans et al., 2005). Yet, in the soil compartment, soil organic matter (SOM) is also a major source of energy fuelling the soil food web (Moore et al., 2004). The local amount and content of SOM is driven by multiple drivers such as plant community composition, climate or parent material (Wiesmeier et al., 2019), and not only by plant productivity.
Considering both solar energy and SOM, hereafter referred as primary and secondary energy, respectively, is thus essential to test the energy-diversity relationship for the soil biota. Therefore, since most soil organisms are thought to be weakly limited by competition due to their limited mobility and the complexity of the soil matrix (Ettema & Wardle, 2002;Wardle, 2006), it could be expected that soil diversity increase monotonously with available energy.
The "physiological tolerance hypothesis" states that favourable environmental conditions support higher biodiversity because a wider range of strategies can persist under such conditions (i.e. tighter niche packing), while only a few well-adapted species can tolerate stressful conditions (Currie et al., 2004;Spasojevic & Suding, 2012). structuring soil biodiversity emerged. We conclude that among the major ecological hypotheses traditionally applied to aboveground organisms, some are particularly relevant to predict the spatial variation in soil biodiversity across the major soil trophic groups.

K E Y W O R D S
environmental DNA metabarcoding, French Alps, macroecology, random forest, soil biodiversity, trophic groups Temperature is one of the most acknowledged factors constraining the "thermal niche" of organisms. Yet, compared to aboveground temperature, soil temperatures are buffered making it more difficult to isolate its effect on soil biodiversity. For example, in mountain environments, soil temperature is strongly regulated by snow cover and duration (Carlson et al., 2015). In the absence of snow, soil frost might impact the structure and activity of soil communities (Schostag et al., 2019;Sulkava & Huhta, 2003). In addition, soil organisms often rely on other abiotic conditions such as water availability, heavy metal content and pH that can generate stressful conditions at extreme values, for example, drought, toxicity, acidity (Gans, 2005;Xu et al., 2012). Indeed, soil pH is recognized as a major driver of soil microorganisms diversity (Fierer & Jackson, 2006).
While the stressful environmental factors may differ, the general response form to stress should be the same for above and belowground diversity.
The "habitat heterogeneity hypothesis" postulates that increasing habitat heterogeneity provides larger niche space or dimensionality that can be finely partitioned and sustain more coexisting species (Stein et al., 2014;Tews et al., 2004). Traditionally, the "habitat heterogeneity hypothesis" is tested at the landscape scale where biodiversity increases with habitat or vegetation diversity (Stein et al., 2014). However, soils can harbour a high degree of heterogeneity at much smaller grains than those considered aboveground (Young & Crawford, 2004), and this partly explains their remarkably high biodiversity (Ettema & Wardle, 2002;Nielsen et al., 2010). On a microscale, habitat heterogeneity can be structural, that is, associated with the size distribution of the pores, which is controlled by soil texture and compaction (i.e. bulk density). Pore size distribution varies within and between soil types, and can influence habitat conditions by modulating nutrient availability, gas diffusion and soil water holding capacity (Ranjard & Richaume, 2001;Six et al., 2004), parameters that may affect the diversity of soil organisms (Nielsen et al., 2010;Xia et al., 2020) and their interactions (Erktan et al., 2020). The effects of soil texture and compaction on the diversity might vary between soil organisms with different sizes or life-history strategies (Seaton et al., 2020) or whether there are ecosystem engineers able to modify the soil structural properties (Decaëns, 2010;Six et al., 2004).
The "resource heterogeneity hypothesis" follows the same rationale as the habitat hypothesis. An increase in resource heterogeneity can lead to an increase in diversity (Steiner, 2001;Heidrich et al., 2020;Dal Bello et al., 2021). We acknowledge that resource heterogeneity can be intrinsically linked to the habitat F I G U R E 1 Overview of the four big ecological hypotheses and theoretical predictions tested in this study within the soil biodiversity context. Each hypothesis is introduced in a coloured box, the predictors used to represent each hypothesis are given at the end of the boxes in a frame heterogeneity, which makes it difficult to separate them. As for aboveground, soil basal resources can take different forms, but their heterogeneity can be well approximated by plant functional diversity since it explains variation in SOM composition, type of potential mycorrhiza, root exudates and direct resources for phytophages (Anderson, 1978;Eviner & Chapin, 2003;Hooper et al., 2000). For higher trophic level groups (secondary and tertiary consumers), the diversity in potential prey might be taken as a proxy for resource heterogeneity.
Here, we tested the above outlined macroecological biodiversity hypotheses and estimated their relative importance in explaining soil biodiversity patterns across most soil trophic groups.
We built on a large-scale observatory network across the French Alps (Orchamp) that provides soil biodiversity measurements from environmental DNA metabarcoding across seventeen elevational gradients ranging from low to very high altitude (280-3160 m), and harbouring very contrasting climatic, vegetation and pedological conditions ( Figure S1). Mountainous systems are well suited to test empirically large-scale drivers of biodiversity as they include wide ranges of environmental conditions and high biotic turnover over a reduced spatial scale (McCain & Grytnes, 2010). Instead of focusing on specific taxonomic orders, we followed a multi-trophic approach to test the above hypotheses on most trophic groups representative of soil ecosystems. After selecting the predictors related to the ecological hypotheses, we used a machine learning approach to account for complex interactions between predictors and soil biodiversity and corrected for remaining spatial dependencies that may originate from processes that have not been considered, such as missing abiotic factors or dispersal limitations. More specifically, we used biodiversity patterns to assess (1) the relative importance of predictors linked to different ecological hypotheses in explaining overall multi-trophic soil biodiversity and (2) the consistency of the response curves across trophic groups.

| Study site and sampling design
The data come from the French Alps long-term observatory, Orchamp (www.orcha mp.osug.fr, Appendix S1), made of multiple elevational gradients distributed across the whole French Alps (ca. 40,500 km 2 ) and representative of the environmental conditions of the region. Each elevational gradient has a homogenous exposure and slope, and consists of four to nine 30 × 30 m plots separated by 200 m of altitude, on average. In this study, we used data gathered from 2016 to 2018, corresponding to 17 gradients ( Figure S1), 90 plots and 540 soil samples. Plant species abundances were quantified at the vegetation peak (mostly in July or August) along a linear transect crossing each plot using the pin-point method (Jonasson, 1988). A second 4-m-wide transect was dedicated to soil sampling at the end of the summer season. Soil was sampled from 3 subplots (2 × 2 m) selected across the transect where we collected around ten soil cores of 5 cm in diameter that were separated into two soil layers, that is, surface (ca. 1-8 cm depth) and subsurface (ca. 8-16 cm depth), which could be differentiated in most cases by a change in the colour. The ten soil cores were pooled together and homogenized by separating the two layers to make a biological sample per soil layer per subplot, to obtain a total of six samples per plot.

| Soil sample processing
Each soil sample was separated into two components. The main part was sieved at 2mm and used to measure soil physicochemical properties (soil pH, SOM content and soil C/N) as described in (Martinez-Almoyna et al., 2020). The other part was used for environmental DNA, where DNA was extracted from a 15 g aliquot and processed in the field using the procedure described in Taberlet et al. (2012), Taberlet et al. (2018). We used six DNA markers to have a complete overview of the soil biota, including two universal markers (euka02 for eukaryote, bact01 for bacteria) and fourth clade-specific markers (fung02 for fungi, inse01 for insect, olig01 for oligochaete and coll02 for collembola). Information on the markers and molecular analyses including PCR, library preparation and sequencing steps are detailed in Appendix S2. A standardized bioinformatic pipeline was then applied (Calderón-Sanou et al., 2020), using the OBITools software (Boyer et al., 2016) and the R package "metabaR" , to remove contaminants and errors and to get the taxonomic composition in terms of Molecular Operational Taxonomic Unit (MOTU) of each sample (Appendix S2).

| Diversity of trophic groups
The obtained MOTUs were classified into 36 trophic groups. We chose to distinguish not only trophic levels but also phylogenetic distant groups of the same trophic level, as they may have different preys/predators or exhibit different resource acquisition strategies  Table S1. The MOTU diversity of each trophic group was estimated per sample using the exponential of the Shannon entropy (i.e. Shannon diversity), which represents "the effective number of MOTUs" as it penalizes rare sequences that could be artefacts in eDNA data. Shannon diversity leads to more robust ecological results and to diversity estimates that are more similar to those assessed from conventional sampling approaches (Calderón-Sanou et al., 2020).

| Environmental predictors
We used two environmental predictors to represent each ecological hypothesis (Figure 1), with the condition of having a final set of weakly correlated predictors (see Figure S2 for a visualization of the correlation between all initially considered parameters).

| Energy hypothesis
It was separated into primary (solar energy) and secondary energy (SOM), and two predictors were selected for each category. Solar radiation and the Normalized Difference Vegetation Index (NDVI) were used to represent the primary energy predictors. Solar radiation directly measures the amount of solar energy arriving into the Earth's surface, while NDVI estimates the amount of solar energy that is transformed by photoautotrophic organisms into available resources accounting for water limitations (Evans et al., 2005). We did not add mean annual temperature as sometimes done to represent energy (Clarke & Gaston, 2006) since it was strongly correlated to NDVI ( Figure S2). Solar radiation was calculated per plot as the sum of the daily surface incident direct and diffuse shortwave radiation accumulated over 10 years, from 2008 to 2018. NDVI was estimated from the surface spectral reflectance at a resolution of 250 m from MODIS (Moderate Resolution Imaging Spectroradiometer), available online: https://lpdaac.usgs.gov/produ cts/mod09 q1v00 6/. Raw NDVI times series were pre-processed following Choler (2015), and we kept the mean yearly sum of NDVI greater than 0.2 over 2009-2019, as the final predictor for the analyses measured at the plot level. To represent secondary energy, we used the SOM content and the C/N ratio, measured from the soil samples. The former indicates the total amount of organic matter available in the soil, while the latter is a proxy for nutrient availability or SOM decomposability (Cleveland & Liptzin, 2007), meaning that soils with low C/N rates have potentially more readily available energy than soils with high C/N ratio, if we account for nutrient stoichiometric constraints.

| Physiological tolerance hypothesis
We used soil pH and the freezing degree days (FDD) to represent potential sources of abiotic physiological stress for soil organisms. The pH has been described as an important limiting physiological factor of soil communities (Fierer & Jackson, 2006;Räty & Huhta, 2003).
The FDD summarizes the duration and intensity of ground freezing events and it has been addressed as a good candidate to model the thermal niches (Choler, 2018). FDD was calculated per plot as the annual sum of average daily degrees below zero, modelled within the first soil horizon (1 cm depth) and averaged over 2008-2018.

| Habitat heterogeneity hypothesis
Clay percentage in soil and bulk density were selected to represent the microscale habitat heterogeneity. Clay percentage characterizes the soil texture and thus reflects the granulometry distribution, the aeration, ability of soil to retain water and more globally the physical properties of the soil (Hao et al., 2007;Seaton et al., 2020). Soil texture might affect diversity differently across trophic groups with different sizes or life-history strategies (Seaton et al., 2020;Vreeken-Buijs et al., 1998). For example, the diversity of mesofauna could be expected to increase in coarse-textured soils (i.e. with low clay percentage), where the higher availability of larger pores provides more different habitats to be potentially colonized by these organisms (Vreeken-Buijs et al., 1998). Bulk density reflects soil compaction and porosity as it accounts for the amount of soil per volume unit when removing water and air spaces (Hao et al., 2007). Compact soils, with higher values of bulk density, have relatively lower total pore space and organic matter content, thus providing a lower heterogeneity of habitats. Both variables were measured from a soil pit carried out next to the plot). Three soil replicates were collected with a volumetric cylinder (100 cm3) from the superficial horizon.
They were dried at 105°C for 24 h and sieved to 2 mm. The mass of dry soil (mS) contained in the cylinder as well as the mass of coarse elements greater than 2 mm (mEG) were measured. The formula applied for the calculation of bulk density is as follows (Equation 1), with V cyl for the volume of the cylinder. The bulk density of the three replicates were averaged.

| Resource heterogeneity hypothesis
For decomposers, detritivores and plant symbionts, we used two metrics of plant functional diversity as predictors, that is, the functional richness and the functional divergence (Villéger et al., 2008), calculated for each plot using the R package "FD" (Laliberté & Legendre, 2010). Functional richness represents the total trait space filled by all the plant species present in the community (here the plot). Functional divergence describes how specie's abundances are distributed within the functional trait volume. To estimate these two metrics, we used our own trait measurement values for species

| Spatial structure
Given the hierarchical sampling design of the data (two soil layers within plots within gradients), we accounted for the overall spatial structure of the samples to avoid having spatial autocorrelation issues (Dray et al., 2012). We defined a set of spatial predictors representing the residual spatial structure (i.e. the left-out spatial structure not explained by the environmental predictors) to include in the models. This approach aims to reduce the spatial autocorrelation that could remain in the residuals and to identify potential spatial structures with a strong influence on soil diversity. We did so using Moran's eigenvector maps (MEM), a method based on computing the principal coordinates of a matrix of geographic neighbours (Dray et al., 2006). The obtained eigenvectors are orthogonal and have a straightforward interpretation as each of them represents a spatial pattern at a given scale that can be ranked from broad spatial structures to fine spatial structures. We identified 18 MEM-variables describing significant spatial autocorrelation (only positive eigenvalues, Dray et al., 2006) based on the Euclidean geographic distances between each subplot using the function dbmem from the R package "adespatial" (Dray et al., 2021). MEM 1 to 8 described broad scale spatial structures, while MEM 9 to 18 represented intermediate to fine spatial structures ( Figure S3). To remove the imprint of the environment on these MEMs, we modelled with a random forest each of the 18 MEMs as a function of our environmental predictors and extracted the residuals of these relationships. These residuals thus represented the spatial structure not explained by our environmental predictors (e.g. missing predictors, dispersal limitations). This approach differs from partialling out the spatial component of diversity and compare the pure effect of environment, the pure effect of space and the shared explained variance (Borcard et al., 1992). Here, we argue that space is likely affecting environment and that environment is then affecting biodiversity. The shared explained variance of space and environment is thus relevant for our hypotheses. We treat the pure effect of space as a statistical nuisance as we cannot link it to ecological processes, given that we jointly analyse taxa with very different dispersal abilities. We made sure that it was properly accounted for to avoid residual spatial autocorrelation (Dray et al., 2012).

| Random forest
To model the diversity of each trophic group as a function of the predictors representing our four hypotheses and the residual spatial structure, we used random forest models (Breiman, 2001), which are particularly suited when nonlinear relationships and complex interactions among predictors are expected. Random forest analyses were run with the R package "party" (Hothorn et al., 2006) with the cforest_unbiased function, which avoids bias introduced by heterogeneity in scale and number of categories among predictors (Strobl et al., 2007). The number of trees was set to 1,000 and the number of variables randomly sampled as candidates at each split (mtry) was tuned using the function train of the R package "caret" (Kuhn, 2020; Table S2). Variable importance was estimated as the mean decrease in accuracy using the function varimp. The method allows assessing relative variable importance, by identifying the covariates which, when removed, ensure a significant drop of prediction power (Strobl et al., 2007). It thus avoids any over-fitting and allows sound inference. Overall explained variance (r-square) was calculated by extracting the coefficient of determination between predictions and observations. The shape of the relationship between the diversity and the predictors was assessed with partial dependent plots obtained from the R package "iml" (Molnar et al., 2018), which estimate the marginal effect of a given predictor while accounting for the average effect of the other predictors in the model. We considered that a relationship was relevant, when the predictor had a predictive importance higher than 25%. The predictive importance was assessed by permuting each predictor one by one and then evaluating how the prediction was affected.
A single random forest model was run for each trophic group with the same set of predictors, that is, solar radiation, NDVI, SOM, C/N ratio, percentage of clay, bulk density, two variables corresponding to resource heterogeneity (variable across trophic groups, and excluded for autotrophs) and the 18 residual spatial structure predictors. All analyses were run in the R statistical environment (R Core Team, 2020).

| RE SULTS
We identified 222,739 bacterial and 50,241 eukaryotic (including 5,467 metazoans and 11,115 protists, Figure A2- The predictors underlying the tested ecological hypotheses explained a significant part of the spatial variation of diversity of most trophic groups. The overall explained variance varied from 29% for detritivorous insects to 79% for arbuscular mycorrhizal fungi ( Figure 2a, Table S2). The residual spatial structure explained much less variance than the environmental predictors, confirming the relevance of the latter to predict soil biodiversity. Only the diversity of predatory and phytophagous insects, and photoautotrophic protists was better explained by pure broad residual spatial structures than by the environment (MEM7, Figure S3).
We found that predictors associated with the energy and the physiological tolerance hypotheses were generally the most im- In general, we found that the partial response curves of diversity to predictors were consistent across most soil trophic groups ( Figure 4) and in agreement with predictions ( Figure 1), with some few exceptions. The diversity of most trophic groups including zooparasites protists and fungi, metazoans consumers and ectomycorrhizal fungi, strongly increased with NDVI, but decreased for photolithoautotroph bacteria, phytophagous protists and earthworms ( Figure 4a). The steepest changes in soil diversity across the NDVI gradient occurred in the transition from forest (high NDVI) to alpine grasslands (low NDVI). Groups for which diversity strongly increased with solar radiation included zooparasite bacteria, phytophagous protists and earthworms. All trophic groups primarily feeding on detritus positively increased in diversity with SOM ( Figure 4b).
The diversity of several groups was also influenced by the C/N ratio: diversity decreased for herbivorous and bacterivorous nematodes, and root endophyte and arbuscular mycorrhizal fungi, but increased for ectomycorrhizal fungi and fungivorous nematodes (Table S3). With the exception of rotifers and tardigrades, all trophic groups responding to pH increased in diversity in more alkaline soils (Figure 4c). This positive relationship had a sigmoid form for all groups, but both the inflection points and associated slopes strongly varied across trophic groups. Saprotrophic, root endophytes and phytoparasitic fungi, and also photolithoautotrophic bacteria were positively affected by the soil clay content, and chemolithoautotrophic bacteria were positively affected by soil bulk density (Table   S3). All phytophagous insects, saprotrophic fungi and bacterivore groups responded positively to resource heterogeneity, that is, plant functional richness and bacteria diversity respectively (Figure 4d).

| DISCUSS ION
Testing ecological hypotheses has largely contributed to our understanding on how biodiversity is structured on Earth. However, generality can only be claimed if a significant part of biodiversity is covered. Here, we add an important missing piece to the general picture by testing several major ecological hypotheses simultaneously on the majority of trophic groups inhabiting the soil and along sharp environmental gradients which allow some generalization to be made. Our results confirm that the main environmental drivers of soil biodiversity are variable across soil trophic groups and depend on their resource or physiological requirements. Yet, we also find major commonalities in the ecological processes structuring soil biodiversity. Overall, the energy and physiological tolerance hypotheses had the strongest support from soil multi-trophic biodiversity.
Our results are in agreement with previous studies finding that an increase in primary energy increases the diversity of soil organisms such as protists (Oliverio et al., 2020), metazoans (Peters et al., 2016), soil predators (Binkenstein et al., 2018) and fungi (Tedersoo et al., 2014, Figure 2b). Our results also reveal that secondary energy, related to soil organic matter, has a positive effect on soil biodiversity, especially for fungivorous and detritivorous animals, in agreement with earlier work (Binkenstein et al., 2018;Canedoli et al., 2020;Caruso et al., 2019). We found that the relative importance between primary and secondary energy varies across trophic groups, with no clear trends across trophic levels, suggesting that both energy channels are at play across the soil food web. However, some groups responded to specific energy predictors in a way that differs from the predictions of the "energy hypothesis" (Figure 1).
For example, the diversity of earthworms, phytophagous fungi and photolithoautotroph bacteria decreased with increasing NDVI.
Part of these divergent trends between diversity and NDVI might be explained by the transition from forest to grassland in the NDVI gradient in our study system, for example, alpine grassland soils are more suitable for autotrophic bacteria adapted to high elevation stressful conditions (Guo et al., 2015). Otherwise, a negative interaction between ectomycorrhizal fungi and phytophagous fungi could explain the decrease in diversity of the latter (Figure 4a). Indeed, ectomycorrhizal fungi can provide protection against pathogens to their plant hosts, thus reducing the incidence of phytophagous fungi and their diversity (Antunes & Koyama, 2017;Wang et al., 2019).
Other divergent, but not unexpected, trends were found along the TA B L E 1 Information on the environmental DNA data characterizing each trophic group, including the DNA marker used to sample each group and the final number of reads, families (orders for protists), MOTUs and Shannon diversity obtained in total across the French alps and per sample (mean ± SD)

DNA marker
Total reads (per sample)

Total families/orders (per sample)
Total MOTUs (per sample)  (Lindahl & Tunlid, 2015). An increase in ectomycorrhizal fungi diversity could presumably cascade on fungivore nematodes diversity. Furthermore, while we show that energy has mainly a positive influence on soil biodiversity, the underlying mechanisms remain to be tested. For example, the more individual hypothesis states that greater energy availability allows a community to contain a larger number of individuals, and hence of a larger number of species with viable population size (Wright, 1983). Quantifying species abundance or biomass would be needed to test this hypothesis, but this information is unfortunately not yet available with eDNA metabarcoding data (Taberlet et al., 2018), and would be extremely challenging to obtain for the wide range of organisms studied here.
Physiological tolerances, mainly to soil pH, were also a strong predictor of the diversity of soil organisms, especially for organisms living in the aqueous phase of the soils. Indeed, in the study system, the diversity of groups of bacteria, protists and microfauna was more constrained by pH-induced stress rather than limited by energy or habitat and resource heterogeneity (Figure 3b), in accordance with previous studies highlighting the importance of pH for soil microbes (Fierer & Jackson, 2006;Karimi et al., 2018) and invertebrates Räty & Huhta, 2003). The sigmoid trend observed between diversity and pH might correspond to the first half of the humpback curve expected from the theory (Figure 1).
Indeed, our sampling had relatively few sites with alkaline soils, and did not include soils with pH >8, levels from which other studies have observed a decrease of diversity (e.g. Fierer & Jackson, 2006).
Our results revealed consistent decreases of diversity in more acidic soils, but also different tolerance thresholds across soil trophic groups. The strong effect of soil pH might also be the sum of multiple linked factors not considered in this study including bedrock type and plant communities (Roy et al., 2013). Contrarily, FDD had a minor effect on soil biodiversity. Limited effect of freezing events on soil biodiversity has previously been reported, and may result from the frost resistance (Männistö et al., 2018;Stres et al., 2010) or the rapid recovery of soil communities (Sulkava & Huhta, 2003).
Theoretically, this low importance could be due to a scale mismatch between the measured soil communities (subplots are 4m 2 large) and the climatic data resolution (~300m). However, between the available in situ temperature HOBOs and the climatic data used here showed very consistent patterns, rendering the scale mismatch hypothesis unprobeable. Otherwise, a change in composition or activity, without changes in local diversity, might also have occurred and remains to be tested (Schostag et al., 2019;Stres et al., 2010).
The "habitat heterogeneity" and the "resource heterogeneity" hypotheses weakly explained the spatial variation in diversity of soil trophic groups compared to "energy" and "physiological tolerance",

TA B L E 1 (Continued)
with notable exceptions. Saprotrophic, root endophytes and phytoparasitic fungi, as well as autotrophic bacteria were highly affected by habitat heterogeneity. We found that these groups tended to be more diverse in fine-textured soils (higher clay percentage), which usually exhibit greater water retention capacity but also more recalcitrant and stable organic matter (Ranjard & Richaume, 2001;Six et al., 2004). Previous studies have shown that soil texture can influence bacterial and fungal diversity, with subgroups of taxa responding differently to the proportion of soil particles (i.e. clay, sand, silt) (Karimi et al., 2018;Seaton et al., 2020;Xia et al., 2020). Our results showed that such differences are also visible when considering different trophic groups of fungi and bacteria.
The importance of "habitat heterogeneity" could be expected to vary across soil trophic groups, as the spatial scale at which heterogeneity is perceived by organisms of different sizes or different lifestyles can be highly variable (Heidrich et al., 2020). Here again, perhaps the scale at which we measured heterogeneity was not relevant for some specific groups. When looking at the effect of resource heterogeneity, prey's diversity was remarkably important for bacterivores. Strong associations between bacterivore protists F I G U R E 2 Relative importance of competing hypotheses in explaining the alpha diversity of soil trophic groups. (a) Total r-squared of the random forest model for each trophic group. Colours represent the relative importance of the environmental versus the spatial predictors. Environmental predictors correspond to all the biotic and abiotic variables used to test the ecological hypotheses, and spatial predictors correspond to the residuals of the spatial structure when removing the effect of the environment. (b) Relative importance of the environmental predictors used to test the ecological hypotheses (colour key). The relative variable importance is the mean decrease in squared error, rescaled to sum the total r2 (a) or 1 (b). Letters correspond to broad taxonomic groups: Bacteria (B.), Protozoa (P.), Metazoa (C.: Collembola, I.: Insects, M.: Mites, N: Nematodes) and Fungi (F.). Symbols indicate the size category for fauna groups and bacteria diversity have been recently reported (Oliverio et al., 2020;Xiong et al., 2021), and could indicate a degree of trophic specialization in bacterivorous protists. Co-variation in diversity might also indicate shared habitat preferences between protists (or nematodes) and bacteria, but our results and previous studies point to noticeable differences in the factors shaping the diversity of these groups (Oliverio et al., 2020;Xiong et al., 2021). Moreover, the strong response of saprotrophic fungi to plant functional diversity could be explained by a trophic specialization, in accordance with a recent study showing a high degree of specialization to specific soil and litter compounds for some saprotrophic fungi (Algora Gallardo et al., 2021). The significant association does not necessarily imply the realization of a trophic interaction, but it is a first step in assessing whether such interactions exist, leave signals in diversity distribution and can give us insights into the degree of food speciation in the focus trophic group.
To conclude, our near-complete coverage of soil biodiversity across trophic groups and across large and steep environmental gradients provides consistent and novel insights on the macroecological rules shaping the distribution of belowground biodiversity. Building on the efficiency of environmental DNA analyses combined with the wealth of existing knowledge on soil organisms, we showed that energy and physiological tolerance are the most plausible hypotheses to explain the spatial distribution of soil diversity at a regional scale. Interestingly, we found strong commonalities between trophic groups in their response to environmental drivers that should be later compared to aboveground organisms living in the same locations (e.g. ground-dwelling arthropods, pollinators or herbivores). Should belowground and aboveground compartments respond differently to environmental drivers, it will complexify their management under humaninduced pressures. Finally, identifying how these patterns in local diversity translate into compositional changes and interaction network structuration in space will be of crucial importance to understand soil biodiversity assembly and how it might be affected by ongoing environmental changes.

F I G U R E 3
Boxplots of the relative importance of ecological hypotheses by trophic position and body size category. Relative importance of the four ecological hypotheses tested in this study across groups categorized by trophic position (a) or body size category (b). The values of relative importance correspond to the mean decrease in squared error from the random forest per trophic group, rescaled to sum the total r-square

ACK N OWLED G EM ENTS
This project was possible thanks to the many people that participated in the field work and in the laboratory procedures, including students, interns, postdoctoral fellows, field assistants, park managers and rangers. We thank Cindy Arnoldi for conducting

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data supporting the results of this study are publicly available in the DRYAD repository https://doi.org/10.5061/dryad.q573n 5tm6.

F I G U R E 4
Predicted diversity of soil trophic groups as a function of the environmental predictors representing the ecological hypotheses. Partial dependence plots showing the marginal effect of the predictors representing the ecological hypotheses on the diversity of soil trophic groups based on the random forest model results. The predictors represented are (a) NDVI, with the colours corresponding to the transition from alpine grassland (yellow) to forest (green), (b) Soil organic matter, (c) pH and (d) plant functional richness and Bacteria MOTUs diversity. Diversity was standardized by the maximum diversity observed per group to have a comparable scale across groups. Only groups for which the predictors had a predictive importance higher than 25% were represented. The predictive importance of the predictor was assessed by permuting each predictor one by one and then evaluating how the prediction was affected. Taxonomic groups are abbreviated as Bacteria (B.), Collembola (C.), Earthworms (Earth), Fungi (F.), Insects (I.), Mites (M.), Nematodes (N.) Protozoa (P.), Rotifera (Rotif) and Tardigrada (Tardi). The rest of the trophic component is abbreviated as arbuscular mycorrhizal (arb), bacterivore (bac), detritivore (det), ectomycorrhizal (ect), epigeic (epi), euedaphic-hemiedaphic (e-h), fungivore (fun), herbivore (her), heterotroph (het), omnivore (omn), photolitoautotroph (pho), phytoparasite or phytophageous (phy), predator (pre), protistivore (pro), saprotroph (sap) and zooparasite (zoo)