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.