Field sites and sampling
The study was based on extensive data collection and analysis across 343 study sites and 15 countries. About half of the data were published before (53.6%) 46–48,62,71–93, while other data were compiled for this study. The dataset comprised 15,893 sample records of paired δ13C and δ15N values in soil animals distributed across four climatic regions: subarctic, temperate, subtropical, and tropical. The investigated ecosystem types included woodlands, agricultural systems and grasslands. Mean annual precipitation and temperature were included as abiotic drivers and were retrieved from WorldClim 94 based on latitude and longitude. The details of study sites are listed in Table S8.
For details on the sampling methods for published data see 46–48,62,71–93. For unpublished data, standard extraction methods were used. Nematodes were sampled by extracting 5 cm diameter soil cores encompassing the litter layer and the top 0–5 cm of the mineral soil. Soil meso- and macrofauna were sampled by using heat Berlese or Kempson extractors 95 and preserved in 70-96% ethanol. Sampling methods deviations are listed directly in the Table S8.
Animals were classified into 26 high-rank taxonomic groups and further into five major functional groups as follows: herbivores (Hemiptera, Orthoptera, Thysanoptera, and Lepidoptera), detritivores (Lumbricina, Diplopoda, Isopoda, Isoptera, Dermaptera, Blattodea, Gastropoda, and Enchytraeidae), microbivores (Collembola, Oribatida, Nematoda, Protura, Prostigmata, Psocoptera, and Symphyla) and predators (Araneae, Chilopoda, Diplura, Formicidae, Mesotigmata, Opiliones, and Pseudoscorpiones) and groups displaying mixed feeding (Diptera and Coleoptera) 15,96. It has been shown that high-rank animal taxa in soil typically are trophically and functionally consistent 96.
Stable isotope analysis
Animals were identified to family- (86.9%), genus- (70.1%) or species- level (58.5%) before being processed for stable isotope analysis. Before stable isotope analysis, animals were dried at 50-60°C for 24 h, then weighed and enclosed in tin capsules; sample weights ranged from 0.01 to 1.0 mg. For small-sized animals, the whole body of individual animals were used for stable isotope analysis, with multiple individuals were bulked when more biomass was required, for large-sized animals we used body parts dominated by muscle tissue (e.g., legs) 97. Stable isotope ratios of 13C/12C and 15N/14N were measured using a system comprising an elemental analyser and a mass spectrometer; for details of the individual set ups see Table S8. Ratios between the heavy isotope and the light isotope (13C/12C, 15N/14N; R) were presented in parts per thousand relative to the standard using the delta notation, denoted as δ13C or δ15N = (Rsample/Rstandard − 1) × 1000 (‰). Vienna PD Belemnite and atmospheric nitrogen served as the standards for 13C and 15N, respectively.
Calculation of trophic diversity and dissimilarity
Trophic diversity of soil animals can be determined by computing the standard ellipse area (SEA) based on the position of soil animals within the δ13C-δ15N biplot of taxonomic and functional groups at each site. We used corrected standard ellipse area (SEAc) instead of SEA in our study, which is more robust in handling small and variable sample sizes than SEA 40. The relationship between SEA and SEAc can be formulated as SEAc = SEA(nsample size -1)(nsample size -2)-1. We visualise some examples of the SEAc of detritivores, microbivores and predators in woodland the agricultural systems by randomly picked 5 sites for each group (Fig. S8). Moreover, to further mitigate potential bias stemming from small sample sizes, ellipses were exclusively computed for taxonomic and functional groups that consisted of five or more samples per site. To assess trophic dissimilarity among taxonomic groups within each functional group, we calculated the mean pairwise distance between the centroids of isotopic-positions of taxonomic groups within each functional group.
We used uncalibrated stable isotope values (δ13C and δ15N) for assessing trophic diversity and trophic dissimilarity of soil animals, as calibration using litter δ13C and δ15N values did not significantly affect SEAc and trophic dissimilarity. We calculated the standard deviation of δ13C and δ15N values within the same functional group at each site. These values served as indicators of the variation in both basal resource use and trophic position among functional groups.
Statistical analyses
All analyses were done in R 4.0.3 98. To assess the effects of land use and climate on the SEAc and trophic dissimilarity of soil animals, we selected the subsets from tropical/temperate and agricultural systems/ woodlands from the whole dataset. We fitted linear mixed-effects models (LMMs) using log-transformed SEAc and trophic dissimilarity, and then applied contrasts between tropical and temperate ecosystems, as well as between agricultural systems and woodlands to estimate effect sizes. We conducted three separate LMMs for log-transformed SEAc of functional groups, SEAc of taxonomic groups and trophic dissimilarity. The models included functional groups (herbivores, detritivores, microbivores and predators), land use (agricultural systems, woodlands), biome (tropical, temperate) and their interactions as fixed effects, with site included as random effect.
For estimating effect sizes of land use and climate, we used the emmeans package to compute the estimated marginal means in the linear models. Then, we used the contrast function from the emmeans package to calculate contrasts between temperate versus tropical and woodland versus agriculture 99.
Additionally, we built another model and included sampling number and family richness as covariates to inspect their effects on the log-transformed SEAc of functional groups. The models included functional groups, land use, climate, sampling number and family richness as fixed effects, with site included as random effect. We checked model assumptions of the most parsimonious models by fitting model residuals versus the results of fitted models.
To elucidate the drivers behind larger functional group SEAc, we employed two separate LMMs. One model included SEAc of taxonomic groups and trophic dissimilarity as explanatory variables, while the other explored variations in δ13C and δ15N as separate explanatory variables. We then used the estimated value of the coefficient for each independent variable to estimate their contribution to SEAc of functional groups.
For exploring whether temperature or precipitation are responsible for the observed differences in trophic diversity among the four latitudinal zones, we used structural equation modelling (SEM) employing maximum likelihood estimation to provide insight into how environmental factors affected trophic diversity of soil animals. The SEM included MTP and MAP as direct climatic effects. Furthermore, we included SEAc of taxonomic groups and trophic dissimilarity as mediators affected by environmental factors. The data was z-score scaled before SEM analysis. To determine the goodness of fit of the model, we used χ2-test, the comparative fit index (CFI) > 0.95 and the standardised root-mean-square residual (SRMR) with values ≤ 0.05 as threshold for significance 100. Our SEM adequately described the data (p = 0.33, df = 3, CFI = 0.999, RMSEA = 0.022, SRMR = 0.007, AIC = 1724.75). We also built a competing model that allows for direct climate effects on SEAc, but the model had higher AIC (1725.32), and both MTP and MAP had no significant effects on SEAc (Fig. S9), thus we used the previous model in the main text.
The SIBER package was used for calculating the trophic diversity of soil animals 40. The lme4 package was used to fit LMMs 101 and the emmeans package to compute the estimated marginal means in the linear models 99. The SEM analysis was performed with the lavaan package 102. All mixed models were visually checked to meet the assumption of residual homogeneity of variance. Results were visualised using the ggplot2 package 103.