Warming alters the energetic structure and function but not resilience of soil food webs.

Climate warming is predicted to alter the structure, stability, and functioning of food webs1-5. Yet, despite the importance of soil food webs for energy and nutrient turnover in terrestrial ecosystems, warming effects on these food webs-particularly in combination with other global change drivers-are largely unknown. Here, we present results from two complementary field experiments testing the interactive effects of warming with forest canopy disturbance and drought on energy fluxes in boreal-temperate ecotonal forest soil food webs. The first experiment applied a simultaneous above- and belowground warming treatment (ambient, +1.7°C, +3.4°C) to closed canopy and recently clear-cut forest, simulating common forest disturbance6. The second experiment crossed warming with a summer drought treatment (-40% rainfall) in the clear-cut habitats. We show that warming reduces energy fluxes to microbes, while forest canopy disturbance and drought facilitates warming-induced increases in energy flux to higher trophic levels and exacerbates reductions in energy flux to microbes, respectively. Contrary to expectations, we find no change in whole-network resilience to perturbations, but significant losses of ecosystem functioning. Warming thus interacts with forest disturbance and drought, shaping the energetic structure of soil food webs and threatening the provisioning of multiple ecosystem functions in boreal-temperate ecotonal forests.

Climate warming modifies consumer-resource interactions5,7,8, consequently altering food web structure and fluxes of energy and matter through ecosystems1–4. The energetic organization of food webs affects their stability9–11 and the provisioning of ecosystem functions12, as energy fluxes among trophic groups describe the ecosystem functions carried out by each consumer. In accordance with the temperature dependence of metabolism13,14, energy fluxes through food webs should increase with warming as energetic demands increase. It is also likely that variation in responses of different trophic groups to warming2, e.g., due to taxon-specific metabolic scaling with temperature14 or warming-induced mismatches of energy demand and consumption7,15, could reorganize the biomass structure of food webs1,4. Thus, warming-induced declines in biomass of certain trophic groups could decrease absolute energy fluxes to these groups, although the energy flux per unit biomass may increase with warming. Other global change drivers acting in concert with climate warming may further influence energy fluxes. For example, recent studies have shown that drought and land-use intensification reduce energy flux through food webs12,16,17. Accounting for possible interactions with other global change drivers may therefore advance our understanding of warming effects on food web structure and the provisioning of multiple ecosystem functions.

In terrestrial ecosystems, most of the energy supplied by primary production directly enters the soil and fuels belowground food webs18. These food webs provide key ecosystem functions described by the flux of energy among trophic levels19,20; for example the mineralization of nutrients and the sequestration of carbon, which fundamentally affect all terrestrial organisms and regulate soil carbon-climate feedbacks21. Despite the inordinate importance of soil food webs, our knowledge is still limited on how climate warming and its interactions with other global change drivers impacts structure, stability and functioning of whole soil food webs. Warming and drought in microcosms were reported to increase consumption rates of predators, thereby reducing microbial-feeding prey and increasing fungal-driven decomposition22. Other studies have found that warming restructures the community composition of soil microbes23, soil microarthropods24, and nematodes25, which potentially impact the provisioning of ecosystem functions and the stability of these communities as determined by their energetic structure10. Yet, longer-term warming experiments have shown such responses to vary over time, likely due to adaptation or compositional shifts towards better adapted species23,26. Furthermore, other global change drivers, such as drought and land use16 were shown to shape the influence of climate change on soil food webs.

Here, we test how warming alters energy fluxes through soil food webs in boreal-temperate ecotonal forest and how these effects interact with forest canopy disturbance and summer drought. We hypothesized that energy fluxes in soil food webs would increase with warming due to the positive temperature-metabolism relationship13,14, but that this increase would be strongly counteracted in treatments crossed with forest disturbance and drought due to losses of species richness and biomass12,17. Additionally, we expected that warming-induced increases in energy fluxes would yield higher interaction strengths at different trophic levels resulting in reduced whole-network resilience to perturbations, with drought and forest disturbance further exacerbating lowered network resilience. Taken together, we expected that the overall energetic structure of the soil food webs would be affected by the combined global change treatments, resulting in potential decreases of ecosystem functions12,20 and the resilience of whole networks to further perturbations.

To test our hypotheses, we performed two complementary field experiments, each replicated at two sites in northern Minnesota, USA6. The first experiment applied a simultaneous above- and belowground warming treatment (ambient, +1.7°C, +3.4°C), both in 40- to 60-year-old forest stands and in recently clear-cut habitats with planted tree seedlings of eleven species (see Methods). We sampled soil in both habitat types after two years of warming to investigate the effect of warming in combination with forest canopy disturbance. The second experiment crossed the warming treatment in the clear-cut habitats with a summer drought treatment, and was sampled after seven years of warming and four years of summer drought.

Soil food webs were assessed by measuring free-living soil nematodes, microarthropods, microbial biomass carbon (C), and basal respiration. Based on literature, we aggregated nematode and microarthropod taxa into feeding guilds and assembled these as well as microbes into a functional group-level food web20,27. We calculated community metabolism for each trophic group using metabolic rates for soil animals derived from scaling relationships of body mass, soil temperature, and phylogeny14 along with measured microbial basal respiration. Accounting for taxa-specific assimilation efficiencies from the literature20,27,28 and food web structure, we calculated the energy flux necessary to support the energetic demands of each feeding guild12 (see Methods for details). These energy fluxes describe the resource uptake of each trophic group in the food web and, consequently, represent functions carried out by these groups. For example, the energy flux to herbivores serves as a measure of herbivory, while the energy fluxes to microbes and detritivores should be directly related to total decomposition rates in the soil20. In order to test this theoretical assumption, we correlated the calculated energy fluxes to microbes and detritivores with data on measured decomposition of a standardized cellulose-based substrate in the field29 (see Methods). We found a significant positive relationship between the combined energy fluxes to microbes and detritivores and the removal rate of cellulose (r2 = 0.16, p-value = 0.009; Supplementary Figure 8), supporting the claim that energy flux to specific trophic groups describes their respective ecosystem functions.

We summed up energy fluxes within the food web to obtain a community-level measure of energy flux as well as per trophic group (microbes, detritivores, herbivores, predators). We calculated the relative contribution of each of these groups to the total community energy flux and also differentiated microbes from total fauna. To analyze the stability of food webs, we calculated interaction strengths between the feeding guilds based on energy fluxes and constructed Jacobian matrices. Stability was assessed as the resilience of the whole network to a small perturbation by determining the minimum intraspecific interaction strength (i.e. the diagonal values of the Jacobian matrix) needed such that all eigenvalues (calculated from the matrix of interaction strengths) have negative real parts30 (see Methods for details).

Warming altered the amount and distribution of energy fluxes in the soil food web in both experiments, with varying responses to warming among different consumer groups (Fig. 1, Supplementary Figures 3 and 6). As expected, however, forest canopy disturbance and summer drought were important in determining the strength of warming effects on the energetic structure of soil food webs. In the first experiment, warming reduced the energy flux to microbes in both undisturbed and disturbed canopy habitats (Fig. 2a, Table 1), contrary to our first hypothesis. In disturbed canopy habitats, microbial mass-specific respiration even declined in response to warming (Supplementary Figure 4, Supplementary Table 1). This may be related to warming-induced decreases in soil moisture (Supplementary Figure 1), which is known to strongly limit the microbial community16 (Supplementary Figure 5). Faunal biomass and energy fluxes were generally greater in disturbed than in undisturbed canopy habitats (Fig. 1a, Fig. 2a, Table 1, Supplementary Figure 4), whereby a marginally significant interaction of warming and canopy disturbance suggests that warming increased faunal fluxes in disturbed canopy habitats (Fig. 2a, Table 1). Similarly, warming increased fluxes to detritivores, especially in disturbed canopy habitats (Fig. 1a, Table 1), indicating differences between both habitats in their faunal communities and how these respond to warming.

Figure 1

Warming effects on soil food web structure.

(a) Effects of warming and forest canopy disturbance on energy flux (arrow width) and biomass C (circle area) distribution among microbes (blue), herbivores (green), detritivores (brown), and predators (red) after two years of warming (first experiment, n=69). Microbial energy flux and biomass are displayed with only 10% and 1% of their actual size. Within treatment combinations, white baselines indicate circle and arrow sizes in the respective ambient temperature treatment. (b) Effects of warming and summer drought on energy flux and biomass distribution in the soil food web after seven years of warming and four years of summer drought (second experiment, n=35).

Figure 2

Warming effects on energy fluxes in the soil food web.

(a) Effects of warming in undisturbed and disturbed canopy habitats on overall energy flux in the whole food web, microbial-driven fluxes, faunal-driven fluxes, and on the relative contributions of microbes and fauna to whole food web energy fluxes (n=69). (b) Effects of warming under ambient and reduced rainfall on overall energy flux in the whole food web, microbial-driven fluxes, faunal-driven fluxes, and on the relative contributions of microbes and fauna to whole food web energy fluxes (n=35). Box plots are based on median (horizontal line), first and third quartile (rectangle), 1.5 × interquartile range (whiskers), and outliers (isolated points).

Table 1

Results of linear mixed effects models testing the effects of warming, canopy disturbance, and their interaction on energy fluxes and biomass of the whole food web, microbes, total fauna, and the three faunal trophic groups (detritivores, predators, herbivores), on the relative contributions of microbes, detritivores, predators, and herbivores to whole food web energy flux as well as on the resilience of the whole food web in the warming × canopy disturbance experiment. Significance of fixed effects was obtained by Wald chi-square tests. Significant effects (p < 0.05) are reported in bold. For models with significant (p < 0.05) and marginally significant (p < 0.1) interaction terms we performed post-hoc Tukey´s HSD tests that are reported in Supplementary Table 2.

	Transformation	Warming			Canopy			Warming x Canopy
		χ2	df	p-value	χ2	df	p-value	χ2	df	p-value
Energy flux (g C m−2 d−1)
Whole food web	log10(x)	2.757	1	0.097	0.062	1	0.803	0.102	1	0.749
Microbes	log10(x)	5.020	1	0.025	0.426	1	0.514	0.734	1	0.392
Fauna	log10(x)	2.138	1	0.144	10.554	1	0.001	3.380	1	0.066
Detritivores	log10(x)	4.230	1	0.040	12.267	1	<0.001	3.769	1	0.052
Predators	log10(x)	1.658	1	0.198	5.152	1	0.023	4.995	1	0.025
Herbivores	log10(x+0.01)	0.790	1	0.374	1.125	1	0.289	0.219	1	0.640
Biomass (g C m−2)
Whole food web	log10(x)	4.564	1	0.033	1.796	1	0.180	0.015	1	0.901
Microbes	log10(x)	3.372	1	0.066	1.913	1	0.167	0.027	1	0.871
Fauna	log10(x)	0.360	1	0.548	24.650	1	<0.001	1.441	1	0.230
Detritivores	log10(x)	2.009	1	0.156	24.515	1	<0.001	0.184	1	0.668
Predators	log10(x)	0.157	1	0.692	6.689	1	0.010	3.433	1	0.064
Herbivores	log10(x+1)	3.811	1	0.051	2.023	1	0.160	0.106	1	0.745
Relative energy flux (proportion)
Microbes	logit(x)	6.726	1	0.010	8.817	1	0.003	5.230	1	0.022
Detritivores	logit(x)	10.249	1	0.001	12.918	1	<0.001	5.665	1	0.017
Predators	logit(x)	4.582	1	0.032	4.485	1	0.034	6.594	1	0.010
Herbivores	logit(x+0.001)	1.615	1	0.204	0.756	1	0.385	0.279	1	0.597
Resilience to perturbations
Whole food web		0.816	1	0.366	3.082	1	0.079	0.369	1	0.544

Warming also had contrasting effects on energy fluxes to predators in both habitat types, with a reduction by 30% in undisturbed canopy habitats and an increase by 110% in disturbed canopy habitats (Fig. 1a, Table 1). Predators may not completely meet their increasing energetic demands under warming, as metabolic rates increase faster with temperature than their feeding rates7,15. In disturbed canopy habitats, the absence of negative warming effects on predators indicates that predators could meet their increasing energetic demands under warming, which likely was facilitated by the considerably higher prey densities in these habitats (Fig. 1, Table 1). These results suggest that warming could have contrasting effects on food web stability among habitat types. In undisturbed canopy habitats, consumer starvation may reduce top-down pressure and thus stabilize population dynamics7. In contrast, the increase in energy flux at higher compared to lower trophic levels as found in the disturbed canopy habitats is predicted to reduce food web stability9,11. Surprisingly, however, our analyses revealed that in both habitat types warming-induced alterations in the energetic structure of soil food webs had no effect on their resilience to perturbations (Fig. 3a, Table 1). This is probably due to simultaneous shifts in fluxes to different trophic groups and opposing trends in consumer densities compared to their metabolic rates. That is, self-damping due to increased metabolic demands in consumers can minimize destabilizing effects resulting from greater energy fluxes to higher trophic levels.

Figure 3

Effects of warming, canopy disturbance, and drought on whole network resilience.

Presented are treatment effects in (a) the first (n=69) and (b) the second experiment (n=34) on the degree of self-damping (minimum s) that is required for all eigenvalues of the Jacobian matrix to have negative real parts. Food webs that required a smaller degree of self-damping were considered to be more resilient to perturbations. Here, differences among treatments are not significant. Box plots are based on median (horizontal line), first and third quartile (rectangle), 1.5 × interquartile range (whiskers), and outliers (isolated points).

While warming appeared to decrease energy fluxes through the whole soil food web by 12% in both disturbed and undisturbed forest canopy habitats, this trend was not significant (Fig. 2a, Table 1). In undisturbed canopy sites, warming also did not affect the relative contribution of the four trophic groups to whole food web energy fluxes (Fig. 2a, Table 1), indicating that in habitats with intact canopy, warming might have relatively modest effects on ecosystem functions driven by the soil food web. In contrast, warming in disturbed canopy habitats shifted the relative contribution of all trophic groups except herbivores, with a decrease in the contribution of microbes (from 93% to 84%), but an overall increase in the contribution of fauna (from 7% to 16%) to whole food web energy flux (Fig. 2a, Table 1).

In the second experiment, summer drought clearly altered the response of microbes to warming. In plots with ambient rainfall, warming had no effect on microbial biomass and energy fluxes, while in plots subjected to drought both measures were severely reduced by warming (Figs. 1b and 2b, Table 2). The emerging negative warming effect on microbial mass-specific respiration in the first experiment was considerably pronounced after seven years of warming, with drought further reducing microbial mass-specific respiration (Supplementary Figure 7, Supplementary Table 4). Thus, warming likely affected the microbial community indirectly via decreasing soil moisture, which was negatively impacted by warming and drought (Supplementary Figure 1).

Table 2

Results of linear mixed effects models testing the effects of warming, drought, and their interaction on energy fluxes and biomass of the whole food web, microbes, total fauna, and the three faunal trophic groups (detritivores, predators, herbivores), on the relative contributions of microbes, detritivores, predators, and herbivores to whole food web energy flux as well as on the resilience of the whole food web in the warming × drought experiment. Significance of fixed effects was obtained by Wald chi-square tests. Significant effects (p < 0.05) are reported in bold. For models with significant (p < 0.05) and marginally significant (p < 0.1) interaction terms we performed post-hoc Tuckey´s HSD tests that are reported in Supplementary Table 6.

	Transformation	Warming			Drought			Warming x Drought
		χ2	df	p-value	χ2	df	p-value	χ2	df	p-value
Energy flux (g C m−2 d−1)
Whole food web	log10(x)	5.622	1	0.018	15.438	1	<0.001	7.892	1	0.005
Microbes	log10(x)	12.400	1	<0.001	11.927	1	0.001	8.455	1	0.004
Fauna	log10(x)	9.054	1	0.003	1.016	1	0.313	0.593	1	0.441
Detritivores	log10(x)	8.935	1	0.003	0.111	1	0.739	0.055	1	0.814
Predators	log10(x)	1.249	1	0.264	1.208	1	0.272	0.038	1	0.846
Herbivores	log10(x)	2.814	1	0.093	0.162	1	0.688	1.355	1	0.244
Biomass (g C m−2)
Whole food web	log10(x)	0.811	1	0.368	0.164	1	0.686	6.048	1	0.014
Microbes	log10(x)	0.872	1	0.351	0.164	1	0.685	6.088	1	0.014
Fauna	log10(x)	1.421	1	0.233	0.191	1	0.662	0.200	1	0.654
Detritivores	log10(x)	0.001	1	0.979	1.039	1	0.308	0.583	1	0.445
Predators	log10(x)	0.059	1	0.808	0.000	1	0.984	0.032	1	0.857
Herbivores	log10(x)	0.762	1	0.383	0.054	1	0.817	6.177	1	0.013
Relative energy flux (proportion)
Microbes	logit(x)	11.359	1	0.001	0.277	1	0.599	1.104	1	0.293
Detritivores	logit(x)	15.457	1	<0.001	0.676	1	0.411	1.147	1	0.284
Predators	logit(x)	1.074	1	0.300	1.025	1	0.311	0.013	1	0.910
Herbivores	logit(x)	2.545	1	0.111	0.095	1	0.758	1.297	1	0.255
Resilience to perturbations
Whole food web		0.098	1	0.754	0.685	1	0.408	1.299	1	0.254

In line with our hypothesis, warming significantly increased faunal-driven energy fluxes, especially to detritivores (Figs. 1b and 2b, Table 2), indicating increased decomposition carried out by soil detritivore fauna after seven years of warming. Drought, however, had no effect on either faunal biomasses or faunal energy fluxes (Fig. 1, Table 2). Whole food web energy flux was also not affected by warming in plots with ambient rainfall, but reduced by warming under summer drought (Fig. 2b, Table 2). Warming-induced shifts in the relative contribution of microbes and fauna to whole food web energy fluxes remained true after seven years of warming and under drought conditions (Fig. 2b, Table 2), confirming results of the first experiment. In disturbed canopy habitats, warming thus increased the energy transfer from microbes to higher trophic levels in both experiments (Fig. 1, Tables 1, 2), suggesting that detritivores became more efficient in consuming microbes at higher temperatures. This notion is supported by a soil microcosm experiment, which revealed that feeding rates of fungal feeders increase with warming and thereby reduce fungal biomass22. Thus, it is likely that elevated feeding rates of detritivores in warmed plots had negative effects on microbial biomass with potentially negative consequences for decomposition processes. As in the first experiment, the energetic restructuring of the soil food webs appeared to have surprisingly little influence on whole-network resilience to perturbations, even under combined warming and drought treatments (Fig. 3b, Table 2).

Our results reveal that energy fluxes at the whole food web level were only slightly affected by warming, yet varying effects on different trophic groups clearly modified the interplay of different ecosystem functions, such as herbivory, decomposition, and predation. These soil functions are crucial for growth and regeneration of forests and the carbon dynamics of these ecosystems. While our results suggest that the structure and functioning of soil food webs in intact boreal-temperate ecotonal forests may be relatively well buffered against climate warming, canopy disturbance is likely to increase the susceptibility of these food webs to warming, which may be even further amplified under drought conditions. Surprisingly, the food webs in warmed soils exhibited a similar resilience as those of colder soils, even when subjected to forest disturbance and drought. Despite the restructuring of these food webs in response to multiple global change drivers, we suspect that compensatory dynamics of different functional groups may have dampened global change effects on food web resilience at the whole-network level. Our study opens new questions about whether the observed shifts in energy fluxes toward higher trophic levels might affect the temporal stability of specific functional groups. In this vein, future investigation of the temporal stability of specific consumer nodes in food webs subjected to similar global change scenarios are needed to shed further light on the resilience of ecological networks and their ecosystem functions in boreal forest soil food webs.

Methods

Field site and experimental design

The B4WarmED experiment (Boreal Forest Warming at an Ecotone in Danger) is located at the ecotone between boreal and temperate forests in northern Minnesota, USA, at two field sites: the Cloquet Forestry Center, Cloquet (46° 40′ 46″ N, 92° 31′ 12″ W, 382 m a.s.l., 4.5°C mean annual air temperature (MAT), 807 mm mean annual precipitation) and ∼150 km further north, the Hubachek Wilderness Research Center, Ely (47° 56′ 46″ N, 91°45′ 29″ W, 415 m a.s.l., 3.0°C MAT, 722 mm mean annual precipitation)6. The soil at Cloquet is a Typic Dystrudept and at Ely it is a Lithic Udorthent31, whereas the soil texture at both sites is classified as sandy loam.

Since 2009, a three-level warming treatment (ambient, +1.7°C, +3.4°C) was applied in 40–60 year old mixed aspen-birch-fir forest habitats with intact canopy (∼5–10% of full sunlight) and in disturbed canopy habitats (recently clear-cut, ∼80% of full sunlight) at both sites. The initial experimental design was a 2 (site) × 2 (canopy disturbance) × 3 (warming treatment) factorial, with six replicates (two per block), resulting in a total of 72 circular 3 m-diameter plots. In 2012, rain shelters were added to half of the plots in disturbed canopy habitats to simulate summer drought (~40% less rainfall in June-September), resulting in a 2 (site) × 2 (drought treatment) × 3 (warming treatment) factorial design with 3 replicates (1 per block).

Within each plot, 121 seedlings of eleven tree species (Abies balsamea, Acer rubrum, A. saccharum, Betula papyrifera, Picea glauca, Pinus banksiana, Pinus strobus, Populus tremuloides, Quercus macrocarpa, Quercus rubra, and Rhamnus cathartica) were planted into the occurring understory plant vegetation in 2008. In 2012, all disturbed canopy plots were harvested and re-planted according to the original planting scheme. Above- and belowground warming was accomplished for on average 8 months per year through the use of 6–8 ceramic heating elements (Mor-Electric, model FTE-1000, MI, USA) and heating cables (Danfoss GX, Devi A/B, Denmark) buried at 10 cm depth6. Simulated drought was accomplished by manually stretching rain-shelters over the plots in advance of rain events, which effectively reduced rainfall by about 40% in each year. Air temperature and soil temperature at 10 cm depth as well as soil volumetric water content were constantly logged6.

Soil sampling

Soil for extracting nematodes and analyzing soil microbial biomass carbon (C) and basal respiration was sampled by taking three soil cores (2 cm in diameter, 7 cm deep) per plot in August 2010 and four soil cores (2 cm in diameter, 7 cm deep) per plot in August 2015. The soil was pooled at the plot-level, homogenized, and stored in plastic bags at 4°C until further processing. The soil was sieved through a 2 mm mesh and subdivided for measurements of microbial biomass C and basal respiration, and extraction of nematodes. The subsample for microbial measurements was frozen until further processing and the subsample for nematode extraction was immediately processed. For soil mesofauna extraction, one soil core (7 cm in diameter, 7 cm deep in August 2010; 4.8 cm in diameter, 7 cm deep in August 2015) was taken from each experimental plot.

Soil bulk density was measured separately for 0 to 5 cm depth and for 5 to 10 cm depth at three locations in each block in 2008 (outside of treatment plots). Three cores (5 cm in diameter) per location were sampled and the dry mass of each was determined. The volume of the soil removed was determined by filling a flexible plastic bag placed in the soil hole with a known volume of water. From this, we calculated a mean bulk density of the upper 7 cm of soil for each block.

Microbial measurements

Soil microbial basal respiration and biomass C were determined by using an O2-microcompensation apparatus32. Microbial respiratory response was measured at hourly intervals for 24 h. Basal respiration (µl O2 h−1 g−1 soil dry weight) was determined without addition of substrate and was measured as mean of the O2-consumption rates of hours 14 to 24 after the start of the measurements. Substrate-induced respiration was determined from the respiratory response to D-glucose33, whereas the mean of the lowest three readings within the first 10 h was taken as maximum initial respiratory response (MIRR; µl O2 h−1 g−1 soil dw). From this, microbial biomass (µg C g−1 soil dw) was calculated as 38 × MIRR (Ref. 34). O2-consumption rates with and without addition of substrate were measured at a standard lab temperature of 20°C in 2010. In 2015, we measured substrate-induced respiration at 20°C and basal respiration at the mean temperature of the respective warming treatment in the sampling month (17.2°C, 18.9°C, and 20.6°C, respectively). To standardize units between microbial and soil animal data, we used soil bulk density data to express soil microbial respiration and biomass as values per m2, allowing for direct comparisons between these groups.

Soil nematodes

Soil nematodes were extracted from approximately 10 g of fresh soil in 2010 and 25 g of fresh soil in 2015 using a modified Baermann method35. Nematodes were preserved in 4% formaldehyde, and individuals were counted and identified to family (juveniles) or genus (adults and most of the juveniles) level wherever possible. In 2010, ca. 46% of nematode identifications could only be done to the subclass level. At least 100 individuals (if available in the sample) were randomly identified following Bongers36. Nematode taxa were then assigned to trophic groups as bacterial feeding, fungal feeding, plant feeding, predators, and omnivores37,38. We assessed taxon-specific data on nematode fresh body mass from the “Nemaplex” database (http://plpnemweb.ucdavis.edu/nemaplex/Ecology/nematode_weights.htm, Supplementary Table 6). Mean fresh body masses for subclasses were calculated as the mean of the fresh body masses of families within their respective subclasses that were also present in samples collected in 2015. Nematode numbers were related to g soil dry weight and subsequently transformed to densities per m2 by using soil bulk density data. We converted mean nematode fresh masses into mean carbon dry masses assuming a dry matter content of 25% (Ref. 39) and a carbon content of 50% of the dry matter40. Based on densities per m2 and mean carbon dry masses, we estimated the biomass of the five nematode feeding guilds per m2.

Soil mesofauna

Soil mesofauna was extracted from intact soil cores by gradually heating the soil cores for seven days from 25°C up to 50°C (Ref. 41). Extracted animals were preserved in ethanol (70%), counted (abundance per m2), and assigned to a taxonomic group following Schäfer and Brohmer (2006) (Ref. 42) and Crotty and Shepherd (2014) (Ref. 43). Animals were identified mainly to the level of orders and families (see Supplementary Table 7 for list of identified taxa).

Sampling and extraction methods were targeted on assessing soil mesofauna44. We therefore excluded taxa from analyses for which sampling and extraction methods were inadequate or which were not tightly connected to the soil food web (e.g., Annelida, Chilopoda, Coleoptera, Diplopoda, adult Diptera, Gastropoda, Hymenoptera, and Isopoda). Body lengths of mesofauna were determined to the nearest 0.01 mm using an ocular micrometer mounted in a dissecting microscope. We measured the body lengths of all individuals sampled in 2015 (Supplementary Table 7b). However, it was not possible to directly assess the body lengths of soil mesofauna sampled in 2010. Thus, we assessed a mean body length per taxon (Supplementary Table 7a) based on animals which were collected in 2012 and 2015 from the same experiment.

We used length-mass regressions from the literature (Supplementary Table 8) to convert body lengths into body masses. We calculated the biomass per taxonomic group and biomass of the total community at the plot-level by multiplying mean body masses by densities (2010) and by summing together individual body masses (2015). We expressed biomasses as carbon dry masses per m2. For most taxonomic groups, we assessed dry masses using length-mass regressions. Dry mass of Acari were assumed to be 43.1% of fresh mass45. Collembola fresh masses were converted into dry masses using the equation given by Mercer et al. (2001) (Ref. 46). Carbon content of dry mass of Acari and Collembola was set to be 48% (Ref. 47), and for other groups we assumed a dry mass carbon content of 50% (Ref 48).

Construction of the food web

Based on published feeding relationships in the soil food web20,27,40, micro- and mesofauna groups were lumped into nine feeding guilds: herbivorous nematodes, bacterivorous nematodes, fungivorous nematodes, predatory nematodes, herbivorous mesofauna, detrivorous mesofauna, Oribatida, predatory mites, and predatory mesofauna, from which the food web was constructed (Supplementary Figure 2).

Nematode taxa were classified as either herbivorous, bacterivorous, fungivorous, or predatory37,38. Predatory nematodes were assumed to feed on all other nematode groups37,40. Taxa typically designated as omnivorous nematodes, which mainly feed on bacteria, fungi, and microfauna37,49,50, were proportionately assigned to bacterivorous, fungivorous, and predatory nematodes according to the assumption that their diet consists of 25% bacteria, 25% fungi, and 50% other nematodes.

Aphidina, Ciccadina, Coccina, Heteroptera nymphs, and Thysanoptera were classified as herbivorous mesofauna. Detritivorous mesofauna comprised microbivorous mites, all Collembola, which generally are treated as microbial (fungal) feeders40,49, as well as Diptera larvae, Pauropoda, and Protura51. Adult Oribatida, although being mainly microbial feeders, were treated as a separate feeding guild, as they are relatively resistant to predation due to their strong sclerotization40,52. It was assumed that only 25% of adult Oribatida are subjected to predation. All mesostigmatid mites (mainly Gamasina), as well as Bdellidae, Cunaxidae, Rhagidiidae, and 50% of other Eupodoidea were assumed to be predators feeding on all nematode groups, microbivorous mesofauna, and Oribatida53–55. Astigmatic and prostigmatic mites not designated as predators were assumed to be microbivorous40,54. For the first experiment (2010), astigmatic and prostigmatic mites were not further identified. As particular prostigmatic mites contain many predators, it was assumed that these groups consist of 75% microbivores and 25% predators following Andrén et al. (1990) (Ref. 54). Araneae, Pseudoscorpionida, and Diplura were assumed to prey on herbivorous mesofauna, detrivorous mesofauna, Oribatida, and predatory mites56 and were thus classified as predatory mesofauna. Symphyla are omnivorous, feeding on a variety of resources51, but mainly on animal prey55. Assuming that their diet consists of 20% plant material, 20% microbes, and 60% animal prey, they were proportionately assigned to herbivorous mesofauna, detritivorous mesofauna, and predatory mesofauna, respectively.

Calculation of metabolic rates

We calculated mean metabolic rates for each identified animal taxon in 2010, while we calculated mean metabolic rates for each identified nematode taxon and individual metabolic rates for all arthropod individuals sampled in 2015. The calculations were based on fresh body masses, temperature, and, if possible, phylogeny14 using the formula where I is the metabolic rate, a is the allometric exponent, M is the fresh body mass, E is the activation energy, k is the Boltzmann constant, T is the temperature, and i is a normalization factor12,14. For i, a, and E, taxon-specific values were used, if available; otherwise, general parameters12,14 were used (Supplementary Table 10). This yielded metabolic rates expressed as J h−1. Using the conversion factor 1 ml O2 = 20.1 J (Ref. 57), we converted microbial basal respiration from µl O2 h−1 to J h−1. Subsequently, all metabolic rates were expressed as W m−2.

Calculation of energy fluxes

For energy flux calculations, we assumed a steady state, which means that the energy flux to a feeding guild in the food web exactly balances the energy losses of that feeding guild due to metabolism and predation. We calculated the energy flux to each feeding guild in the food web as where F is the total flux of energy into the feeding guild, e is the specific assimilation efficiency20,27,28, X is the summed metabolic rates of all individuals within the feeding guild, and L is the energy loss to predation12. To account for potential variations in microbial assimilation efficiency, we used microbial mass-specific respiration as an estimate of the inverse microbial efficiency58 and rescaled the values to vary between 0.35 to 0.45. The difference of these values from one was taken as microbial assimilation efficiency, which accordingly varied in the range from 0.55 to 0.65, as was suggested by Adu and Oades (1978) (Ref. 28).

We started by calculating the energy flux to the top predators in the food web, for which energy loss to predation was assumed to be zero. We proceeded downwards to the guilds with the lowest trophic position, where the loss to predation of these lower feeding guilds was synonymous with energy fluxes to their consumer feeding guilds. We assumed that polyphagous feeding guilds consumed their resources according to the relative numerical abundance of the resources. Predatory mites feed on microarthropods and nematodes, which differ substantially in their densities. We thus applied density-dependent feeding only within these groups and assumed that predatory mites consume these groups based on their relative biomass and a 2:1 preference of microarthropods over nematodes27. Symphyla were assumed to feed on all trophic levels. We assumed that they constantly derive 20% of their energy from plant material, 20% from microbes, and 60% from animal prey. We assumed that Symphyla consume mesofauna and nematodes depending on the relative biomass of these groups and that within these groups they consume prey guilds based on their density. For energy flux calculations, Symphyla were treated as an independent feeding guild. For further analyses, however, energy fluxes to Symphyla were lumped with fluxes to detritivorous, herbivorous, and predatory mesofauna according to the above-mentioned composition of their diet. Adult Oribatida are relatively resistant to predation and, thus, only 25% of their density was used to weight the energy flux to predator guilds. We assumed all detritivorous and microbivorous feeding guilds to derive their energy from the microbial biomass pool.

Energy fluxes were summed up per trophic level (microbes, detritivores, predators, herbivores), and we differentiated microbes from total fauna. Finally, total energy flux was calculated by summing up all fluxes within the food web. In addition, the relative contributions of the different trophic levels, and of microbes and total fauna to total energy flux were calculated. Energy fluxes were expressed as g C m−2 d−1 using the conversion factor 1 g C = 46 kJ (Ref. 59).

Bait-lamina tests

During the second experiment in 2015, bait-lamina tests29,60,61 were conducted in all plots in disturbed canopy habitats to regularly (~every two weeks) assess the decomposition of a standardized cellulose-based substrate under field conditions. We used bait-lamina strips—16 cm long PVC strips with 16 holes of 2 mm in diameter—and filled their holes with a substrate consisting of 65% cellulose (micro granular), 15% agar (pulverized), 10% loess, and 10% wheat bran (finely ground and sieved). In each plot, six strips were vertically inserted into the soil (to a depth of ~10 cm; at a distance of >2 cm), and after two weeks, decomposer activity was assessed by counting the number of holes where the substrate was completely removed or perforated. Empty holes were scored as 1, whereas holes containing perforated substrate were scored as 0.5. The summed scores per strip were averaged on the plot-level for further analyses. To relate these data to energy fluxes based on soil food web data recorded in August 2015, we used the mean of two bait-lamina measurements taken in August 2015.

Assessment of network resilience

Based on energy fluxes between the feeding guilds, we calculated interaction strengths and constructed Jacobian matrices to analyze food web stability. The strength of the effect of resource i on consumer j was taken as where e is the assimilation efficiency of consumer j, F is the energy flux from resource i to consumer j, and M is the biomass of resource i. The strength of the effect of consumer j on resource i was taken as where F is the energy flux from resource i to consumer j, and M is the biomass of consumer j. Diagonal elements of the Jacobian matrix corresponded to the strength of the effect of population i on population i (intraguild interaction strength) and was taken as where X is the metabolism of feeding guild i, M is the biomass of feeding guild i., and s is a free parameter (but identical for all feeding guilds i) between 0 and 1. Increasing s generally decreases the real parts of the eigenvalues of the Jacobian matrix, whereas the matrix is considered stable when all eigenvalues have negative real parts. The stability of a food web was assessed by determining the minimum intraspecific interaction strength (the minimum s) required for all eigenvalues of the Jacobian matrix to have negative real parts, whereby food webs that required less intraspecific interaction strength (smaller s) to achieve negative real parts were considered more stable30,62.

Statistical analyses

Across both experiments (in 2010 and 2015), we tested treatment effects on a range of response variables: energy flux, biomass, and metabolism for the whole food web, the microbes, and the total fauna, as well as separately for the three faunal trophic groups (detritivores, predators, herbivores); mass-specific metabolism of microbes, detritivores, predators, and herbivores; the relative contributions of these groups to whole food web energy fluxes; and the degree of self-damping as a measure of inverse network resilience. To test for treatment effects on all response variables, linear mixed effects models with Gaussian error and the restricted estimates maximum likelihood method (REML) were fitted with block nested in site as random factors, using the “nlme” package in R 3.1.3 (RDevelopmentCoreTeam, http://www.R-project.org). For the first experiment, warming (numeric variable), canopy disturbance (categorical variable), and their interaction were analyzed as fixed effects. In the same way, warming (numeric variable), summer drought (categorical variable), and their interaction were used as fixed effects in the second experiment. We obtained chi-square values using Type II Wald chi-square tests63 implemented in the “car” package in R 3.1.3. To meet the assumptions of normality and ensure homogeneity of variance, the logit-transformation implemented in the “car” package was used for response variables from proportion data, while all other response variables except the degree of self-damping were log10-transformed. To transform variables containing zeros, the constant 10− was added to the respective variables, where k was a positive integer that was chosen separately for each variable to best meet model assumptions. All response variables were standardized according to Gelman (2008) (Ref. 64) using the “arm” package, to be able to compare the effects of predictor variables on different response variables. We also standardized numeric predictor variables to reduce multicollinearity. To interpret models with significant interaction terms, post-hoc Tukey’s HSD tests were performed using the “multcomp” package65 in R 3.1.3.

In the first experiment, additional mixed effects models were fitted with soil moisture as a co-variable of the treatment variables to test whether microbial measures were affected by water availability. As effects of soil moisture may vary depending on different treatments, initial models also included all possible interactions between soil moisture and the treatment variables as fixed effects. Subsequently, best-fit models were selected based on AIC in a stepwise backward selection procedure eliminating interaction terms but not main effects. For the model selection procedure, models were fitted using maximum likelihood (ML) estimation. The best-fit model was tested using the restricted estimates maximum likelihood method (REML)66.

In order to test the theoretical assumption that energy fluxes describe ecosystem functions, we correlated the calculated energy fluxes to microbes and detritivores with the mean decomposition of a standardized cellulose-based substrate in August 2015. We used ranged major axis regression as implemented in the “lmodel2” package in R 3.1.3.