Network analyses
With the aim of analysing spatial variation in soil food web structure, we used beta diversity metrics that account for compositional changes in the food web across samples. We used a set of network diversity indices based on the generalisation of Hill numbers proposed by (Ohlmannet al. 2019) These indices allow to quantify diversity in trophic groups and trophic interactions, varying the weight of the relative abundance of groups and links. Here, we focussed on Shannon diversity (q=1, see Calderón‐Sanou et al. 2020) but also replicated the analyses using species richness (q=0) and Simpson diversity (q=2) to assess the extent to which the observed changes are compositional versus structural. Since interactions in the metaweb are binary (either presence or absent), we approximated the abundance of interactions as the product of the relative abundances of the interacting groups. These diversity indices can be further decomposed into alpha, gamma, and ß-diversity (ßT), and be used to calculate the pairwise dissimilarities for both groups (ßpP) and links (ßpL). ß-diversity (ßT) quantifies the effective number of different communities, based on shared groups (ßTP) or interactions (ßTL), ranging from 1 (indicating identical group or interaction abundance distribution) to the total number of communities, when networks do not share any common group or interaction. The diversity metrics were computed using the R package ’econetwork’ (Mieleet al. 2021).
To quantify changes in the spatial structure of soil food webs across habitats, we calculated pairwise dissimilarity measures for both groups (ßpP) and interactions (ßpL) among all pairs of local food webs. Subsequently, we used UMAP, a non-linear dimension reduction algorithm, to visualize the local food webs in a two-dimensional space based on their group and interaction dissimilarities (McInnes et al.2020). Given the strong structural differences observed in local soil food webs between forests and grasslands (Fig.2, S3), and the high variability in food web composition in shrublands due their intermediate state along the elevational gradient, our analysis focused exclusively on grasslands and forests for further comparison. We compared the relative abundances of trophic groups and classes between forests and grasslands with a Wilcoxon test and used the ‘group-TL-tsne’ network layout from the R package metanetwork (Ohlmann et al.2023) together with edge bundling provided by edgebundle R package (Schoch 2022) to represent the difference network between average food webs in forests vs. grasslands.
We tested the imprint of trophic interactions on community assembly with null models. Firstly, we constructed an average soil food web for each habitat by using the average relative abundances of trophic groups per habitat. We then calculated the ß-diversity (ßT) for both group and link interactions between the two habitats. Next, we built the null model. Since group and interaction diversities are strongly correlated, we aimed at keeping group diversity constant in the null model while modifying interaction diversity. This was achieved by permuting the node labels within the metaweb. Consequently, group ß-diversity remains unchanged as this diversity index is invariant to label permutation. However, interaction ß-diversity was affected since it relies on the product of relative abundances of interacting groups. If dominant interacting groups differ between the two habitats, the observed interaction ß-diversity would be higher than expected under a random distribution of group abundances. To evaluate the significance of our findings, we conducted 3000 permutations and compared the observed diversities with the corresponding null distribution.
Finally, to quantify the relative importance of environmental and spatial distances in explaining food web dissimilarity within habitats (i.e., grasslands and forests), we used Generalised Dissimilarity Models (GDM) (Ferrier et al. 2007). We built two GDM per habitat using trophic group turnover (ßpP) and interaction turnover (ßpL) as response variables, with the R package gdm (Fitzpatrick et al.2022). For the environmental predictors we selected a set of weakly correlated variables (Pearson’s r < 0.5, Fig.S6) representing climatic, soil and vegetation conditions that could influence soil organisms (Table 1). The spatial coordinates of the samples and all the selected environmental variables were used as predictors. The environmental distance between samples was directly calculated by the gdm package, but for plant composition, we used our own dissimilarity indices to represent changes in vegetation composition from both a taxonomic and functional point of view (Table 1). Plant taxonomic dissimilarity between plots was calculated using thebeta.pair function from the R package ‘betapart’ (Baselgaet al. 2022). To estimate plant functional dissimilarity we retrieved trait values for species of which cumulative abundances represented at least 90% of the plot coverage. Missing values (<15% of the total data) were estimated using the imputation method offered by the R package ‘mice’ (van Buuren & Groothuis-Oudshoorn 2011). Plant functional dissimilarity was estimated using the beta.fd.multidim function from the R package ‘mFD’ (Magneville et al. 2021). In the analysis, soil samples from the same subplot collected in different years (i.e., 95 of the 418 samples) were treated as separate samples as they differed in environmental conditions, except for climate and NDVI that were averaged across a period of 10 years. Their spatial dependency was considered through the spatial coordinates of the plot.