Soil fauna diversity increases CO2 but suppresses N2 O emissions from soil.

Soil faunal activity can be a major control of greenhouse gas (GHG) emissions from soil. Effects of single faunal species, genera or families have been investigated, but it is unknown how soil fauna diversity may influence emissions of both carbon dioxide (CO2 , end product of decomposition of organic matter) and nitrous oxide (N2 O, an intermediate product of N transformation processes, in particular denitrification). Here, we studied how CO2 and N2 O emissions are affected by species and species mixtures of up to eight species of detritivorous/fungivorous soil fauna from four different taxonomic groups (earthworms, potworms, mites, springtails) using a microcosm set-up. We found that higher species richness and increased functional dissimilarity of species mixtures led to increased faunal-induced CO2 emission (up to 10%), but decreased N2 O emission (up to 62%). Large ecosystem engineers such as earthworms were key drivers of both CO2 and N2 O emissions. Interestingly, increased biodiversity of other soil fauna in the presence of earthworms decreased faunal-induced N2 O emission despite enhanced C cycling. We conclude that higher soil fauna functional diversity enhanced the intensity of belowground processes, leading to more complete litter decomposition and increased CO2 emission, but concurrently also resulting in more complete denitrification and reduced N2 O emission. Our results suggest that increased soil fauna species diversity has the potential to mitigate emissions of N2 O from soil ecosystems. Given the loss of soil biodiversity in managed soils, our findings call for adoption of management practices that enhance soil biodiversity and stimulate a functionally diverse faunal community to reduce N2 O emissions from managed soils.

INTRODUCTION

There is considerable interest in how soils influence global climate change, and through which feedback mechanisms soil quality may be deteriorating. The biodiversity within an active soil biome is a large contributor to the soil function–climate feedback (Wall, 2012). However, empirical research on the connection between soil biodiversity and climate regulation is lacking. Specifically, the close relationship between soil‐derived greenhouse gas (GHG) emissions and soil processes such as biogeochemical cycling of C and N may result in feedback mechanisms that can either amplify or diminish the initial climate forcing (Crowther et al., 2016; Van Nes et al., 2015). As part of the soil biodiversity, soil fauna species can have a major impact on this relationship through their influence on biogeochemical cycling of nutrients (Coulis et al., 2015; Filser et al., 2016; Sauvadet, Chauvat, Brunet, & Bertrand, 2017; Wall et al., 2008), plant productivity (Liiri, Setala, Haimi, Pennanen, & Fritze, 2002; Van Groenigen et al., 2014), and soil‐derived GHG emissions (Lubbers, Van Groenigen, et al., 2013).

Agroecosystems are a known focal point of problems related to soil biodiversity, soil‐derived GHG emissions and the biogeochemical cycling of nutrients. Feedback mechanisms that may influence global climate change are especially relevant in these agroecosystems, where soil biodiversity is mostly negatively affected by agricultural intensification (Ponge et al., 2013; Postma‐Blaauw, Goede, Bloem, Faber, & Brussaard, 2012; Tsiafouli et al., 2015), where inputs of nitrogen (N) fertilizers are increasing globally (FAO, 2016), and where soil‐derived GHG emissions are high (IPCC, 2014). Agroecosystems also have the best potential for C storage in the soil by restoring soil organic C that was previously lost due to agricultural practices (Lal, 2004; Paustian, Cole, Sauerbeck, & Sampson, 1998; Smith, 2016; Stockmann et al., 2013). Yet, increasing soil C stocks can be negated by elevated soil‐derived emissions of the potent GHG N2O (Lubbers, Groenigen, Brussaard, & Groenigen, 2015; Powlson et al., 2014; Powlson, Whitmore, & Goulding, 2011). Global change models that underpin predictions of future climate scenarios still lack the potentially intricate and non‐linear interactions between soil fauna and microbes regulating soil biogeochemical cycles (Grandy, Wieder, Wickings, & Kyker‐Snowman, 2016). The lack of sufficient empirical data addressing ecological feedbacks to ecosystem processes such as GHG emissions precludes the parameterization necessary for further optimization of such models (Blagodatsky & Smith, 2012; Crowther et al., 2015; Grandy et al., 2016). To further unravel how soils and soil disturbance contribute to GHG emissions, mechanistic research is needed on functional diversity of soil fauna affecting the soil processes underlying climate regulation.

Belowground, gas fluxes of CO2, methane (CH4) and N2O are the result of a variety of (micro)biotic processes, occurring in ‘hotspots’ and during ‘hot moments’ (Kuzyakov & Blagodatskaya, 2015): CO2 is produced by soil respiration (including root, faunal and microbial respiration; Rastogi, Singh, & Pathak, 2002), CH4 through methanogenesis (Le Mer & Roger, 2001), while N2O is produced by a combination of microbial transformation processes, most importantly denitrification (Opdyke, Ostrom, & Ostrom, 2009), as well as nitrification and nitrifier–denitrification (Wrage‐Monnig et al., 2018). All these soil (micro)biotic processes resulting in the production of CO2, CH4 and N2O are controlled by numerous interacting conditions, characterized by key physico‐chemical factors such as the availability of labile C, mineral N and oxygen (O2), soil pH, soil moisture, temperature and diffusivity (Granli & Bøckman, 1994). Whereas CO2 is an end product of soil respiration, the net efflux of CH4 is the balance between anaerobic decomposition by methanogens and oxidation by methanotrophic bacteria (Dalal & Allen, 2008). Nitrous oxide can be both an intermediate and a by‐product of N transformation processes, and N2O gas can therefore also be consumed during (upward) diffusion through the soil profile (Klefoth, Oenema, & Willem Van Groenigen, 2012). This makes the net emission rates of N2O from soils highly variable and notoriously difficult to predict: N2O emission is dependent on the prevailing soil physico‐chemical factors, as well as on soil faunal activity (Kuiper, Deyn, Thakur, & Groenigen, 2013; Lubbers, Van Groenigen, et al., 2013).

The only soil faunal group so‐far acknowledged to have a pronounced influence on soil‐derived N2O emissions are earthworms (Lubbers, Van Groenigen, et al., 2013). These macro‐detritivores have been well studied with respect to N2O emissions and they can increase N2O emissions considerably through interacting with the microbial community (Drake & Horn, 2006; Nebert, Bloem, Lubbers, & Groenigen, 2011), by increasing the mineralization of nutrients and altering the soil structure (Lubbers, Brussaard, Otten, & Groenigen, 2011; Lubbers, Van Groenigen, et al., 2013). Other macrofauna, such as isopods and grubs, and mesofaunal groups, such as potworms, fungivorous mites and springtails, can also influence N2O emissions (Kuiper et al., 2013; Majeed et al., 2014; Tianxiang, Huixin, Tong, & Feng, 2008; Van Vliet, Beare, Coleman, & Hendrix, 2004). However, their (mostly enhancing) effects on N2O emission have mainly been determined as single‐species effects. The few studies on species combinations of different soil fauna feeding groups with respect to N2O emissions show diverging results: from a strong increase of soil‐derived N2O emissions for a combination of mesofaunal species (Thakur, Groenigen, Kuiper, & Deyn, 2014) to the absence of any interactive effects between earthworms and mesofauna (Schorpp et al., 2016; Wu, Lu, Lu, Guan, & He, 2015). Furthermore, most soil fauna–GHG emission studies to date have focussed on combinations of very few (up to four) different species under controlled laboratory conditions. Based on modelling studies (De Vries et al., 2013), it has been argued that changes in soil biodiversity in the field have substantial consequences for ecosystem processes, such as soil‐derived GHG emissions. However, no controlled laboratory studies so far have looked into the potential role of community composition and functional diversity of soil fauna in increasing or decreasing net GHG emissions from soil. In order to establish that changes in soil biodiversity can indeed influence the production and emission of GHG emissions, a proof‐of‐principle incubation study is needed to focus on causal relationships between soil fauna diversity and GHG emissions from soil.

Here, we quantified biodiversity effects on soil‐derived emissions of N2O and CO2 under controlled laboratory conditions. Methane was not measured because the experimental conditions representative of temperate managed soils (aerobic, relatively low ambient temperature) were not favourable to significant CH4 production (Levy et al., 2012). We experimentally varied species composition in hay‐amended soil microcosms with an increasing number of detritivorous/fungivorous soil mesofauna and macrofauna species: none, one, two, four and eight species per microcosm. We randomly assigned eight species belonging to four different taxonomic groups (earthworms, potworms, mites, springtails) to multispecies treatments following a ‘constrained random sampling design’ (Heemsbergen et al., 2004). Given the collective evidence that a higher species richness enhances fundamental ecological soil processes (Cardinale et al., 2006; Nielsen, Ayres, Wall, & Bardgett, 2011; Tresch et al., 2019; Wall, 2012), we hypothesized increasing soil fauna biodiversity to enhance the microbial (transformation) processes that underlie soil‐derived emissions of N2O and CO2. To investigate this hypothesis, we measured N2O and CO2 emissions over an experimental period of 120 days and calculated net biodiversity effects (Heemsbergen et al., 2004; Loreau & Hector, 2001). We also determined soil biogeochemical properties following the destructive sampling of microcosm soil and calculated the functional dissimilarity among species (Heemsbergen et al., 2004).

MATERIALS AND METHODS

Experimental design

We manipulated faunal species composition in hay‐amended soil microcosms at two levels: species and taxa level. We set the number of soil mesofauna and macrofauna at none, one, two, four and eight species per microcosm treatment (Table 1). The eight selected detritivorous/fungivorous soil fauna species belonged to four different taxonomic groups (earthworms, potworms, mites, springtails; Table 1) and were either obtained from cultures (Vrije Universiteit, Amsterdam, the Netherlands; Wageningen University, Wageningen, the Netherlands) or from grass fields on campus at Wageningen University. To determine compositional effects on soil‐derived GHG fluxes (N2O and CO2), different two‐ and four‐species combinations were included following a ‘constrained random sampling design’ (Heemsbergen et al., 2004). This means that each species was assigned randomly to multispecies treatments with two constraints: (a) each species to be represented equally (six times) and (b) two‐ and four‐species combinations containing similar taxonomic group diversity, that is, one versus two taxonomic groups and two versus four taxonomic groups respectively. Single‐species treatments of all eight species were included to quantify their per capita (mg−1 dry weight animal) effect on N2O and CO2 emissions, and these were used for net biodiversity effect and functional dissimilarity calculations (see Section 2.6). Two control treatments, ‘Control Soil’ and ‘Control Hay’, completed the number of treatments to 23 (Table 1).

Table 1

Overview of the fauna treatments, including treatment codes, # species, # taxa and faunal density (# individuals microcosm−1) per treatment

Abbr.	Fauna	Species	# Species	1			2			2				4				4			8
			# Taxa	1			1			2				2				4			4
Pw1	Potworm	Enchytraeus albidus	Faunal combinations	C	Pw1	50	K	Pw1	25	O	Pw2	25		S	Mi1	100	U		Pw1	12	W	Pw1	6
Pw2	Potworm	Enchytraeus crypticus		D	Pw2	50		Pw2	25		Ew1	2			Mi2	100			Sp2	65		Pw2	6
Mi1	Mite	Rhizoglyphus robini		E	Mi1	400	L	Mi1	200	P	Pw1	25			Ew1	1			Ew1	1		Mi1	50
Mi2	Mite	Oppia nitens		F	Mi2	400		Mi2	200		Sp2	130			Ew2	1			Mi1	100		Mi2	50
Sp1	Springtail	Sinella curviseta		G	Sp1	260	M	Sp1	130	Q	Mi1	200		T	Pw1	13	V		Ew2	1		Sp1	33
Sp2	Springtail	Folsomia candida		H	Sp2	260		Sp2	130		Sp1	130			Pw2	13			Mi2	100		Sp2	33
Ew1	Earthworm	Aporrectodea caliginosa		I	Ew1	4	N	Ew1	2	R	Mi2	200			Sp1	65			Sp1	65		Ew1	1
Ew2	Earthworm	Lumbricus rubellus		J	Ew2	2		Ew2	1		Ew2	1			Sp2	65			Pw2	12		Ew2	1

Codes for treatments are alphabetically: A up to W, for, respectively, Control Soil, Control Hay, Pw1, Mi1, Mi2, Sp1, Sp2, Ew1, Ew2, Pw1 & Pw2, etc. Two control treatments were included in the experimental design: Control Soil (treatment code = A) and Control Hay (treatment code = B), completing the number of treatments to 23.

Microcosm set‐up

Polypropylene microcosms (diameter = 6.7 cm, height = 15 cm, volume = 500 cm3) were filled with 250 g dry soil. The sandy soil (typic Endoaquoll, 75% sand, 23% silt and 2% clay) was collected from ‘Droevendaal Agricultural Farm’ (51°59′N, 5°39′E) from 0 to 25 cm soil depth and contained 14.6 g total C/kg, 1.3 g total N/kg and had a pH (measured in 0.01 M CaCl2) of 5.4. The soil was sieved over an 8 mm mesh and was heated for 24 hr at 70°C to eliminate mesofauna and macrofauna while maintaining a viable soil microbial population (Kaneda & Kaneko, 2011). Distilled water was added at a rate of 220 ml/kg to reach a water‐filled pore space (WFPS) of 62% after gently pressing the soil to a bulk density of 1.4 g/cm3. The soil depth was 6 cm.

Hay residue (32.8 g N/kg dry matter, 448.5 g C/kg dry matter, C:N ratio = 13.8) was chopped in pieces of ~1 cm and wet autoclaved. To represent a realistic N fertilization of 125 kg N/ha, 1.34 g (of dry hay equivalent) was placed on top of the soil surface. After 2 weeks of pre‐incubation, the experiment was started by introducing the soil fauna species to the microcosms according to densities shown in Table 1. Both potworm species were added in an equal ratio of adults to juveniles. Mites (Rhizoglyphus robini) could only be added with some yeast flakes from their cultures; therefore, all microcosms received some yeast flakes. Earthworms were prepared following the wet filter paper method, that is, earthworms were moved to damp filter paper to void gut contents for 48 hr before weighing (Dalby, Baker, & Smith, 1996). All added earthworms, potworms, mites and springtails were added individually by hand with tweezers or a brush.

To avoid intraspecific competition (Dias et al., 2013), the total abundance of the soil mesofauna was kept constant across all treatments, according to established soil food web numbers (De Ruiter, Van Veen, Moore, Brussaard, & Hunt, 1993). For individual mesofauna species, this meant that we added abundances corresponding to their share in the treatment (100% for single‐species treatment, 50% for two‐species combinations, 25% for the four‐species combination and 12.5% for the eight‐species combination, as indicated in Table 1). The total abundance of the macrofauna, that is, the earthworms, was also kept constant across treatments. However, because of their large body size and low abundance, we assigned earthworms on the basis of their body mass rather than their abundance (Table 1, Table S1).

The microcosms were covered tightly with black cloth to prevent fauna from escaping, facilitate gaseous exchange, but minimize moisture loss. Soil moisture was resupplied gravimetrically every 2–3 days after taking gas flux measurements.

The microcosms were arranged in a randomized block design consisting of five replicates divided over five separate blocks; replicates six, seven and eight were placed in block one, two and three, respectively, to be destructively sampled earlier in the experimental period than replicates 1–5. During the experimental period of 120 days, there were two destructive sampling moments: on day 47 (n = 3, replicates six, seven and eight) and on day 120 (n = 5, replicates 1–5). The microcosms were placed in the dark at 16°C and a relative humidity of 60%.

N2O and CO2 flux measurements

Nitrous oxide and CO2 fluxes were measured 34 times over the experimental period of 120 days. Using a static closed chamber technique, gaseous concentrations of N2O and CO2 in the microcosm headspace were measured with an Innova 1302 photo‐acoustic infrared multi gas analyser (Brüel & Kjaer; Velthof, Kuikman, & Oenema, 2002). For N2O flux measurements, a soda‐lime filter was used to minimize interference by CO2; CO2 headspace concentration was therefore measured separately from N2O headspace concentration.

Cumulative N2O and CO2 emissions were calculated assuming a linear change in emission rates between subsequent flux measurements, which was checked occasionally. Values for N2O and CO2 were transformed to CO2 equivalents (CO2‐eq; IPCC, 2013), using a 100 year time horizon as in the Kyoto Protocol. The conversion of N2O emission rates to CO2 equivalents was done as follows: 1 kg N2O = 298 kg CO2 (Lubbers, Van Groenigen, et al., 2013).

Hay decomposition rates were determined by calculating the difference in CO2 production (mg CO2‐C/microcosm) between the control treatment with (Control Hay: ‘B’) and without hay residues (Control Soil: ‘A’). The cumulative hay‐derived decomposition was calculated to confirm that the experiment had proceeded within realistic respiration rates for laboratory conditions (Heemsbergen et al., 2004; Kuiper et al., 2013; Porre, Groenigen, Deyn, Goede, & Lubbers, 2016; Thakur et al., 2014). Hay decomposition rates have not been used in further analyses of faunal effects, since soil fauna diversity effects on these decomposition rates would be similar to effects on cumulative CO2 emissions (see Section 2.6 for the determination of soil fauna diversity effects).

Soil analyses

Destructive sampling for soil property analyses was done on day 1 (n = 4), day 47 (n = 3) and on day 120 (n = 4, randomly block 2–5). Each microcosm was divided vertically in two halves, one of which was used for faunal extractions and the other for soil chemical and microbial analyses. Because of the extraction of the soil fauna from the hay–soil interface, we did not remove and weigh the leftover hay residues from the microcosms on day 47 or day 120. Since in many cases the residues were much stuck to the soil and were partly or entirely mixed in, we carefully removed any left‐over residues from the vertical halves that were kept aside for soil analyses. After thoroughly mixing this latter half, a subsample was used to determine microbial biomass nitrogen (MBN), following the chloroform fumigation and extraction technique (using a correction factor of 0.54; Brookes, Landman, Pruden, & Jenkinson, 1985). Subsequently, total dissolved N (Nts), ammonia (), nitrate and nitrite ( + ) concentrations were measured colorimetrically in a 0.5 M K2SO4 extract, at 1:4 w/v, by segmented flow analysis (3x Segmented Flow Analyzer; Klug, 1999); the dissolved organic N (DON) content was calculated by subtracting  + ( + ) from Nts. The pH and dissolved organic C (DOC) were measured in a 0.01 M CaCl2 extract, at 1:10 w/v, by segmented flow analysis (3x Segmented Flow Analyzer; Houba, Temminghoff, Gaikhorst, & Vark, 2000).

Fauna abundances

Soil fauna was extracted from the other vertical half of the microcosm on day 47 (n = 3) and on day 120 (n = 5). When both potworms and mites were present in the microcosms, we subdivided one microcosm halve further in two, following Thakur et al. (2014). This was necessary because of the different extraction techniques for potworms (wet extraction with a Bearmann funnel; Petersen & Luxton, 1982) and mites (dry extraction with a Berlese funnel; Tullgren funnel; Petersen & Luxton, 1982). Earthworms were carefully collected by hand from the entire microcosm. The numbers of live earthworms were recorded and fresh weights were determined following the wet filter paper method (Dalby et al., 1996). Microcosms belonging to the ‘Control Hay’ treatment were treated in the same way as when potworms and mites were present, to ascertain that no faunal infestations had occurred. Fauna abundances at two destructive sampling dates are shown in Table S1.

More information on the determination of fauna dry weight is in Supporting Information.

Calculations of net biodiversity and functional dissimilarity

The net biodiversity effect was calculated following Heemsbergen et al. (2004) and Loreau and Hector (2001). First, for each single species, its per capita effect (g/DW animal) on the cumulative N2O and CO2 emissions after 120 days was determined by subtracting the Control Hay from the respective single‐species treatments; the average dry weight animal over 120 days was used, that is, the average dry weight at the start of the experiment and at the final harvest after 120 days (Table S1). Second, the expected effect of a multispecies treatment (two‐, four‐ and eight‐species combinations) on either N2O or CO2 was calculated as the sum of the per capita effects of its component species. The difference between the expected effect (Y E) and the observed GHG flux (Y O) of a multispecies treatment was then the net biodiversity effect (ΔY): ΔY = Y O − Y E. Positive net biodiversity effects point to facilitative interactions and negative effects point to inhibitive interactions.

Apart from the influence of species richness on net biodiversity effects, it is known that communities composed of functionally dissimilar species are likely to have stronger effects on process rates (Heemsbergen et al., 2004). Therefore, we also calculated net biodiversity effects in relation to mean functional dissimilarity of species in multispecies mixtures in terms of the key controlling factors of N2O and CO2 production. To analyse the role of biodiversity in soil‐derived GHG emissions, we considered the key driving factors for N2O production and emission (concentrations of ,  + , DOC and pH; Granli & Bøckman, 1994; these factors also impact CO2 production and emission) to shape an ‘attribute space’ in which the functional dissimilarity between species could be quantified as the standardized ‘ecological distance’ (Walker, Kinzig, & Langridge, 1999). The variables used to define the ‘attribute space’ were the per capita concentrations of ,  + , DOC and pH in the single‐species treatments (Table S2). Since these variables do not necessarily have the same units, we defined a linear scale between the smallest and largest value for each variable, and the single‐species treatments were placed along the non‐dimensionalized scales. The use of these normalized scales to shape our ‘four‐dimensional attribute space’ allowed an estimate of functional dissimilarity (sum of pairwise ‘ecological distances’ between species in the four‐dimensional space, divided by the number of interspecific distances between fauna species present in a combination) in which all variables had an equal weight in the calculation.

When including the mean pore size as a soil structural element in a new ‘attribute space’, we decided to interchange pH for mean pore size. This was done for two reasons: (a) keep the ‘attribute space’ four‐dimensional because using four variables (instead of more) is optimal for our number of samples; and (b) pH was selected over the other variables because pH was significantly correlated with most of the other variables used for the ‘attribute space’.

Statistical analyses

All statistical analyses were conducted in SPSS, version 22.0 (SPSS Inc.). Our data of the soil‐derived GHG fluxes were much skewed by the presence of earthworms, and therefore, treatment imposed differences were first established by a nonparametric statistical test: the Kruskal–Wallis test was used to establish the difference between treatments with and without earthworms. Subsequently, the with‐ and without‐earthworm data sets were checked for normality and homoscedasticity using the Shapiro–Wilk test and Levene's test respectively. If needed, data were transformed to normality (with the inverse square root in the case of cumulative N2O fluxes) before using standard one‐way analysis of variance (ANOVA) with blocking and post hoc Tukey HSD to quantify the significance of the effects of soil fauna. Our data of the soil properties (, + , MBN, DOC and pH) were also much skewed. The same procedure was repeated; when no transformations to normality appeared to be possible on day 120, the Kruskal–Wallis Test was again used with subsequent pairwise comparisons (Mann–Whitney U Test with a Bonferroni correction) to test for differences between fauna treatments. The treatment ‘Control Soil’ was never part of statistical analysis.

To test for differences in the net biodiversity effect (on cumulative N2O and CO2 emissions after 120 days) between two‐species, four‐species and eight‐species combinations, an independent samples t test assuming unequal variances was used. A one‐sample t test was used to test the difference between net biodiversity effects and the constant zero; not significantly different from zero indicates a neutral net biodiversity effect.

To identify the relationship between the mean functional dissimilarity (predictor variable) and the net biodiversity effect on cumulative N2O and CO2 (response variables), we used a fitted linear regression model.

For all analyses, a p‐value of .05 or smaller was considered significant.

RESULTS

Net biodiversity effects on faunal‐induced GHG emissions

Figure 1 shows the net biodiversity effect on per capita soil‐derived emissions of N2O and CO2 in relation to species number. The 0‐line indicates a neutral net biodiversity effect. Deviations from this 0‐line (i.e. positive or negative interactions among species) could be used to explain the observed variation in soil‐derived emissions of N2O and CO2 within multispecies treatments. The net biodiversity effects on per capita soil‐derived emissions of N2O and CO2 showed contrasting patterns: on average, species richness reduced cumulative faunal‐induced N2O emissions, whereas it increased cumulative faunal‐induced CO2 emissions (Figure 1). For N2O, we found significant differences between the two‐ and eight‐species treatments (p = .013) and the four‐ and eight‐species treatments (p = .001). For CO2, there were no significant differences between the two‐, four‐ and eight‐species groups. However, for both N2O and CO2, the net biodiversity effects in the four‐ and eight‐species groups were significantly different from the constant 0: net biodiversity effects for N2O in the four‐ and eight‐species treatments were negative, pointing to inhibitive interactions; net biodiversity effects for CO2 in the four‐ and eight‐species treatments were positive, pointing to facilitative interactions (Figure 1).

Figure 1

Net biodiversity effect on per capita soil‐derived emissions of N2O and CO2 (mg N2O‐N m−2 mg−1 DW and g CO2‐C m−2 mg−1 DW, respectively) in relation to species number. The 0‐line indicates a neutral net biodiversity effect. Each box plot shows the 5th and 95th percentiles and the mean (red line) of all treatments for two species (n = 40), four species (n = 20) and eight species (n = 5), respectively. Data points that lie outside the 5th and 95th percentiles are shown as dots. For N2O, significant differences between the two‐ and eight‐species treatments (two‐sample t test assuming unequal variances, t statistic = 2.61, p two‐tail = .013) and the four‐ and eight‐species treatments (two‐sample t test assuming unequal variances, t statistic = 3.70, p two‐tail = .001) are indicated by different letters. For CO2, there are no significant differences between the two‐, four‐ and eight‐species groups. For both N2O and CO2, the four‐ and eight‐species groups are significantly different from the constant 0 (established by a one‐sample t test); levels of significance are *<0.05; **<0.01; ***<0.001; ns = not significant

Figure 2 shows the net biodiversity effect on per capita soil‐derived emissions of N2O and CO2 in relation to mean functional dissimilarity. Similar to the effects of net biodiversity in relation to species number (Figure 1), we found contrasting patterns for net biodiversity effects in relation to mean functional dissimilarity: a negative relationship for soil‐derived N2O emission and a positive relationship for soil‐derived CO2 emission with mean functional dissimilarity (both significant; Figure 2). For N2O, we found a significant negative regression between mean functional dissimilarity and the net biodiversity effect (p = .030). For CO2, we found a significant positive regression between mean functional dissimilarity and the biodiversity effect (p < .001). Since the functional dissimilarity of species mixtures was neither related to species number nor number of taxa, we could assume that it is the divergent influences of the different detritivorous/fungivorous species on the controlling factors of N2O and CO2 production that explain the negative net biodiversity effects for N2O and the positive net biodiversity effects for CO2 emission.

Figure 2

Net biodiversity effect on per capita soil‐derived emissions of N2O and CO2 (mg N2O‐N m−2 mg−1 DW and g CO2‐C m−2 mg−1 DW, respectively) in relation to mean functional dissimilarity (sum of effects on NH4, NO3, DOC and pH (Granli & Bøckman, 1994) of species in the community). The 0‐line indicates a neutral net biodiversity effect. Each vertical series of dots represent a treatment (n = 5 replicates per treatment; some dots overlap). For N2O, a significant negative regression between mean functional dissimilarity and the net biodiversity effect (linear regression, F 64, 63 = 4.95, p = .030) indicates that communities with functionally dissimilar species are more likely to have negative net biodiversity effects. For CO2, a significant positive regression between mean functional dissimilarity and the net biodiversity effect (linear regression, F 64, 63 = 14.5, p < .001) indicates that communities with functionally dissimilar species are more likely to have positive biodiversity effects. Uppercase letters at the top of the figure refer to the species combination given in Table 1. Colours indicate the number of species present: blue = 2 species present, pink = 4 species present and green = 8 species present. Underlined uppercase letters refer to species combinations that include an earthworm species

Community composition effects on faunal‐induced GHG emissions

Figure 3 shows the observed cumulative emissions of N2O and CO2 after 120 days. Faunal‐induced N2O and CO2 emissions in our experiment were dominantly triggered by the presence of large ecosystem engineering detritivores, the two earthworm species, Aporrectodea caliginosa and Lumbricus rubellus. We observed considerably increased rates of N2O emissions for treatments that included earthworms, especially the endogeic earthworm A. caliginosa (Figure 3, the black bars; and Table S3). Increased rates of CO2 emission in single‐species earthworm treatments (treatments ‘I’ and ‘N’ for A. caliginosa and L. rubellus, respectively) were similar for both earthworm species, and much smaller than the increased rates of N2O. The extremely increased N2O emission corresponded with equally increased concentrations of mineral N ( and/or ); after the experimental period of 120 days, there were no significant earthworm (or any other faunal species) effects on dissolved organic C (Table S2).

Figure 3

Observed cumulative emissions of N2O and CO2 after 120 days (in mg CO2‐eq per microcosm). A nonparametric test (Independent‐samples Mann–Whitney U test) established the difference in cumulative emissions from treatments with earthworms present (black bars ‘+EW’) compared to treatments without earthworms present (grey bars ‘−EW’; black bars vs. grey bars: p < .001 for both N2O and CO2). Lowercase letter denotes differences between treatments with earthworms; any other differences occurring are reported in Table S3. Capital letters placed underneath each bar refer to the species combination given in Table 1. Colours indicate the number of species present: blue = 2 species present, pink = 4 species present and green = 8 species present. Underlined uppercase letters refer to species combinations that include an earthworm species

DISCUSSION

Increased species richness and functional dissimilarity result in contrasting faunal‐induced emissions of N2O and CO2

Considering the major contribution of soil fauna species to belowground processes, our findings emphasize the importance of soil fauna functional diversity in the soil function–climate feedback. We found an unexpected contrasting pattern for net biodiversity effects on soil‐derived emissions of N2O and CO2: N2O emissions became smaller than expected when species number increased, whereas CO2 emissions became larger. Consequences of soil fauna diversity for total soil‐derived GHG emissions, therefore, can only be evaluated when understanding the net effects of reduced N2O and enhanced CO2 emission.

For CO2 emissions, the results are corresponding to what is generally known about the effects of soil fauna species richness on ecosystem functioning: increasing species richness generally speeds up the rates of ecosystem processes, such as litter decomposition and/or soil respiration and nutrient cycling (Bardgett & Van der Putten, 2014; Bender, Wagg, & Heijden, 2016; Handa et al., 2014; Nielsen et al., 2011), especially at the low end of the diversity spectrum, as is the case in our controlled experimental set‐up. However, the decreasing N2O emissions with increasing species richness are remarkably opposite to this well‐studied pattern. The explanation might lie in the fundamentally different nature of N2O and CO2 emissions. Whereas CO2 is produced as the end product of respiration and decomposition processes, N2O is produced as an intermediate or by‐product in heterotrophic microbial N transformation processes (Wrage‐Monnig et al., 2018). Therefore, any produced N2O gas could be consumed again and further converted to elemental N2 (Klefoth et al., 2012). Mesofauna and macrofauna activity could potentially have both facilitative and inhibiting effects on the soil microbial transformation processes underlying N2O emission. These constitute direct effects by, for example, stimulating the soil denitrifier community (Drake & Horn, 2006), or indirect effects by, for example, altering the in situ (an)aerobic conditions and influencing the C and N substrate availability (De Ruiter, Moore, et al., 1993; Schroter, Wolters, & Ruiter, 2003). For instance, net biodiversity effects on per capita MBN in relation to species number showed a similar pattern compared to the net biodiversity effects on N2O emission (Figure S1); this could mean that the changes in MBN resulting from increasing species numbers are directly underlying the biodiversity effects found for the per capita N2O emissions. Whether the changes in MBN can be related to increased/decreased activity of N‐related microbes, or whether communities of N‐related microbes were altered by increasing species richness is impossible to deduce from Figure S1.

Soil engineering species such as earthworms and potworms may also cause indirect effects by changing soil structural properties that affect gas diffusion in and out of the soil (Lee & Foster, 1991; Lubbers et al., 2011; Porre et al., 2016). All these faunal interactions with multiple microbial transformations converge in ‘hotspots’ and ‘hot moments’ of microbial activity (Kuzyakov & Blagodatskaya, 2015), and may shift the complex balance of N2O production and consumption towards less N2O production or more complete denitrification, and therefore to reduced rates of faunal‐induced N2O emission.

Next to establishing the influence of species richness on faunal‐induced emissions of CO2 and N2O, we determined the effect of functional dissimilarity of species on CO2 and N2O in our multispecies mixtures. Production of CO2 was most dominant in combinations including ecosystem engineering species (earthworms and potworms; Figure 2 combinations K, N, O, R, S, U, V, W) and therewith corresponded with previous findings by Heemsbergen et al. (2004). The observed inhibition of N2O was likewise mostly influenced by species combinations that included ecosystem engineers, foremost the two earthworm species. As in our experiment the net biodiversity effects were explained most by soil ecosystem engineering traits (i.e. measurable effects on soil structural properties, Capowiez, Sammartino, & Michel, 2011; Lubbers et al., 2011; Porre et al., 2016), we included a soil structural element in our calculations of mean functional dissimilarity. Within the functional dissimilarity matrix, we interchanged pH for mean pore size (as measured by soil X‐ray micro computed tomography; see Porre et al., 2016 for information about the X‐ray tomography procedures and used software and programming, and see Sections 2 and 2.6 for information on calculations of mean functional dissimilarity). Especially the regression between the net biodiversity effect on N2O and the mean functional dissimilarity of species in the community improved when mean pore size was included (F 64, 63 = 8,56, p = .005, Figure S2). For CO2, the regression model did not improve. This indicates that ecosystem engineering traits in soil animals are more important for N2O emission than for CO2 emission.

A similar point can be made when designing the dissimilarity matrix with MBN instead of pH and DON instead of NO3, thereby shifting the matrix more towards a focus on N cycling: the regression between the net biodiversity effect on N2O and the mean functional dissimilarity of species in the community became statistically stronger for N2O when including MBN and DON in the dissimilarity matrix (F 64, 63 = 8.79, p = .004, Figure S3). For CO2, the regression model did not improve. The effects of soil fauna on N cycling properties thus also add explanatory power to the relationship between net biodiversity and mean functional dissimilarity. The reason why this relationship became stronger for N2O and not for CO2 is likely because the MBN varied mostly between the single‐taxum earthworm treatments and the other non‐earthworm fauna treatments (see Table S2), and this was also especially true for N2O emissions (and to a lesser degree for CO2; see Figure 3; Table S3). However, in what ways (size and activity of the N‐related microbial community, as well as its community composition) faunal effects on MBN in the bulk soil influences N2O emission remains elusive in this study.

To establish the influence of species richness on MBN as a response variable, we determined the effect of functional dissimilarity of species on MBN in our multispecies mixtures (see Figure S4 with a dissimilarity matrix of NH4, NO3, DOC and pH; and Figure S5 with a dissimilarity matrix of NH4, NO3, DOC and pore size). The net biodiversity effect on MBN in relation to mean functional dissimilarity shows a significant negative regression in both cases (F 51, 50 = 4.81, p = .033 in Figure S4 and F 51, 50 = 4.48, p = .039 in Figure S5). Like the net biodiversity effect on MBN in relation to species number, this pattern of reduced MBN with increasing functional dissimilarity follows the responses of N2O emission in relation to species number and functional dissimilarity. However, when the soil structural element ‘pore size’ is included in the dissimilarity matrix, the statistical relationship between net biodiversity effect on MBN in relation to mean functional dissimilarity does not become stronger. This corroborates our expectation that ecosystem engineering traits (i.e. effects on soil structure such as pore size) in soil animals are especially important to consider in N2O emissions.

Community composition is key to explaining faunal‐induced GHG emissions

Increased GHG emissions (N2O and CO2) by earthworm activity from intensively managed soils have been well established (Lubbers, Van Groenigen, et al., 2013). Mechanistic explanations for earthworm‐induced N2O emission by previous studies focus on increased N mineralization, nitrification and dentrification, as well as on a changed pore size distribution and an improved connectivity of soil pores with the atmosphere in the presence of earthworms (Bradley, Whalen, Chagnon, Lanoix, & Alves, 2011; Lubbers et al., 2011; Lubbers, González, Hummelink, & Groenigen, 2013; Nebert et al., 2011; Postma‐Blaauw et al., 2006). Here, we show that earthworms increased N mineralization and activated N transformation processes, while C substrate to facilitate heterotrophic microbial processes did not become limiting. However, when we increased biodiversity of other fauna in the presence of earthworms (by increasing species richness as well as soil functional dissimilarity), faunal‐induced N2O emissions decreased despite enhanced C cycling.

Implications for biodiversity effects on soil GHG mitigation

Different hierarchical levels of soil biodiversity may have differential influences on key soil ecosystem processes (Kardol, Throop, Adkins, & Graaff, 2016). However, body size of soil fauna may span the entire hierarchical biodiversity effect, and large soil ecosystem engineers such as earthworm species can dominate soil processes and have disproportional impact on GHG emissions. Therefore, we considered the earthworm effects on soil‐derived N2O and CO2 apart from biodiversity effects by relating observed cumulative GHG emissions to earthworm share (in %, based on earthworm biomass) in the species combinations (Figure S6). The earthworm share is positively related to the effect of earthworm activity on CO2 emissions (linear regression, p < .001, see Figure S7). This is in accordance with the mass ratio hypothesis (Grime, 1998), which states that a given ecosystem process is strongly influenced by the characteristics of a dominant species—proportional to its relative dominance in the community. The earthworm share is also positively related to the effect of earthworm activity on N2O emissions (quadratic regression, p < .001, see Figure S7); however, in the case of N2O emissions, the effects of earthworm presence are nonadditive and are not predictable from the sum of single species. This is indicative of a dominant effect of a unique trait. In our study, it could mean that the diminishing of earthworm‐induced N2O emissions is enhanced by increased variation in species trait values in the community; so in fact by increased functional biodiversity.

Bílá et al. (2014) concluded that dominant trait values in the community play a chief role in driving ecosystem processes such as C cycling processes (so corroborating the mass ratio hypothesis); it is yet to be determined whether the generality of their findings also hold true for other ecosystem processes. Here, we find that ecosystem engineers like our two earthworm species indeed possess dominant trait values concerning emissions of N2O and CO2. However, while for CO2, the mass ratio hypothesis was supported; for N2O, we have shown that functional diversity (indicative of the variation in species trait values in the community) promotes a non‐linear decrease of faunal‐induced N2O emissions.

POTENTIAL CAVEATS AND OUTLOOK

The extreme earthworm‐N2O effects observed in our experiment may partly be a consequence of our microcosm set‐up. Earthworm numbers in the microcosms ranged from 1,140/m2 to 570/m2 for A. caliginosa and L. rubellus, respectively, in single‐species treatments, to 285/m2 for both earthworm species in the eight‐species treatment. It is possible to have very high abundances of earthworms in both natural and cultivated environments; almost as high as in the single‐species treatments of our microcosm set‐up (Schon, Mackay, Gray, Koten, & Dodd, 2017). However, especially when A. caliginosa was present in a single‐species treatment (treatment ‘I’), or combined with L. rubellus (treatment ‘N’), its presence could result in the collapse of soil structure. The microcosm soil could subsequently become compacted throughout, and partly anaerobic (likely going hand‐in‐hand with dramatic increases in the WFPS locally), but with the earthworm burrows mostly conserved. With this severe impact on soil structure, A. caliginosa may have caused its extreme N2O effect by creating perfect in situ conditions for N2O production, and simultaneously, a surface‐connected pore network by which N2O could easily escape to the atmosphere.

High earthworm abundances are ultimately associated with grassland systems (Schon et al., 2017). However, our microcosm set‐up did not allow for plant growth (in the dark at 16°C and covered with cloth to keep the introduced soil animals alive and inside) and our results are therefore valid for bulk soil without an active rhizosphere. When an active rhizosphere is present, it is plausible that it has a substantial influence on earthworm‐induced N2O and CO2 emissions. On the one hand, active plant roots take up nutrients such as N‐ and N‐ for growth (Kraiser, Gras, Gutiérrez, González, & Gutiérrez, 2011), and thereby (partly) remove the substrate for the manifold nitrate‐reducing bacteria and denitrifiers that carry out the N transformation processes underpinning N2O production. On the other hand, it is also to be expected that more easily available C is excreted in root exudates that are produced by active plant roots in the rhizosphere (Jones, Nguyen, & Finlay, 2009). Related consequences of an active rhizosphere may be enhanced microbial activity (increased C substrate/increased C availability), changing soil water and O2 dynamics (give‐and‐take of water by roots as well as O2 diffusion from the soil–atmosphere interface), and likely an interchanging use of the (bio)pore network by roots and soil fauna (perhaps slowing down diffusion rates of N2O to the atmosphere and thereby creating more possibilities for N2O to be reduced to N2). Such consequences may especially cascade down in soil fauna diversity effects on the microbial processes (nitrification, denitrification, nitrifier‐denitrification, decomposition) that underlie N2O and CO2 emissions. It is then to be expected that fauna diversity effects on N2O and CO2 are smaller or even occluded in an active rhizosphere. Yet, previous experimental work by Lubbers et al. (2011; Lubbers, González, et al., 2013) confirmed the occurrence of earthworm‐induced effects on N2O emissions in the presence of growing plants (Lolium perenne). These previously found effects could vary enormously, especially under field conditions: depending on the season and on earthworm species, earthworms could increase N2O by 0%–394% (Lubbers, González, et al., 2013).

A potential unexpected source of crucial substrates for N2O production is the extra N and C that entered the microcosm soil following the death of a number of earthworms. However, their dead biomass or emaciation during the experimental period (Table S1) was unlikely to have been the cause of either extremely increased or less increased rates of N2O emission. Not only did few earthworms die but also there was no relationship between replicate microcosms in which earthworms had died and the N2O or CO2 emissions from the corresponding microcosms. In a previous study, Lubbers et al. (2011) calculated the highest possible contribution of N released from decaying earthworm tissue and concluded that earthworm‐induced N2O emissions could not have resulted from dead earthworms.

These considerations in relation to potential caveats accompanying microcosm studies such as ours call for experimental work under realistic conditions including growing plants, prevailing weather conditions, and up‐to‐date agricultural management practices. Microcosm incubation studies that confirm ‘proof of principle’, namely, in our study effects of soil fauna diversity on emissions of CO2 and N2O, are a prerequisite for extended investigations that can be successfully extrapolated to realistic field conditions.

To conclude, we found that higher species richness and functional dissimilarity led to increased faunal‐induced CO2 emission but decreased N2O emission. These contrasting dynamics are best explained by the fundamentally different processes underlying CO2 and N2O emissions: CO2 is produced as the end product of respiration and decomposition processes, while N2O is produced as an intermediate or by‐product of microbial N transformation processes; during upward diffusion, any produced N2O can be consumed again and be converted to elemental N2. Community composition appeared to determine the observed cumulative GHG emissions over an experimental period of 120 days. The presence of keystone species such as earthworms greatly enhanced soil‐derived emissions of N2O, but increasing species richness and functional dissimilarity reduced these elevated earthworm‐induced N2O emissions. These effects in our study call for further investigations that may contribute to the adoption of management practices that enhance soil biodiversity and stimulate a functionally diverse faunal community (rather than stimulating e.g. only target species such as earthworms) as to reduce N2O emissions from managed soils.

Taken together, our findings recognize the importance of soil biodiversity conservation and call for incorporation of methods to maintain and stimulate the diversity of the living biota that provide critical processes and ecosystem services. New research may be directed towards a better mechanistic understanding of biodiversity effects on in situ soil structure including an active rhizosphere, with implications for gas diffusion and microbial diversity and activity. Within the contexts of ‘soil ecological engineering’ (Bender et al., 2016) and ‘climate‐smart soils’ (Paustian et al., 2016), we need knowledge on faunal diversity effects connected to the potential of soils to mitigate GHG emissions that can be applied to achieve sustainable agriculture and viable agricultural systems.

AUTHOR CONTRIBUTIONS

IML and JWG conceived and designed the study, with suggestions and input from MPB and GBDD; IML and GBDD set up and destructively sampled the experiment; IML measured all GHG fluxes, maintained the experiment, did laboratory analyses and performed the statistical analyses; IML, JWG, MPB, WHP and GBDD interpreted and discussed the results; IML, JWG and WHP wrote the paper, with substantial contributions from all co‐authors; JWG had the overall supervision of the project.

Click here for additional data file.

Click here for additional data file.