# Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment

Pierre Ganault, Johanne Nahmani, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, Isabelle Bertrand, Bruno Buatois, Alexandru Milcu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Ganault, Johanne Nahmani, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, et al., Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment. PLoS ONE, 2024, 19 (2), pp.e0289859. 10.1371/journal.pone.0289859. [hal-04552879](https://hal.inrae.fr/hal-04552879)

HAL Id: hal-04552879

https://hal.inrae.fr/hal-04552879

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment

Pierre Ganault, Johanne Nahmani, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, Isabelle Bertrand, Bruno Buatois, Alexandru Milcu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Ganault, Johanne Nahmani, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, et al., Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment. PLoS ONE, 2024, 19 (2), pp.e0289859. 10.1371/journal.pone.0289859. [hal-04552879v2](https://hal.inrae.fr/hal-04552879v2)

HAL Id: hal-04552879

https://hal.inrae.fr/hal-04552879v2

Submitted on 19 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment

Pierre Ganault<sup>1,2,3,4*</sup>, Johanne Nahmani<sup>2</sup>, Yvan Capowiez<sup>5</sup>, Nathalie Fromin<sup>2*</sup>, Ammar Shihan<sup>2</sup>, Isabelle Bertrand<sup>1,6</sup>, Bruno Buatois<sup>2</sup>, Alexandru Milcu<sup>2,7</sup>

1 ECODIV, INRAE, Normandie Université, UNIROUEN, Rouen, France, 2 CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France, 3 German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany, 4 Institute of Biology, Leipzig University, Leipzig, Germany, 5 INRAE, UMR 1114 EMMAH, INRAE/Université d'Avignon, Site Agroparc, Avignon, France, 6 UMR Eco&Sols, CIRAD, INRAE, IRD, Montpellier SupAgro, Université de Montpellier, Montpellier, France, 7 Montpellier European Ecotron, Univ Montpellier, CNRS, Campus Baillarguet, Montferrier-sur-Lez, France

<sup>¤</sup> Current address: CNRS–7 rue du Four Solaire, Odeillo, France

* pierre.ganault@univ-rouen.fr

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

## OPEN ACCESS

**Citation:** Ganault P, Nahmani J, Capowiez Y, Fromin N, Shihan A, Bertrand I, et al. (2024) Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment. PLoS ONE 19(2): e0289859. https://doi.org/10.1371/journal.pone.0289859

**Editor:** Farhan Hafeez, COMSATS University Islamabad Abbottabad Campus, PAKISTAN

**Received:** July 27, 2023

**Accepted:** November 23, 2023

**Published:** February 15, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0289859

**Copyright:** © 2024 Ganault et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The datasets and R code of the current study is currently available on

## Abstract

Earthworms can stimulate microbial activity and hence greenhouse gas (GHG) emissions from soils. However, the extent of this effect in the presence of plants and soil moisture fluctuations, which are influenced by earthworm burrowing activity, remains uncertain. Here, we report the effects of earthworms (without, anecic, endogeic, both) and plants (with, without) on GHG (CO₂, N₂O) emissions in a 3-month greenhouse mesocosm experiment simulating a simplified agricultural context. The mesocosms allowed for water drainage at the bottom to account for the earthworm engineering effect on water flow during two drying-wetting cycles. N₂O cumulative emissions were 34.6% and 44.8% lower when both earthworm species and only endogeic species were present, respectively, and 19.8% lower in the presence of plants. The presence of the endogeic species alone or in combination with the anecic species slightly reduced CO₂ emissions by 5.9% and 11.4%, respectively, and the presence of plants increased emissions by 6%. Earthworms, plants and soil water content interactively affected weekly N₂O emissions, an effect controlled by increased soil dryness due to drainage via earthworm burrows and mesocosm evapotranspiration. Soil macroporosity (measured by X-ray tomography) was affected by earthworm species-specific burrowing activity. Both GHG emissions decreased with topsoil macropore volume, presumably due to reduced moisture and microbial activity. N₂O emissions decreased with macropore volume in the deepest layer, likely due to the presence of fewer anaerobic microsites. Our results indicate that, under experimental conditions allowing for plant and earthworm engineering effects on soil moisture, earthworms do not increase GHG emissions, and endogeic earthworms may even reduce N₂O emissions.
RESEARCH ARTICLE

# Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment

Pierre Ganault<sup>1,2,3,4*</sup>, Johanne Nahmani<sup>2</sup>, Yvan Capowiez<sup>5</sup>, Nathalie Fromin<sup>2*</sup>, Ammar Shihan<sup>2</sup>, Isabelle Bertrand<sup>1,6</sup>, Bruno Buatois<sup>2</sup>, Alexandru Milcu<sup>2,7</sup>

1 ECODIV, INRAE, Normandie Université, UNIROUEN, Rouen, France, 2 CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France, 3 German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany, 4 Institute of Biology, Leipzig University, Leipzig, Germany, 5 INRAE, UMR 1114 EMMAH, INRAE/Université d'Avignon, Site Agroparc, Avignon, France, 6 UMR Eco&Sols, CIRAD, INRAE, IRD, Montpellier SupAgro, Université de Montpellier, Montpellier, France, 7 Montpellier European Ecotron, Univ Montpellier, CNRS, Campus Baillarguet, Montferrier-sur-Lez, France

<sup>¤</sup> Current address: CNRS–7 rue du Four Solaire, Odeillo, France

* pierre.ganault@univ-rouen.fr

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Ganault P, Nahmani J, Capowiez Y, Fromin N, Shihan A, Bertrand I, et al. (2024) Earthworms and plants can decrease soil greenhouse gas emissions by modulating soil moisture fluctuations and soil macroporosity in a mesocosm experiment. PLoS ONE 19(2): e0289859. https://doi.org/10.1371/journal.pone.0289859

**Editor:** Farhan Hafeez, COMSATS University Islamabad Abbottabad Campus, PAKISTAN

**Received:** July 27, 2023

**Accepted:** November 23, 2023

**Published:** February 15, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0289859

**Copyright:** © 2024 Ganault et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The datasets and R code of the current study is currently available on

## Abstract

Earthworms can stimulate microbial activity and hence greenhouse gas (GHG) emissions from soils. However, the extent of this effect in the presence of plants and soil moisture fluctuations, which are influenced by earthworm burrowing activity, remains uncertain. Here, we report the effects of earthworms (without, anecic, endogeic, both) and plants (with, without) on GHG (CO₂, N₂O) emissions in a 3-month greenhouse mesocosm experiment simulating a simplified agricultural context. The mesocosms allowed for water drainage at the bottom to account for the earthworm engineering effect on water flow during two drying-wetting cycles. N₂O cumulative emissions were 34.6% and 44.8% lower when both earthworm species and only endogeic species were present, respectively, and 19.8% lower in the presence of plants. The presence of the endogeic species alone or in combination with the anecic species slightly reduced CO₂ emissions by 5.9% and 11.4%, respectively, and the presence of plants increased emissions by 6%. Earthworms, plants and soil water content interactively affected weekly N₂O emissions, an effect controlled by increased soil dryness due to drainage via earthworm burrows and mesocosm evapotranspiration. Soil macroporosity (measured by X-ray tomography) was affected by earthworm species-specific burrowing activity. Both GHG emissions decreased with topsoil macropore volume, presumably due to reduced moisture and microbial activity. N₂O emissions decreased with macropore volume in the deepest layer, likely due to the presence of fewer anaerobic microsites. Our results indicate that, under experimental conditions allowing for plant and earthworm engineering effects on soil moisture, earthworms do not increase GHG emissions, and endogeic earthworms may even reduce N₂O emissions.
the first author GitHub repository at the following link https://github.com/PGanault/EarthwormPlantGHG.

**Funding:** The author(s) received no specific funding for this work

**Competing interests:** The authors have declared that no competing interests exist.

Introduction

Soil invertebrates strongly influence organic matter dynamics, nutrient cycling, and water and gas fluxes through their feeding activity, movement through soil layers, and interactions with other organisms [1–3]. Among them, earthworms are essential and promote soil fertility and plant productivity [4, 5], resulting in a 25% increase in crop yield [6]. Conversely, the presence of earthworms may increase CO₂ emissions by 33% and N₂O emissions by 42% [7], which is notable because N₂O has 265 times the global warming potential of CO₂ [8]. However, the microbial processes causing greenhouse gas (GHG) emissions are extremely complex and involve numerous interactions between earthworms, microbial communities, plants, and soil water and aeration status, determining C and N mineralization or stabilization [9–14]. Simultaneously, exploring the complexity of these interactions and mechanisms mentioned above in experimental settings poses a challenge. Most existing studies have focused on scenarios without plants and at constant soil moisture levels that favor earthworm and microbial activity [7]. This limits our understanding of how soil moisture fluctuations, modulated by earthworms and plants and their interactions, affect N₂O and CO₂ emissions [7, 15].

Soil water content (SWC) is a long recognized key factor that explains up to 95% of GHG emissions [16, 17], as it drives microbial processes such as respiration, denitrification, and nitrification that produce GHG [9, 18, 19]. Indeed, SWC determines gas and nutrient diffusion and hence the availability of oxygen, nitrate, ammonium and carbon to microorganisms, thereby modulating their activity. Under anoxic conditions, at high SWC, N₂O emissions are the highest, mostly due to denitrification, while aerobic conditions favor N₂O emissions by nitrification [20, 21]. Similarly, a substantial body of evidence has shown that carbon substrate limitation occurs in drier conditions, while oxygen limitation occurs under conditions close to water saturation, with optimal conditions for respiration and hence CO₂ emissions at intermediate levels of SWC [22, 23]. Soil moisture fluctuations (e.g., drying–rewetting cycles) can also affect the proportion of nitrogen denitrified into N₂O or N₂, thus modulating the N₂O/N₂ ratio that will be emitted into the atmosphere [24, 25]. Therefore, keeping soil moisture constant, as in most existing experiments, limits the occurrence of these mechanisms.

Earthworm feeding and burrowing activity strongly influence carbon and nutrient dynamics and gas and water fluxes in the soil profile. Earthworms' fresh casts are richer in organic matter and water than bulk soil and can promote microbial respiration, while the anoxic conditions in earthworm digestive tracts promote the growth of denitrifying bacteria, two mechanisms leading to higher CO₂ and N₂O emissions, respectively [7, 26, 27]. Moreover, earthworms' aging casts can constitute carbon stocks by the physical protection of particulate organic matter when mixed with the mineral fraction and the increase in microbial necromass, which constitutes another form of stable carbon stock [12, 26, 28–30]. Earthworms also have indirect impacts by burrowing into the soil, changing soil macroporosity and affecting air and water fluxes, soil moisture, soil compaction, and CO₂ and N₂O diffusivity [31]. The effect of earthworms on GHG emissions therefore depends on the balance between these mechanisms and whether they may occur in the experimental set-up.

As primary producers, plants control organic matter quantity and quality in soils [32]. Root exudates contribute more to carbon stabilization than plant litter because the production of labile compounds (e.g., simple sugars) in the mineral soil favors the microbially driven formation of mineral-associated organic matter [14, 33, 34]. Inputs of root-derived C substrates can lead to high transient O₂ demand and can cause suboxic microsites in the rhizosphere, thus favoring denitrification [9, 35]. Conversely, plants compete with microbes for nitrogen acquisition, reduce SWC by transpiration, modify soil porosity by root growth [36, 37], and thus can change the preponderance of the controlling N₂O emission processes (nitrification,
the first author GitHub repository at the following link https://github.com/PGanault/EarthwormPlantGHG.

**Funding:** The author(s) received no specific funding for this work

**Competing interests:** The authors have declared that no competing interests exist.

Introduction

Soil invertebrates strongly influence organic matter dynamics, nutrient cycling, and water and gas fluxes through their feeding activity, movement through soil layers, and interactions with other organisms [1–3]. Among them, earthworms are essential and promote soil fertility and plant productivity [4, 5], resulting in a 25% increase in crop yield [6]. Conversely, the presence of earthworms may increase CO₂ emissions by 33% and N₂O emissions by 42% [7], which is notable because N₂O has 265 times the global warming potential of CO₂ [8]. However, the microbial processes causing greenhouse gas (GHG) emissions are extremely complex and involve numerous interactions between earthworms, microbial communities, plants, and soil water and aeration status, determining C and N mineralization or stabilization [9–14]. Simultaneously, exploring the complexity of these interactions and mechanisms mentioned above in experimental settings poses a challenge. Most existing studies have focused on scenarios without plants and at constant soil moisture levels that favor earthworm and microbial activity [7]. This limits our understanding of how soil moisture fluctuations, modulated by earthworms and plants and their interactions, affect N₂O and CO₂ emissions [7, 15].

Soil water content (SWC) is a long recognized key factor that explains up to 95% of GHG emissions [16, 17], as it drives microbial processes such as respiration, denitrification, and nitrification that produce GHG [9, 18, 19]. Indeed, SWC determines gas and nutrient diffusion and hence the availability of oxygen, nitrate, ammonium and carbon to microorganisms, thereby modulating their activity. Under anoxic conditions, at high SWC, N₂O emissions are the highest, mostly due to denitrification, while aerobic conditions favor N₂O emissions by nitrification [20, 21]. Similarly, a substantial body of evidence has shown that carbon substrate limitation occurs in drier conditions, while oxygen limitation occurs under conditions close to water saturation, with optimal conditions for respiration and hence CO₂ emissions at intermediate levels of SWC [22, 23]. Soil moisture fluctuations (e.g., drying–rewetting cycles) can also affect the proportion of nitrogen denitrified into N₂O or N₂, thus modulating the N₂O/N₂ ratio that will be emitted into the atmosphere [24, 25]. Therefore, keeping soil moisture constant, as in most existing experiments, limits the occurrence of these mechanisms.

Earthworm feeding and burrowing activity strongly influence carbon and nutrient dynamics and gas and water fluxes in the soil profile. Earthworms' fresh casts are richer in organic matter and water than bulk soil and can promote microbial respiration, while the anoxic conditions in earthworm digestive tracts promote the growth of denitrifying bacteria, two mechanisms leading to higher CO₂ and N₂O emissions, respectively [7, 26, 27]. Moreover, earthworms' aging casts can constitute carbon stocks by the physical protection of particulate organic matter when mixed with the mineral fraction and the increase in microbial necromass, which constitutes another form of stable carbon stock [12, 26, 28–30]. Earthworms also have indirect impacts by burrowing into the soil, changing soil macroporosity and affecting air and water fluxes, soil moisture, soil compaction, and CO₂ and N₂O diffusivity [31]. The effect of earthworms on GHG emissions therefore depends on the balance between these mechanisms and whether they may occur in the experimental set-up.

As primary producers, plants control organic matter quantity and quality in soils [32]. Root exudates contribute more to carbon stabilization than plant litter because the production of labile compounds (e.g., simple sugars) in the mineral soil favors the microbially driven formation of mineral-associated organic matter [14, 33, 34]. Inputs of root-derived C substrates can lead to high transient O₂ demand and can cause suboxic microsites in the rhizosphere, thus favoring denitrification [9, 35]. Conversely, plants compete with microbes for nitrogen acquisition, reduce SWC by transpiration, modify soil porosity by root growth [36, 37], and thus can change the preponderance of the controlling N₂O emission processes (nitrification,
denitrification) [9, 38]. Simultaneously, the action of earthworms burrowing through the soil profile increases soil aeration, water drainage and possibly soil drying [39–42], thus potentially creating conditions that are less favorable for denitrification and N₂O emission [15, 43]. The above mechanisms all interact, as earthworms promote plant growth by increasing nutrient availability to plants and by consuming roots [4–6]. To the best of our knowledge, no study has investigated earthworm-mediated soil moisture variation effects on GHG emissions. Furthermore, the vast majority of existing experimental studies have used micro/mesocosms with sealed bottoms, which impedes water drainage and the modulation of water infiltration via earthworm burrows and root growth that would typically occur in realistic field conditions. Two studies that evaluated the impacts of soil moisture fluctuations on GHG in the presence of earthworms, although without plants, showed that cumulative N₂O and CO₂ emissions were reduced in the presence of earthworms [15, 43].

Allowing earthworm burrowing activity to influence SWC could aid in understanding the varying effects of different earthworm species and ecological categories on earthworm-mediated GHG emissions. Indeed, lumbricid earthworms are broadly classified into three main ecological categories based on their feeding and burrowing characteristics: (1) anecic species that feed on surface litter by pulling it into permanent vertical burrows and creating surface casts, (2) epigeic species that live and feed in surface litter, making very few nonpermanent burrows and (3) endogeic species that live in the soil, feed on roots and soil organic matter, and make numerous nonpermanent burrows [44, 45]. Hence, anecic earthworms are likely to have a stronger impact on organic matter redistribution and water fluxes [38]. Conversely, endogeic and epigeic earthworms may primarily influence organic matter redistribution in the soil and at the surface, respectively, with epigeic earthworms having a lesser impact on water fluxes [46]. There is evidence that CO₂ and N₂O emissions depend on the earthworm ecological category, with significantly higher emissions for the anecic group [7], but the net balance between mineralization and stabilization over time was reported to be highly variable within each category [47]. Additionally, whether this finding holds in the presence of plants and soil moisture fluctuations remains to be tested.

In this study, we assessed the impact of four levels of earthworm treatments (one endogeic species, one anecic species, both anecic and endogeic species, and a control without earthworms) and two levels of plant treatments (with or without a model grass species). We used a full factorial design (4 × 2 = 8 treatment combinations, with 7 replicates each) during a three-month greenhouse mesocosm experiment simulating a simplified agricultural setup. The experiment involved simulating two drying–wetting cycles, and to facilitate the earthworm engineering effect on soil water infiltration and status, the mesocosms were designed to enable effective water drainage and escape via percolation. We measured weekly CO₂ and N₂O fluxes, aboveground plant biomass, litter cover and multiple soil parameters representing potentially relevant predictors of GHG emissions, including soil nitrogen and water status, microbial biomass and respiration, denitrification potential and multiple metrics of soil macroporosity, using X-ray tomography. We hypothesized that 1) CO₂ and N₂O emissions will be lower in the presence of earthworms relative to controls, as increased carbon and nitrogen mineralization will be offset by the drier and more aerated conditions due to the earthworm soil engineering effect (burrowing) on water drainage, 2) plant presence will reduce N₂O emissions due to nitrogen and water uptake but will increase CO₂ emissions due to increased carbon substrates entering the soil via rhizodeposition, and 3) differences in N₂O and CO₂ fluxes among the two earthworm ecological categories will be mediated by the burrowing patterns affecting soil (macro) porosity and water status, with lower emissions with higher microporosity expected, as these conditions increase the volumes of aerobic and dryer sites.
denitrification) [9, 38]. Simultaneously, the action of earthworms burrowing through the soil profile increases soil aeration, water drainage and possibly soil drying [39–42], thus potentially creating conditions that are less favorable for denitrification and N₂O emission [15, 43]. The above mechanisms all interact, as earthworms promote plant growth by increasing nutrient availability to plants and by consuming roots [4–6]. To the best of our knowledge, no study has investigated earthworm-mediated soil moisture variation effects on GHG emissions. Furthermore, the vast majority of existing experimental studies have used micro/mesocosms with sealed bottoms, which impedes water drainage and the modulation of water infiltration via earthworm burrows and root growth that would typically occur in realistic field conditions. Two studies that evaluated the impacts of soil moisture fluctuations on GHG in the presence of earthworms, although without plants, showed that cumulative N₂O and CO₂ emissions were reduced in the presence of earthworms [15, 43].

Allowing earthworm burrowing activity to influence SWC could aid in understanding the varying effects of different earthworm species and ecological categories on earthworm-mediated GHG emissions. Indeed, lumbricid earthworms are broadly classified into three main ecological categories based on their feeding and burrowing characteristics: (1) anecic species that feed on surface litter by pulling it into permanent vertical burrows and creating surface casts, (2) epigeic species that live and feed in surface litter, making very few nonpermanent burrows and (3) endogeic species that live in the soil, feed on roots and soil organic matter, and make numerous nonpermanent burrows [44, 45]. Hence, anecic earthworms are likely to have a stronger impact on organic matter redistribution and water fluxes [38]. Conversely, endogeic and epigeic earthworms may primarily influence organic matter redistribution in the soil and at the surface, respectively, with epigeic earthworms having a lesser impact on water fluxes [46]. There is evidence that CO₂ and N₂O emissions depend on the earthworm ecological category, with significantly higher emissions for the anecic group [7], but the net balance between mineralization and stabilization over time was reported to be highly variable within each category [47]. Additionally, whether this finding holds in the presence of plants and soil moisture fluctuations remains to be tested.

In this study, we assessed the impact of four levels of earthworm treatments (one endogeic species, one anecic species, both anecic and endogeic species, and a control without earthworms) and two levels of plant treatments (with or without a model grass species). We used a full factorial design (4 × 2 = 8 treatment combinations, with 7 replicates each) during a three-month greenhouse mesocosm experiment simulating a simplified agricultural setup. The experiment involved simulating two drying–wetting cycles, and to facilitate the earthworm engineering effect on soil water infiltration and status, the mesocosms were designed to enable effective water drainage and escape via percolation. We measured weekly CO₂ and N₂O fluxes, aboveground plant biomass, litter cover and multiple soil parameters representing potentially relevant predictors of GHG emissions, including soil nitrogen and water status, microbial biomass and respiration, denitrification potential and multiple metrics of soil macroporosity, using X-ray tomography. We hypothesized that 1) CO₂ and N₂O emissions will be lower in the presence of earthworms relative to controls, as increased carbon and nitrogen mineralization will be offset by the drier and more aerated conditions due to the earthworm soil engineering effect (burrowing) on water drainage, 2) plant presence will reduce N₂O emissions due to nitrogen and water uptake but will increase CO₂ emissions due to increased carbon substrates entering the soil via rhizodeposition, and 3) differences in N₂O and CO₂ fluxes among the two earthworm ecological categories will be mediated by the burrowing patterns affecting soil (macro) porosity and water status, with lower emissions with higher microporosity expected, as these conditions increase the volumes of aerobic and dryer sites.
# Materials and methods

## Soil and biological material

The soil, classified as a gleyic luvisol, was excavated from a field margin adjacent to a wheat-corn-alfalfa rotation at the EFELE experimental site (Northwest of France, 8°05'35.9"N, 1°48'53.1"W) belonging to the long-term observatories SOERE-PRO-network. Only soil from the upper 0–30 cm layer was used in this experiment. The soil was composed of 14.6% clay, 72.1% silt and 13.3% sand, with a pH of 6.14 and a volumetric water content at field capacity of 39.2% (Soil Analysis Laboratory, INRA Arras, France). The soil contained 1.5% total organic matter, 0.84% carbon, and 0.1% nitrogen, with a C:N ratio of 8.4. The mesocosm setup aimed at reproducing the agricultural field from which it was sampled by using its soil and adding locally present earthworms and a plant species commonly used as a model system for cereal grass.

Adult individuals of *Lumbricus terrestris* L. were supplied by Wurmwelten Company (Dassel, Germany) and weighed 4.8 ± 1.3 g fresh weight on average. Adult individuals of *Aporrectodea icterica* Savigny were harvested from a pesticide-free orchard in Avignon by manual digging and weighed 0.4 ± 0.2 g fresh weight on average. The earthworms were kept in their original soil for 3 days at a temperature of 14 ± 2°C and then placed in a mixture of the original soil and the experimental soil for one week at 8°C in the dark before the onset of the experiment.

*Brachypodium distachyon* L. was the plant species selected for this study due to its small size and short life cycle of less than 3 months [48] and because it is frequently used in controlled environment experiments [49]. The seeds of the wild-type variety (Bd 21 WT) were supplied by Observatoire du Végétal, INRAE Versailles (Paris, France). After one week of germination in seedling trays in vermiculite, four seedlings were planted in each mesocosm outside the central cylinder that was introduced as a base for flux measurements of greenhouse gases (see S1 Fig). During the first 3 weeks, any dead seedlings were replaced. The experiment ended with a final destructive harvest.

## Mesocosm design and experimental treatments

The mesocosms consisted of PVC tubes 16 cm in diameter and 37 cm in height (S1 Fig). Each mesocosm was filled up to 3 cm from the brim with 9.2 kg of soil already containing 10% gravimetric water content, sieved to 2 mm, and compacted to a bulk density of 1.21 g cm⁻³. The mesocosms were sealed at their base with a 1 mm mesh followed by a PVC lid pierced with 5 holes (1 cm in dia.), which allowed the drainage of surplus water out of the mesocosm. A transparent plastic film 10 cm in height was attached around the top perimeter of the mesocosms to prevent the earthworms from leaving them. As a previous meta-analysis indicated that earthworm effects on plant growth are more prevalent in the presence of crop residues, which serve as a food resource for soil biota [6], the soil surfaces of all mesocosms were covered with a 4 g litter mixture (2.2% N, C/N = 24) consisting of 1.3 g dry weight *Medicago truncatula* Gaertn. shoots and 2.7 g dry weight of *Zea mays* L. leaves, the equivalent of organic residue inputs of 1060 kg C ha⁻¹ and 44 kg N ha⁻¹.

The mesocosm experiment presented in this study included four levels of earthworm treatments (henceforth Ew): a control without earthworms, an anecic earthworm species (*L. terrestris*) with two individuals weighing 9.6 ± 1 g fresh weight (FW) on average per replicate, an endogeic earthworm species (*A. icterica*) with 7 ± 1.1 individuals weighing 2.9 ± 0.1 g FW on average per replicate, and a mixture of both species with one *L. terrestris* individual (4.9 ± 0.9 g FW biomass) and 5 ± 1.6 *A. icterica* individuals (1.7 ± 0.5 g FW biomass) per replicate. The
# Materials and methods

## Soil and biological material

The soil, classified as a gleyic luvisol, was excavated from a field margin adjacent to a wheat-corn-alfalfa rotation at the EFELE experimental site (Northwest of France, 8°05'35.9"N, 1°48'53.1"W) belonging to the long-term observatories SOERE-PRO-network. Only soil from the upper 0–30 cm layer was used in this experiment. The soil was composed of 14.6% clay, 72.1% silt and 13.3% sand, with a pH of 6.14 and a volumetric water content at field capacity of 39.2% (Soil Analysis Laboratory, INRA Arras, France). The soil contained 1.5% total organic matter, 0.84% carbon, and 0.1% nitrogen, with a C:N ratio of 8.4. The mesocosm setup aimed at reproducing the agricultural field from which it was sampled by using its soil and adding locally present earthworms and a plant species commonly used as a model system for cereal grass.

Adult individuals of *Lumbricus terrestris* L. were supplied by Wurmwelten Company (Dassel, Germany) and weighed 4.8 ± 1.3 g fresh weight on average. Adult individuals of *Aporrectodea icterica* Savigny were harvested from a pesticide-free orchard in Avignon by manual digging and weighed 0.4 ± 0.2 g fresh weight on average. The earthworms were kept in their original soil for 3 days at a temperature of 14 ± 2°C and then placed in a mixture of the original soil and the experimental soil for one week at 8°C in the dark before the onset of the experiment.

*Brachypodium distachyon* L. was the plant species selected for this study due to its small size and short life cycle of less than 3 months [48] and because it is frequently used in controlled environment experiments [49]. The seeds of the wild-type variety (Bd 21 WT) were supplied by Observatoire du Végétal, INRAE Versailles (Paris, France). After one week of germination in seedling trays in vermiculite, four seedlings were planted in each mesocosm outside the central cylinder that was introduced as a base for flux measurements of greenhouse gases (see S1 Fig). During the first 3 weeks, any dead seedlings were replaced. The experiment ended with a final destructive harvest.

## Mesocosm design and experimental treatments

The mesocosms consisted of PVC tubes 16 cm in diameter and 37 cm in height (S1 Fig). Each mesocosm was filled up to 3 cm from the brim with 9.2 kg of soil already containing 10% gravimetric water content, sieved to 2 mm, and compacted to a bulk density of 1.21 g cm⁻³. The mesocosms were sealed at their base with a 1 mm mesh followed by a PVC lid pierced with 5 holes (1 cm in dia.), which allowed the drainage of surplus water out of the mesocosm. A transparent plastic film 10 cm in height was attached around the top perimeter of the mesocosms to prevent the earthworms from leaving them. As a previous meta-analysis indicated that earthworm effects on plant growth are more prevalent in the presence of crop residues, which serve as a food resource for soil biota [6], the soil surfaces of all mesocosms were covered with a 4 g litter mixture (2.2% N, C/N = 24) consisting of 1.3 g dry weight *Medicago truncatula* Gaertn. shoots and 2.7 g dry weight of *Zea mays* L. leaves, the equivalent of organic residue inputs of 1060 kg C ha⁻¹ and 44 kg N ha⁻¹.

The mesocosm experiment presented in this study included four levels of earthworm treatments (henceforth Ew): a control without earthworms, an anecic earthworm species (*L. terrestris*) with two individuals weighing 9.6 ± 1 g fresh weight (FW) on average per replicate, an endogeic earthworm species (*A. icterica*) with 7 ± 1.1 individuals weighing 2.9 ± 0.1 g FW on average per replicate, and a mixture of both species with one *L. terrestris* individual (4.9 ± 0.9 g FW biomass) and 5 ± 1.6 *A. icterica* individuals (1.7 ± 0.5 g FW biomass) per replicate. The
earthworm FW biomasses were the equivalent of 480, 145 and 330 g m⁻² for the anecic, endogeic and both earthworm treatment levels, respectively, which were 2- to 3-fold higher than that in the field of origin of the soil, where earthworm total biomass varied from 98 to 135 g m⁻². For *L. terrestris*, adding the proportional equivalent of the field biomass would mean introducing just one individual. Given the associated risk that the escape or death of this single individual could jeopardize the treatment, we chose to add two adult individuals of *L. terrestris*. As a result, the overall biomass was roughly double what is typically found per m^2 in agroecosystems. Note that *L. terrestris* was more recently reclassified as a species displaying traits belonging to both anecic and epigeic species [36], but in the vast majority of literature, it is considered an anecic species. The earthworm treatments were factorially crossed with two levels of plant (*B. distachyon* L.) treatments, i.e., with and without plants, with 7 replicates per treatment combination for a total of 56 mesocosms.

The experiment was conducted in a greenhouse kept at temperatures ranging between 20 and 23°C during the day and 18–20°C during the night with an air relative humidity of 80%. The natural light was supplemented during daytime by artificial lighting for 12 hours per day using high-pressure sodium lamps. The mesocosms were divided into two blocks corresponding to their position on the north or south bench of the greenhouse, and their position within the block was randomly changed twice a week, limiting the bias of the position within the block. The experiment ran for 12 weeks between March and May 2017.

## Watering protocol

The watering protocol was specifically designed to include soil moisture fluctuations (analogous to what occurs in natural conditions) and to allow earthworm burrowing to affect SWC. At the beginning of the experiment, the mesocosms were watered with 1.7 L of reverse osmosis water using a laboratory dispenser (two sessions of 850 ml each), a volume sufficient to observe water draining out of the mesocosms from the pierced lids at the bottoms of mesocosms. Measurements of weight changes after 24 h were used to calculate the weights of the mesocosms at field capacity (knowing that the soil already contained 10% gravimetric water, i.e., ~0.9 L). Changes in total mesocosm weight hence allowed for the estimation of the SWC during the experiment and its expression in terms of % of field capacity. At field capacity, the mesocosm water-filled pore space (WFPS) can be estimated at 71.0 ± 2.5% on average (calculated as WFPS = water content/porosity; porosity = 1 - bulk density/particle density; particle density = 2.7 Mg m⁻³). The mesocosms were exposed to a drying phase until one of them reached almost 50% of field capacity, which occurred after 6 weeks. The volume of water lost by the driest mesocosm was determined and then added (second watering in two sessions again) to all mesocosms to set them back to 100% of field capacity. A second drying phase was imposed until the penultimate week, when a third watering was performed with the same amount of water as supplied in the second watering. Following this method, all mesocosms experienced two drying-rewetting cycles (Fig 1A and 1B).

## Response variables

CO₂ and N₂O emissions measurements were carried out weekly and during the weeks with watering, always 24 h after watering events, to capture eventual emission peaks. A static sampling chamber approach was used, following the recommendations of Rochette (2011) [50]. The sampling chamber (9 cm dia., 6 cm height, 370 ± 1 mL, S1 Fig) was equipped with a bung of silicone rubber for gas sampling at the top. During sampling, the chamber was placed on a circular collar (S1 Fig) that was inserted in the center of the mesocosms during the setup of the mesocosms, which allowed measuring soil N₂O and CO₂ fluxes without disturbing the plants.
earthworm FW biomasses were the equivalent of 480, 145 and 330 g m⁻² for the anecic, endogeic and both earthworm treatment levels, respectively, which were 2- to 3-fold higher than that in the field of origin of the soil, where earthworm total biomass varied from 98 to 135 g m⁻². For *L. terrestris*, adding the proportional equivalent of the field biomass would mean introducing just one individual. Given the associated risk that the escape or death of this single individual could jeopardize the treatment, we chose to add two adult individuals of *L. terrestris*. As a result, the overall biomass was roughly double what is typically found per m^2 in agroecosystems. Note that *L. terrestris* was more recently reclassified as a species displaying traits belonging to both anecic and epigeic species [36], but in the vast majority of literature, it is considered an anecic species. The earthworm treatments were factorially crossed with two levels of plant (*B. distachyon* L.) treatments, i.e., with and without plants, with 7 replicates per treatment combination for a total of 56 mesocosms.

The experiment was conducted in a greenhouse kept at temperatures ranging between 20 and 23°C during the day and 18–20°C during the night with an air relative humidity of 80%. The natural light was supplemented during daytime by artificial lighting for 12 hours per day using high-pressure sodium lamps. The mesocosms were divided into two blocks corresponding to their position on the north or south bench of the greenhouse, and their position within the block was randomly changed twice a week, limiting the bias of the position within the block. The experiment ran for 12 weeks between March and May 2017.

## Watering protocol

The watering protocol was specifically designed to include soil moisture fluctuations (analogous to what occurs in natural conditions) and to allow earthworm burrowing to affect SWC. At the beginning of the experiment, the mesocosms were watered with 1.7 L of reverse osmosis water using a laboratory dispenser (two sessions of 850 ml each), a volume sufficient to observe water draining out of the mesocosms from the pierced lids at the bottoms of mesocosms. Measurements of weight changes after 24 h were used to calculate the weights of the mesocosms at field capacity (knowing that the soil already contained 10% gravimetric water, i.e., ~0.9 L). Changes in total mesocosm weight hence allowed for the estimation of the SWC during the experiment and its expression in terms of % of field capacity. At field capacity, the mesocosm water-filled pore space (WFPS) can be estimated at 71.0 ± 2.5% on average (calculated as WFPS = water content/porosity; porosity = 1 - bulk density/particle density; particle density = 2.7 Mg m⁻³). The mesocosms were exposed to a drying phase until one of them reached almost 50% of field capacity, which occurred after 6 weeks. The volume of water lost by the driest mesocosm was determined and then added (second watering in two sessions again) to all mesocosms to set them back to 100% of field capacity. A second drying phase was imposed until the penultimate week, when a third watering was performed with the same amount of water as supplied in the second watering. Following this method, all mesocosms experienced two drying-rewetting cycles (Fig 1A and 1B).

## Response variables

CO₂ and N₂O emissions measurements were carried out weekly and during the weeks with watering, always 24 h after watering events, to capture eventual emission peaks. A static sampling chamber approach was used, following the recommendations of Rochette (2011) [50]. The sampling chamber (9 cm dia., 6 cm height, 370 ± 1 mL, S1 Fig) was equipped with a bung of silicone rubber for gas sampling at the top. During sampling, the chamber was placed on a circular collar (S1 Fig) that was inserted in the center of the mesocosms during the setup of the mesocosms, which allowed measuring soil N₂O and CO₂ fluxes without disturbing the plants.
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Temporal dynamics of soil water content and litter cover. (A,B) Temporal dynamics of soil water content (SWC) expressed as percentage of field capacity (C) and percentage of surface covered by litter as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatments independent of plant treatments, and vice versa, are displayed when the Ew×Plant interaction is not significant (A,B).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g001

The collar was inserted into the soil down to 3 cm and consisted of a frame that provided support but allowed access by earthworms and roots to the inner soil core thanks to two open windows. The aboveground part of the collar contained a gutter-like double-walled section/groove where the static chamber was placed during sampling. On the day of gas sampling, the
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Temporal dynamics of soil water content and litter cover. (A,B) Temporal dynamics of soil water content (SWC) expressed as percentage of field capacity (C) and percentage of surface covered by litter as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatments independent of plant treatments, and vice versa, are displayed when the Ew×Plant interaction is not significant (A,B).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g001

The collar was inserted into the soil down to 3 cm and consisted of a frame that provided support but allowed access by earthworms and roots to the inner soil core thanks to two open windows. The aboveground part of the collar contained a gutter-like double-walled section/groove where the static chamber was placed during sampling. On the day of gas sampling, the
mesocosms were transported to the gas sampling laboratory by trolley, and the two blocks were measured over two days. Prior to sampling, 20 ml of distilled water was added to the groove to provide airtight sealing when placing the static chamber on the collar. CO₂ and N₂O fluxes were measured at the Platform for Chemical Analysis in Ecology (LabEx CeMEB, Montpellier, France). CO₂ concentrations were measured with a gas chromatograph (MicroGC S-Series, SRA XXX Instruments, Marcy l'Etoile, France) using a catharometric detector, quantifying the gases on the basis of their thermal conductivity. N₂O concentrations were measured by a gas chromatograph equipped with an electron capture detector (Varian CP-3800, Varian Inc., Palo Alto, USA). Air samples were taken at T0 (immediately after placing the chamber on the collar) and after 2 h to assess the changes in CO₂ and N₂O concentrations. Previous tests were conducted after 1, 2, 3 and 4 hours, which revealed that gas accumulation was linear during this short time period. A volume of 0.2 mL was sequentially sampled for gas measurements via the silicone bung using a plastic syringe equipped with a 25G needle and was injected immediately into the gas chromatograph via a 1/32"PFA line. Concentration changes in the sampling chamber between T0 and T0+2 h were used to estimate the greenhouse gas emission rates and converted to g C-CO₂ and N-N₂O m⁻² day⁻¹.

At the end of the experiment, the mesocosms were transported to the INRAE center of Nancy to analyze soil macroporosity by X-ray tomography using a medical scanner (Bright-Speed Exel 4, General Electric), with settings of 120 kV and 50 mA for the current and 0.625 mm width for each image. Images were transformed into 16-bit images and binarized (i.e., converted into black and white) using a fixed threshold value [51] because the different peaks (for the soil matrix and the porosity) were well separated [52]. Roots and associated pores could not be included in the analysis due to their smaller average size compared to the resolution of the scanner (0.4 mm per pixel). The burrow system was then characterized by computing the volume and number of burrows in four soil layers (L1 for 0–8.5 cm, L2 for 8.5–17 cm, L3 for 17–25.5 cm and L4 for 25.5–34 cm depth, S2 Fig) using ImageJ [41]. Drying–wetting cycles contributed to the formation of cracks, i.e., macropores resulting from physical processes (shrinkage, swelling [41]) notably in the topsoil layer (S2 Fig). We attempted to differentiate cracks from burrows according to the macropore circularity in 2D images because earthworm burrows are more circular than cracks. However, this method still identified burrows in the control mesocosms (without earthworms), with 25.4% of pores misidentified as burrows in the whole mesocosm and higher error in the first layer (32.2%) than in the bottom layer L4 (20.1%) (S2 Fig). Therefore, assuming that gas fluxes were influenced by the total porosity, regardless of its biological or physical origin, and due to the high correlations between the 15 different porosity variables (S3 Fig), although we investigated how the treatment affected the different types of soil pores (pores, burrows, and cracks), we decided to only use the total pore volume data (burrows and cracks) as a predictor in the models.

After the X-ray scan in Nancy, the mesocosms were transported back to Montpellier for the final destructive harvest. The proportions of earthworms found at the final harvest were 62% and 90% for *L. terrestris* and *A. icterica*, respectively. The proportions of recovered earthworms were likely affected by the mortality that occurred during the days of transport and storage for the X-ray scans (during which the temperatures and vibrations were not controlled), as several *L. terrestris* individuals were found freshly dead at harvest. However, the X-ray scans together with the litter mass loss dynamics (Fig 1C) provided strong evidence that the earthworms were active during the entire duration of the experiment. Litter cover was assessed weekly by a non-destructive visual estimation method with 5% intervals, as commonly performed for ground cover estimation [53, 54]. The observer bias of this method [55] was handled by having only one observer perform the estimations throughout the entire experiment. Other additional soil- and plant-related response variables were measured at the end of the experiment (S3 Fig). Soil
mesocosms were transported to the gas sampling laboratory by trolley, and the two blocks were measured over two days. Prior to sampling, 20 ml of distilled water was added to the groove to provide airtight sealing when placing the static chamber on the collar. CO₂ and N₂O fluxes were measured at the Platform for Chemical Analysis in Ecology (LabEx CeMEB, Montpellier, France). CO₂ concentrations were measured with a gas chromatograph (MicroGC S-Series, SRA XXX Instruments, Marcy l'Etoile, France) using a catharometric detector, quantifying the gases on the basis of their thermal conductivity. N₂O concentrations were measured by a gas chromatograph equipped with an electron capture detector (Varian CP-3800, Varian Inc., Palo Alto, USA). Air samples were taken at T0 (immediately after placing the chamber on the collar) and after 2 h to assess the changes in CO₂ and N₂O concentrations. Previous tests were conducted after 1, 2, 3 and 4 hours, which revealed that gas accumulation was linear during this short time period. A volume of 0.2 mL was sequentially sampled for gas measurements via the silicone bung using a plastic syringe equipped with a 25G needle and was injected immediately into the gas chromatograph via a 1/32"PFA line. Concentration changes in the sampling chamber between T0 and T0+2 h were used to estimate the greenhouse gas emission rates and converted to g C-CO₂ and N-N₂O m⁻² day⁻¹.

At the end of the experiment, the mesocosms were transported to the INRAE center of Nancy to analyze soil macroporosity by X-ray tomography using a medical scanner (Bright-Speed Exel 4, General Electric), with settings of 120 kV and 50 mA for the current and 0.625 mm width for each image. Images were transformed into 16-bit images and binarized (i.e., converted into black and white) using a fixed threshold value [51] because the different peaks (for the soil matrix and the porosity) were well separated [52]. Roots and associated pores could not be included in the analysis due to their smaller average size compared to the resolution of the scanner (0.4 mm per pixel). The burrow system was then characterized by computing the volume and number of burrows in four soil layers (L1 for 0–8.5 cm, L2 for 8.5–17 cm, L3 for 17–25.5 cm and L4 for 25.5–34 cm depth, S2 Fig) using ImageJ [41]. Drying–wetting cycles contributed to the formation of cracks, i.e., macropores resulting from physical processes (shrinkage, swelling [41]) notably in the topsoil layer (S2 Fig). We attempted to differentiate cracks from burrows according to the macropore circularity in 2D images because earthworm burrows are more circular than cracks. However, this method still identified burrows in the control mesocosms (without earthworms), with 25.4% of pores misidentified as burrows in the whole mesocosm and higher error in the first layer (32.2%) than in the bottom layer L4 (20.1%) (S2 Fig). Therefore, assuming that gas fluxes were influenced by the total porosity, regardless of its biological or physical origin, and due to the high correlations between the 15 different porosity variables (S3 Fig), although we investigated how the treatment affected the different types of soil pores (pores, burrows, and cracks), we decided to only use the total pore volume data (burrows and cracks) as a predictor in the models.

After the X-ray scan in Nancy, the mesocosms were transported back to Montpellier for the final destructive harvest. The proportions of earthworms found at the final harvest were 62% and 90% for *L. terrestris* and *A. icterica*, respectively. The proportions of recovered earthworms were likely affected by the mortality that occurred during the days of transport and storage for the X-ray scans (during which the temperatures and vibrations were not controlled), as several *L. terrestris* individuals were found freshly dead at harvest. However, the X-ray scans together with the litter mass loss dynamics (Fig 1C) provided strong evidence that the earthworms were active during the entire duration of the experiment. Litter cover was assessed weekly by a non-destructive visual estimation method with 5% intervals, as commonly performed for ground cover estimation [53, 54]. The observer bias of this method [55] was handled by having only one observer perform the estimations throughout the entire experiment. Other additional soil- and plant-related response variables were measured at the end of the experiment (S3 Fig). Soil
analyses were performed on a homogenized soil sample from the upper 10 cm of the meso-
cosms inside the collar and sieved at 2 mm. Potential soil microbial denitrification enzymatic
activity (DEA) was measured using the acetylene inhibition method, which measures total
potential denitrification (as N₂O and N₂) [56]. This is a complementary method to the fluxes
measured during the experiment, which only measure the N₂O emissions. The MicroResp™
method was used to determine the microbial metabolic quotient [57]. Approximately 0.39 g
dry weight of soil was incubated in six replication wells with a solution of D-glucose (1.5 mg C
g⁻¹ soil) and six replication wells with deionized water (for basal respiration) to reach 80% of
the field capacity in 96-DeepWell Microplates (Fisher Scientific E39199). Cresol red gel detec-
tion plates were prepared as recommended by the manufacturer. After an initial two-hour pre-
incubation at 25°C in the dark, each DeepWell microplate was covered with a CO₂-trap
microplate detection plate using a silicone gasket (MicroResp™, Aberdeen, UK). The assembly
was secured with a clamp and incubated for four additional hours. The optical density at 590
nm (OD590) was measured for each detection well before and after incubation using a Victor
1420 Multilabel Counter (Perkin Elmer, Massachusetts, USA). Calibration relying on absor-
bance (OD590 readings) and CO₂ concentrations was performed using the gas chromatograph
previously described. The final OD590 values were normalized using preincubation OD590
and converted as respiration rates expressed in µg C-CO₂ respired per g⁻¹ of soil per h⁻¹. The
glucose-induced respiration rate was used to estimate the soil microbial C (Cmic, µg Cmicrobial
g⁻¹ dry soil) biomass [58]. Finally, the metabolic quotient (Met_Q) was determined as the ratio
between basal respiration rates measured in the wells with water only and no C substrate (as a
proxy of the microbial basal respiration) and Cmic. Soil mineral nitrogen was extracted from
10 g of freshly sampled soil with 40 mL of 1 M KCl solution. Nitrate and ammonium concen-
trations were measured by continuous flow spectrophotometry (SKALAR 3000 auto analyzer,
Breda, The Netherlands). The plant shoot biomass was weighed after drying at 60°C for three
days. As the roots of *B. dystachyon* are extremely thin and fragile, it was not feasible to sample
root biomass.

Statistical analyses

Statistical analyses were performed using R version 4.0.2 (R Development Core Team, 2015) in RStudio version 1.3.959 (RStudio Team, 2015). Weekly time series of CO₂ and N₂O emissions and SWC were analyzed with the “nmle” package version 3.1–145 [59] to perform repeated measures analyses using a generalized mixed-effects model to test the effects of earthworms, plants, SWC (for GHG emissions only) and sampling week and their interactions on gas fluxes. The identity (ID) of the mesocosm and its position in the blocks were used as random factors to account for temporal pseudoreplication and the effect of the position in the north or south bench in the greenhouse (“random = ~ 1 | Block / ID”). To reach models that respected the assumption of homoscedasticity of the residuals, we tested the model fit with varIdent (for plant and earthworm experimental treatments), varPower (for SWC) and varExp (for SWC) weighting functions [60] and selected the most appropriate models based on maximum likelihood (ML) model comparison tests. A similar approach was used for the cumulative CO₂ and N₂O data at the end of the experiment (estimated assuming constant emission rates between the weekly measurements), but without the sampling week among the fixed effects and the mesocosm ID in the random effects. For the later analysis, we used the mean of weekly SWC values, as this variable is arguably more relevant to the cumulative fluxes. The “r.squaredGLMM” function from the MuMln package [61] was used to derive the proportion of the variation explained by the fixed factors (i.e., marginal r², mr²) in mixed-effects models. The “multcom” package was used to perform Tukey’s HSD (honestly significant difference)
analyses were performed on a homogenized soil sample from the upper 10 cm of the meso-
cosms inside the collar and sieved at 2 mm. Potential soil microbial denitrification enzymatic
activity (DEA) was measured using the acetylene inhibition method, which measures total
potential denitrification (as N₂O and N₂) [56]. This is a complementary method to the fluxes
measured during the experiment, which only measure the N₂O emissions. The MicroResp™
method was used to determine the microbial metabolic quotient [57]. Approximately 0.39 g
dry weight of soil was incubated in six replication wells with a solution of D-glucose (1.5 mg C
g⁻¹ soil) and six replication wells with deionized water (for basal respiration) to reach 80% of
the field capacity in 96-DeepWell Microplates (Fisher Scientific E39199). Cresol red gel detec-
tion plates were prepared as recommended by the manufacturer. After an initial two-hour pre-
incubation at 25°C in the dark, each DeepWell microplate was covered with a CO₂-trap
microplate detection plate using a silicone gasket (MicroResp™, Aberdeen, UK). The assembly
was secured with a clamp and incubated for four additional hours. The optical density at 590
nm (OD590) was measured for each detection well before and after incubation using a Victor
1420 Multilabel Counter (Perkin Elmer, Massachusetts, USA). Calibration relying on absor-
bance (OD590 readings) and CO₂ concentrations was performed using the gas chromatograph
previously described. The final OD590 values were normalized using preincubation OD590
and converted as respiration rates expressed in µg C-CO₂ respired per g⁻¹ of soil per h⁻¹. The
glucose-induced respiration rate was used to estimate the soil microbial C (Cmic, µg Cmicrobial
g⁻¹ dry soil) biomass [58]. Finally, the metabolic quotient (Met_Q) was determined as the ratio
between basal respiration rates measured in the wells with water only and no C substrate (as a
proxy of the microbial basal respiration) and Cmic. Soil mineral nitrogen was extracted from
10 g of freshly sampled soil with 40 mL of 1 M KCl solution. Nitrate and ammonium concen-
trations were measured by continuous flow spectrophotometry (SKALAR 3000 auto analyzer,
Breda, The Netherlands). The plant shoot biomass was weighed after drying at 60°C for three
days. As the roots of *B. dystachyon* are extremely thin and fragile, it was not feasible to sample
root biomass.

Statistical analyses

Statistical analyses were performed using R version 4.0.2 (R Development Core Team, 2015) in RStudio version 1.3.959 (RStudio Team, 2015). Weekly time series of CO₂ and N₂O emissions and SWC were analyzed with the “nmle” package version 3.1–145 [59] to perform repeated measures analyses using a generalized mixed-effects model to test the effects of earthworms, plants, SWC (for GHG emissions only) and sampling week and their interactions on gas fluxes. The identity (ID) of the mesocosm and its position in the blocks were used as random factors to account for temporal pseudoreplication and the effect of the position in the north or south bench in the greenhouse (“random = ~ 1 | Block / ID”). To reach models that respected the assumption of homoscedasticity of the residuals, we tested the model fit with varIdent (for plant and earthworm experimental treatments), varPower (for SWC) and varExp (for SWC) weighting functions [60] and selected the most appropriate models based on maximum likelihood (ML) model comparison tests. A similar approach was used for the cumulative CO₂ and N₂O data at the end of the experiment (estimated assuming constant emission rates between the weekly measurements), but without the sampling week among the fixed effects and the mesocosm ID in the random effects. For the later analysis, we used the mean of weekly SWC values, as this variable is arguably more relevant to the cumulative fluxes. The “r.squaredGLMM” function from the MuMln package [61] was used to derive the proportion of the variation explained by the fixed factors (i.e., marginal r², mr²) in mixed-effects models. The “multcom” package was used to perform Tukey’s HSD (honestly significant difference)
multicomparison posthoc test, but note that this test does not include the random effects, and
the results are occasionally not entirely in line with the fitted coefficients from the mixed-
effects models.

Additional analyses were conducted aiming to link the multiple potential predictors (S3 Fig) measured at the end of the experiment and the CO2 and N2O fluxes from week 12 (just before the experiment was stopped). As multiple response variables were measured to explore potential predictors, a method of best subset selection that penalizes model complexity (i.e., regularization) during estimation was needed. Regularization aims to significantly reduce the variance of the model as well as model overfitting by varying the lambda (λ) parameter, which tunes the level of penalization for the complexity of the model. This approach has been proven to be a viable option for estimating parameters in scenarios with small sample sizes and many collinear/correlated predictors. Here, we used a penalized regression method based on the minimax concave penalty (MCP) to select the best subsets [62] using the *ncvreg* package 3.11-1 [63]. This approach was combined with a 10-fold cross-validation procedure to derive the lambda parameter (also called the regularization rate), which minimizes the cross-validation error. We report the fitted coefficients and the coefficients of determination (r²) at lambda values that minimize the cross-validation error. The subset variables with retained nonzero coefficients were then tested in the generalized mixed-effects models, which have the advantage of including random effects alongside the treatment factors.

## Results

### Soil water content

Over time, the SWC, expressed as a percentage of field capacity, exhibited variations due to the drying-rewetting cycles and the treatments involving plants and earthworms, which were reflected in the Plant×Week and Ew×Week interactions (Table 1 and Fig 1A and 1B). The SWC was significantly lower in the presence of plants during the last three weeks of both drying cycles, with 5% lower values in the absence of plants across the whole experiment (averaged over the earthworm treatments). The presence of anecic earthworms (with or without endogeic earthworms) led to significantly lower SWC values relative to the control during weeks three to six and eight to twelve. Averaged over the twelve weeks, the SWC values were 81.3% of field capacity in the presence of anecic earthworms, 81.7% in the presence of both earthworm species, 85.9% with endogeic earthworms and 87.0% in the control.

### Weekly and cumulative N₂O and CO₂ fluxes

Weekly N₂O emissions were significantly affected by all possible three-way interactions between earthworms, plants, SWC and time (Table 1 and Fig 2A). N₂O emissions were higher after each watering (with the measurements always being performed 24 h after watering) and were the highest in the second week of the experiment (0.10 g N-N₂O m⁻² day⁻¹, Fig 2A). The intensity and duration of these emission peaks depended on the earthworm and plant treatments but varied with time, as indicated by the significant Ew×Plant×Week interaction. For example, N₂O emissions peaked ten days after the first watering in the presence of anecic earthworms compared to the presence of endogeic earthworms and reached a maximum in the absence of plants. Cumulative N₂O emissions over the 12 weeks of the experiment were the highest in the control without earthworms or plants (0.135 g m⁻², Fig 2B and 2C). Relative to the control, the cumulative N₂O emissions were 17.0%, 34.6% and 44.8% lower in the anecic, both and endogeic treatments, respectively, and 19.8% lower in the presence of plants (Fig 2C). The Ew×SWC interaction indicates that cumulative N₂O emissions increased with average SWC for the control, anecic and both earthworm treatments but decreased with SWC
multicomparison posthoc test, but note that this test does not include the random effects, and
the results are occasionally not entirely in line with the fitted coefficients from the mixed-
effects models.

Additional analyses were conducted aiming to link the multiple potential predictors (S3 Fig) measured at the end of the experiment and the CO2 and N2O fluxes from week 12 (just before the experiment was stopped). As multiple response variables were measured to explore potential predictors, a method of best subset selection that penalizes model complexity (i.e., regularization) during estimation was needed. Regularization aims to significantly reduce the variance of the model as well as model overfitting by varying the lambda (λ) parameter, which tunes the level of penalization for the complexity of the model. This approach has been proven to be a viable option for estimating parameters in scenarios with small sample sizes and many collinear/correlated predictors. Here, we used a penalized regression method based on the minimax concave penalty (MCP) to select the best subsets [62] using the *ncvreg* package 3.11-1 [63]. This approach was combined with a 10-fold cross-validation procedure to derive the lambda parameter (also called the regularization rate), which minimizes the cross-validation error. We report the fitted coefficients and the coefficients of determination (r²) at lambda values that minimize the cross-validation error. The subset variables with retained nonzero coefficients were then tested in the generalized mixed-effects models, which have the advantage of including random effects alongside the treatment factors.

Results

Soil water content

Over time, the SWC, expressed as a percentage of field capacity, exhibited variations due to the drying-rewetting cycles and the treatments involving plants and earthworms, which were reflected in the Plant×Week and Ew×Week interactions (Table 1 and Fig 1A and 1B). The SWC was significantly lower in the presence of plants during the last three weeks of both drying cycles, with 5% lower values in the absence of plants across the whole experiment (averaged over the earthworm treatments). The presence of anecic earthworms (with or without endogeic earthworms) led to significantly lower SWC values relative to the control during weeks three to six and eight to twelve. Averaged over the twelve weeks, the SWC values were 81.3% of field capacity in the presence of anecic earthworms, 81.7% in the presence of both earthworm species, 85.9% with endogeic earthworms and 87.0% in the control.

Weekly and cumulative N₂O and CO₂ fluxes

Weekly N₂O emissions were significantly affected by all possible three-way interactions
between earthworms, plants, SWC and time (Table 1 and Fig 2A). N₂O emissions were higher
after each watering (with the measurements always being performed 24 h after watering) and
were the highest in the second week of the experiment (0.10 g N-N₂O m⁻² day⁻¹, Fig 2A). The
intensity and duration of these emission peaks depended on the earthworm and plant treat-
ments but varied with time, as indicated by the significant Ew×Plant×Week interaction. For
example, N₂O emissions peaked ten days after the first watering in the presence of anecic
earthworms compared to the presence of endogeic earthworms and reached a maximum in
the absence of plants. Cumulative N₂O emissions over the 12 weeks of the experiment were
the highest in the control without earthworms or plants (0.135 g m⁻², Fig 2B and 2C). Relative
to the control, the cumulative N₂O emissions were 17.0%, 34.6% and 44.8% lower in the
anecic, both and endogeic treatments, respectively, and 19.8% lower in the presence of plants
(Fig 2C). The Ew×SWC interaction indicates that cumulative N₂O emissions increased with
average SWC for the control, anecic and both earthworm treatments but decreased with SWC
<table><caption>Table 1. Time series and cumulative statistic table.</caption><thead><tr><th>Source</th><th>SWC</th><th>Litter</th><th>N<sub>2</sub>O</th><th>cN<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>cCO<sub>2</sub></th></tr></thead><tbody><tr><td>Ew</td><td>46.68***</td><td>182.89***</td><td>82.97***</td><td>7.07***</td><td>2.53<sup>+</sup></td><td>1.25</td></tr><tr><td>Plant</td><td>173.83***</td><td>2.91<sup>+</sup></td><td>34.58***</td><td>3.82<sup>+</sup></td><td>5.98<sup>*</sup></td><td>ns</td></tr><tr><td>Week</td><td>2758.70***</td><td>68.85***</td><td>112.9***</td><td>NA</td><td>230.28***</td><td>NA</td></tr><tr><td>SWC</td><td>NA</td><td>NA</td><td>0.04</td><td>6.35<sup>*</sup></td><td>128.85***</td><td>7.54<sup>**</sup></td></tr><tr><td>Ew×Plant</td><td>2.63<sup>+</sup></td><td>3.65<sup>*</sup></td><td>9.52***</td><td>2.05</td><td>2.01</td><td>ns</td></tr><tr><td>Ew×Week</td><td>9.40***</td><td>198.63***</td><td>5.61***</td><td>NA</td><td>1.97<sup>**</sup></td><td>NA</td></tr><tr><td>Plant×Week</td><td>24.01***</td><td>ns</td><td>9.32***</td><td>NA</td><td>4.93***</td><td>NA</td></tr><tr><td>Ew×SWC</td><td>NA</td><td>NA</td><td>14.84***</td><td>4.45<sup>**</sup></td><td>5.16<sup>**</sup></td><td>3.16<sup>*</sup></td></tr><tr><td>Plant×SWC</td><td>NA</td><td>NA</td><td>1.43</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td>Week×SWC</td><td>NA</td><td>NA</td><td>8.34***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×Week</td><td>ns</td><td>ns</td><td>1.57<sup>*</sup></td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×SWC</td><td>NA</td><td>NA</td><td>9.51***</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td>Ew×Week×SWC</td><td>NA</td><td>NA</td><td>3.81***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Plant×Week×SWC</td><td>NA</td><td>NA</td><td>5.52***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×Week×SWC</td><td>NA</td><td>NA</td><td>ns</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>m<sup>2</sup>r<sup>2</sup></td><td>0.81</td><td>0.91</td><td>0.67</td><td>0.49</td><td>0.43</td><td>0.27</td></tr></tbody></table>

Minimal adequate models for weekly time series (SWC, Litter cover, N₂O, and CO₂) and cumulative emissions (cN₂O and cCO₂) as affected by the earthworm (Ew), plant (Plant), sampling week (Week), soil water content (SWC) and their interactions. “NA” stands for non-applicable, “ns” stands for variables that were not significant (P > 0.1) and were not retained in the minimal adequate models whereas m² represents the marginal coefficient of determination. Figures are F-values.

***P < 0.001

**P < 0.01

*P < 0.05

+P < 0.1.

https://doi.org/10.1371/journal.pone.0289859.t001

in the presence of endogeic earthworms (Table 1 and Fig 4A). CO₂ weekly emissions were significantly affected by the Ew×Week, Plant×Week and Ew×SWC interactions and were higher after rewatering at higher SWC (Table 1 and Fig 3A). The cumulative CO₂ emissions showed no significant response to earthworm or plant presence (Fig 3B and 3C). However, a significant interactive effect of earthworms with SWC was found (Table 1 and Fig 4E), indicating that when SWC was relatively high (> 85.9% SWC on average), as in the endogeic and earthworm-control treatments, cumulative CO₂ emissions generally decreased with increasing SWC. The opposite was true in the presence of anecic species, while no relationship was found in the presence of both species. Additionally, when the SWC effect was not included as a predictor, the cumulative CO₂ emissions relative to the control were 5.9% and 11.4% lower in the endogeic and both earthworm treatments, respectively, as indicated by Tukey's HSD post hoc test (Fig 3B). Conversely, plant presence increased the cumulative CO₂ emissions by 6%.

## Plant, litter, microbial activity, nutrients and porosity

Although not significant, the presence of earthworms led to higher aboveground plant biomass compared to the control (+57%, +25% and +41% for anecic, endogeic and both species, respectively, Table 2 and S4A Fig). Litter cover depended on the presence of earthworms and varied with time and plant presence, as indicated by the Ew×Week and Ew×Plant interactions, respectively (Table 1 and Fig 1C). After 12 weeks, anecic earthworms reduced litter cover to 5% and 6.4% in the absence and presence of plants, respectively. Moreover, when both species were present, litter cover was 40% without plants and 15% with plants (S4B Fig). At week 12, SWC was positively correlated with litter cover (r = 0.55, t = 4.86, P < 0.001, n = 56; S3 Fig),
<table><caption>Table 1. Time series and cumulative statistic table.</caption><thead><tr><th>Source</th><th>SWC</th><th>Litter</th><th>N<sub>2</sub>O</th><th>cN<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>cCO<sub>2</sub></th></tr></thead><tbody><tr><td>Ew</td><td>46.68***</td><td>182.89***</td><td>82.97***</td><td>7.07***</td><td>2.53<sup>+</sup></td><td>1.25</td></tr><tr><td>Plant</td><td>173.83***</td><td>2.91<sup>+</sup></td><td>34.58***</td><td>3.82<sup>+</sup></td><td>5.98<sup>*</sup></td><td>ns</td></tr><tr><td>Week</td><td>2758.70***</td><td>68.85***</td><td>112.9***</td><td>NA</td><td>230.28***</td><td>NA</td></tr><tr><td>SWC</td><td>NA</td><td>NA</td><td>0.04</td><td>6.35<sup>*</sup></td><td>128.85***</td><td>7.54<sup>**</sup></td></tr><tr><td>Ew×Plant</td><td>2.63<sup>+</sup></td><td>3.65<sup>*</sup></td><td>9.52***</td><td>2.05</td><td>2.01</td><td>ns</td></tr><tr><td>Ew×Week</td><td>9.40***</td><td>198.63***</td><td>5.61***</td><td>NA</td><td>1.97<sup>**</sup></td><td>NA</td></tr><tr><td>Plant×Week</td><td>24.01***</td><td>ns</td><td>9.32***</td><td>NA</td><td>4.93***</td><td>NA</td></tr><tr><td>Ew×SWC</td><td>NA</td><td>NA</td><td>14.84***</td><td>4.45<sup>**</sup></td><td>5.16<sup>**</sup></td><td>3.16<sup>*</sup></td></tr><tr><td>Plant×SWC</td><td>NA</td><td>NA</td><td>1.43</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td>Week×SWC</td><td>NA</td><td>NA</td><td>8.34***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×Week</td><td>ns</td><td>ns</td><td>1.57<sup>*</sup></td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×SWC</td><td>NA</td><td>NA</td><td>9.51***</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td>Ew×Week×SWC</td><td>NA</td><td>NA</td><td>3.81***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Plant×Week×SWC</td><td>NA</td><td>NA</td><td>5.52***</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>Ew×Plant×Week×SWC</td><td>NA</td><td>NA</td><td>ns</td><td>NA</td><td>ns</td><td>NA</td></tr><tr><td>m<sup>2</sup>r<sup>2</sup></td><td>0.81</td><td>0.91</td><td>0.67</td><td>0.49</td><td>0.43</td><td>0.27</td></tr></tbody></table>

Minimal adequate models for weekly time series (SWC, Litter cover, N₂O, and CO₂) and cumulative emissions (cN₂O and cCO₂) as affected by the earthworm (Ew), plant (Plant), sampling week (Week), soil water content (SWC) and their interactions. “NA” stands for non-applicable, “ns” stands for variables that were not significant (P > 0.1) and were not retained in the minimal adequate models whereas m² represents the marginal coefficient of determination. Figures are F-values.

***P < 0.001

**P < 0.01

*P < 0.05

+P < 0.1.

https://doi.org/10.1371/journal.pone.0289859.t001

in the presence of endogeic earthworms (Table 1 and Fig 4A). CO₂ weekly emissions were significantly affected by the Ew×Week, Plant×Week and Ew×SWC interactions and were higher after rewatering at higher SWC (Table 1 and Fig 3A). The cumulative CO₂ emissions showed no significant response to earthworm or plant presence (Fig 3B and 3C). However, a significant interactive effect of earthworms with SWC was found (Table 1 and Fig 4E), indicating that when SWC was relatively high (> 85.9% SWC on average), as in the endogeic and earthworm-control treatments, cumulative CO₂ emissions generally decreased with increasing SWC. The opposite was true in the presence of anecic species, while no relationship was found in the presence of both species. Additionally, when the SWC effect was not included as a predictor, the cumulative CO₂ emissions relative to the control were 5.9% and 11.4% lower in the endogeic and both earthworm treatments, respectively, as indicated by Tukey's HSD post hoc test (Fig 3B). Conversely, plant presence increased the cumulative CO₂ emissions by 6%.

## Plant, litter, microbial activity, nutrients and porosity

Although not significant, the presence of earthworms led to higher aboveground plant biomass compared to the control (+57%, +25% and +41% for anecic, endogeic and both species, respectively, Table 2 and S4A Fig). Litter cover depended on the presence of earthworms and varied with time and plant presence, as indicated by the Ew×Week and Ew×Plant interactions, respectively (Table 1 and Fig 1C). After 12 weeks, anecic earthworms reduced litter cover to 5% and 6.4% in the absence and presence of plants, respectively. Moreover, when both species were present, litter cover was 40% without plants and 15% with plants (S4B Fig). At week 12, SWC was positively correlated with litter cover (r = 0.55, t = 4.86, P < 0.001, n = 56; S3 Fig),
<figure><img src="image_6.png" /><figcaption>**Fig 2.** N₂O emissions during the 12-week experiment. (A) Weekly (B, C) and cumulative N–N₂O emissions as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatment independent of plant treatment, and vice-versa, are displayed when the Ew×Plant interaction is not significant (B, C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g002

suggesting that, at least in part, the presence of earthworms (especially the anecic L. terrestris) also affected SWC via higher evaporation from bare soil due to litter burial. We found marginally higher denitrification potential in the presence of anecic earthworms (+22%) relative to the control but no plant effects (Table 2 and S4C Fig). The soil nitrate content (NO₃⁻) was
<figure><img src="image_6.png" /><figcaption>**Fig 2.** N₂O emissions during the 12-week experiment. (A) Weekly (B, C) and cumulative N–N₂O emissions as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatment independent of plant treatment, and vice-versa, are displayed when the Ew×Plant interaction is not significant (B, C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g002

suggesting that, at least in part, the presence of earthworms (especially the anecic L. terrestris) also affected SWC via higher evaporation from bare soil due to litter burial. We found marginally higher denitrification potential in the presence of anecic earthworms (+22%) relative to the control but no plant effects (Table 2 and S4C Fig). The soil nitrate content (NO₃⁻) was
<figure><img src="image_7.png" /><figcaption>**Fig 3. CO₂ emissions during the 12-week experiment.** (A) Weekly (B, C) and cumulative C-CO₂ emissions as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatment independent of plant treatment, and vice-versa, are displayed when the Ew×Plant interaction is not significant (B, C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g003

always lower (-80%) in the presence of plants and increased in the presence of earthworms, with synergistic effects in the presence of both species (Table 2 and S4G Fig). The macropore volume in the topsoil layer (L1) was not affected by any experimental treatment (Table 3 and S2 Fig). We found lower macropore volume in the presence of plants in the other three layers
A

<figure><img src="image_7.png" /></figure>

B

<figure><img src="image_8.png" /></figure>

C

<figure><img src="image_9.png" /></figure>

**Fig 3. CO₂ emissions during the 12-week experiment.** (A) Weekly (B, C) and cumulative C-CO₂ emissions as affected by the earthworm and plant treatments. Blue arrows represent watering events. Error bars represent ± 1 SEM. Different letters represent significantly different levels as estimated by Tukey's HSD post hoc test. The effects of earthworm treatment independent of plant treatment, and vice-versa, are displayed when the Ew×Plant interaction is not significant (B, C).

https://doi.org/10.1371/journal.pone.0289859.g003

always lower (-80%) in the presence of plants and increased in the presence of earthworms, with synergistic effects in the presence of both species (Table 2 and S4G Fig). The macropore volume in the topsoil layer (L1) was not affected by any experimental treatment (Table 3 and S2 Fig). We found lower macropore volume in the presence of plants in the other three layers
<figure><img src="image_8.png" /><figcaption>**Fig 4.** N₂O and CO₂ emissions as affected by soil water content and macroporosity. (A, E) Total cumulative emissions as affected by the EwxSWC interaction, where SWC represents the 3-month average SWC. (B, C, D, F, G, H) Relationships between N₂O and CO₂ emissions measured in the last experimental week (week 12), (C, F) soil water contents and (D,G) X-ray tomography estimated volumes of macropores in the topsoil layer (L1) and (D, H) in the bottom layer for N₂O (L4). When significant, the linear regression line and 95% confidence intervals are displayed along with the regression equation, coefficient of determination (R²) and p value. Note that this relationship may differ from the mixed-effect model results (Tables 1 and 4).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g004

(-26.0% in L2, -23.8% in L3, -18.1% in L4). Macropore volumes in L2–L4 were affected by the earthworm treatment, with the highest volume in the mesocosms with endogeic earthworms, followed by both the anecic and control treatments (S2 Fig). Similar trends were observed when macropores were differentiated into burrows and cracks (S1 Table and S2 Fig). Burrow volume was largely driven by the earthworm treatment, with high coefficients of determination in L2, L3, and L4 and in total (i.e., the sum of L1 to L4). Plant presence only affected the total burrow volume and that in L2 (S1 Table).

Exploration of multiple predictors for final N₂O and CO₂ fluxes

Out of the 16 tested potential predictors (S3 Fig), the MCP-penalized multiple regression for N₂O emissions indicates that, in addition to a retained positive coefficient for SWC (4.663e⁻⁰⁴), the macropore volumes in the first and fourth soil layers (Vpores_L1 and Vpores_L4) were also retained with negative coefficients (-1.01e⁻⁰⁵ and -4.18e⁻⁰⁵, respectively) at minimum cross-validation error with lambda = 0.001 (Fig 4A-4D). CO₂ emissions were influenced by macropore volume in the topsoil (Vpores_L1), with a negative coefficient (-0.034) at minimum cross-validation error with lambda = 0.321 (Fig 4E-4H). Notably, these selected predictors also ranked among those with the highest correlation coefficients with CO₂ and N₂O fluxes, as shown by the univariate correlations (S3 Fig). In the final step, we examined how the inclusion of porosity metrics influenced the model performance in explaining the emissions from the last experimental week (Table 3). For N₂O emissions, Vpores_L1 and Vpores_L4 could not be incorporated together due to overfitting and lack of convergence, prompting us to run separate models for each porosity variable.

Minimal adequate models for N₂O emissions that included Vpores_L1 explained more variation than without the porosity metric (r²m = 0.86, vs. r²m = 0.69, respectively) and retained two additional interactions, notably the EwxSWC×Vpores_L1 three-way interaction with a
<figure><img src="image_10.png" /><figcaption>**Fig 4.** N₂O and CO₂ emissions as affected by soil water content and macroporosity. (A, E) Total cumulative emissions as affected by the EwxSWC interaction, where SWC represents the 3-month average SWC. (B, C, D, F, G, H) Relationships between N₂O and CO₂ emissions measured in the last experimental week (week 12), (C, F) soil water contents and (D,G) X-ray tomography estimated volumes of macropores in the topsoil layer (L1) and (D, H) in the bottom layer for N₂O (L4). When significant, the linear regression line and 95% confidence intervals are displayed along with the regression equation, coefficient of determination (R²) and p value. Note that this relationship may differ from the mixed-effect model results (Tables 1 and 4).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0289859.g004

(-26.0% in L2, -23.8% in L3, -18.1% in L4). Macropore volumes in L2–L4 were affected by the earthworm treatment, with the highest volume in the mesocosms with endogeic earthworms, followed by both the anecic and control treatments (S2 Fig). Similar trends were observed when macropores were differentiated into burrows and cracks (S1 Table and S2 Fig). Burrow volume was largely driven by the earthworm treatment, with high coefficients of determination in L2, L3, and L4 and in total (i.e., the sum of L1 to L4). Plant presence only affected the total burrow volume and that in L2 (S1 Table).

Exploration of multiple predictors for final N₂O and CO₂ fluxes

Out of the 16 tested potential predictors (S3 Fig), the MCP-penalized multiple regression for N₂O emissions indicates that, in addition to a retained positive coefficient for SWC (4.663e⁻⁰⁴), the macropore volumes in the first and fourth soil layers (Vpores_L1 and Vpores_L4) were also retained with negative coefficients (-1.01e⁻⁰⁵ and -4.18e⁻⁰⁵, respectively) at minimum cross-validation error with lambda = 0.001 (Fig 4A-4D). CO₂ emissions were influenced by macropore volume in the topsoil (Vpores_L1), with a negative coefficient (-0.034) at minimum cross-validation error with lambda = 0.321 (Fig 4E-4H). Notably, these selected predictors also ranked among those with the highest correlation coefficients with CO₂ and N₂O fluxes, as shown by the univariate correlations (S3 Fig). In the final step, we examined how the inclusion of porosity metrics influenced the model performance in explaining the emissions from the last experimental week (Table 3). For N₂O emissions, Vpores_L1 and Vpores_L4 could not be incorporated together due to overfitting and lack of convergence, prompting us to run separate models for each porosity variable.

Minimal adequate models for N₂O emissions that included Vpores_L1 explained more variation than without the porosity metric (r²m = 0.86, vs. r²m = 0.69, respectively) and retained two additional interactions, notably the EwxSWC×Vpores_L1 three-way interaction with a
<table><caption>Table 2. Covariables summary statistics table.</caption><thead><tr><th>Plant</th><th>Earthworm</th><th>Total shoot biomass<br>dry g</th><th>Litter cover<br>%</th><th>DEA<br>g N g soil<sup>-1</sup> h<sup>-1</sup></th><th>Cmic<br>g C-CO<sub>2</sub> g soil<sup>-1</sup> h<sup>-1</sup></th><th>BR<br>g C-CO<sub>2</sub> g soil<sup>-1</sup> h<sup>-1</sup></th><th>Met_Q<br>Ratio</th><th>NO<sub>3</sub><sup>-</sup><br>mg kg<sup>-1</sup></th><th>NH<sub>4</sub><sup>+</sup><br>mg kg<sup>-1</sup></th></tr></thead><tbody><tr><td>- plant</td><td>without</td><td></td><td>100 ± 5.6 c</td><td>0.04 ± 0.01 ab</td><td>56.9 ± 2.5 a</td><td>0.8 ± 0 a</td><td>0.5 ± 0 a</td><td>31.7 ± 5.9 abc</td><td>0.4 ± 0.1 a</td></tr><tr><td>+ plant</td><td>without</td><td>3.6 ± 0.8 a</td><td>100 ± 6.5 c</td><td>0.05 ± 0.01 ab</td><td>62.3 ± 2.9 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>12.9 ± 6.8 a</td><td>0.3 ± 0.1 a</td></tr><tr><td>- plant</td><td>anecic</td><td></td><td>5 ± 6 a</td><td>0.07 ± 0.01 ab</td><td>63 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>54.6 ± 6.3 cd</td><td>0.4 ± 0.1 a</td></tr><tr><td>+ plant</td><td>anecic</td><td>5.6 ± 0.8 a</td><td>6.4 ± 6 a</td><td>0.05 ± 0.01 ab</td><td>58.9 ± 2.6 a</td><td>0.9 ± 0.1 a</td><td>0.5 ± 0 a</td><td>32.4 ± 6.3 abc</td><td>0.2 ± 0.1 a</td></tr><tr><td>- plant</td><td>endogeic</td><td></td><td>100 ± 6 c</td><td>0.04 ± 0.01 a</td><td>54.4 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.6 ± 0 a</td><td>49.1 ± 6.3 bcd</td><td>0.2 ± 0.1 a</td></tr><tr><td>+ plant</td><td>endogeic</td><td>4.5 ± 0.8 a</td><td>100 ± 6 c</td><td>0.05 ± 0.01 ab</td><td>59 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>23.1 ± 6.3 ab</td><td>0.2 ± 0.1 a</td></tr><tr><td>- plant</td><td>both</td><td></td><td>40 ± 6.5 b</td><td>0.06 ± 0.01 ab</td><td>59.7 ± 2.9 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>65.9 ± 6.8 d</td><td>0.2 ± 0.1 a</td></tr><tr><td>+ plant</td><td>both</td><td>5 ± 0.7 a</td><td>15 ± 5.6 ab</td><td>0.07 ± 0.01 b</td><td>61.3 ± 2.5 a</td><td>0.9 ± 0 a</td><td>0.5 ± 0 a</td><td>39.2 ± 5.9 abcd</td><td>0.4 ± 0.1 a</td></tr><tr><td colspan="10">Source Statistical results (minimal adequate models)</td></tr><tr><td></td><td>Ew</td><td>ns</td><td>74.39***</td><td>2.75<sup>+</sup></td><td>ns</td><td>ns</td><td>ns</td><td>15.95***</td><td>2.22<sup>+</sup></td></tr><tr><td></td><td>Plant</td><td>NA</td><td>0.93</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>46.41***</td><td>ns</td></tr><tr><td></td><td>Ew×Plant</td><td>NA</td><td>3.08*</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td></td><td>m<sup>2</sup>r<sup>2</sup></td><td>0.0</td><td>0.91</td><td>0.14</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.73</td><td>0.11</td></tr></tbody></table>

Summary statistics (mean ± standard error) and minimal adequate models for the covariables considered as predictors for the N₂O and CO₂ emissions from the last week of the experiment in addition to the experimental treatments. Different letters represent significantly different treatments according to Tukey's HSD post hoc test. F-values are shown with significance levels:

***P < 0.001

**P < 0.01

*P < 0.05

†P < 0.1. m<sup>2</sup>r<sup>2</sup> = marginal coefficient of determination, Ew = the earthworm treatment, DEA = denitrifying enzyme activity, BR = basal respiration, Met_Q = metabolic quotient.

https://doi.org/10.1371/journal.pone.0289859.t002

positive fitted coefficient (Table 4). The second model that included Vpores_L4 had an intermediate amount of explained variation (r²m = 0.72) and only detected one additional marginally significant SWC×Vpores_L4 interaction (p value = 0.054, Table 4). Regarding CO₂ emissions, the inclusion of Vpores_L1 largely increased the amount of variance explained (r²m = 0.83 vs. r²m = 0.18, Table 4). The four-way interaction Ew×Plant×SWC×Vpores_L1 was significant with a positive fitted coefficient, indicating higher CO₂ emissions in the presence of earthworms (all treatment combinations containing earthworms) and plants under high SWC levels and high macropore volume in L1.

Discussion

To further advance our understanding of the effects of earthworms on GHG emissions, our study was designed to simultaneously investigate the effects of earthworms, plants, soil moisture fluctuations, and their interactions, with an experimental setup allowing earthworms and plants to affect soil water status and macroporosity. In line with our first hypothesis, we found not only that the presence of earthworms did not increase the CO₂ and N₂O cumulative emissions over 12 weeks but also that the presence of the endogeic species *A. icterica* (alone or with the anecic *L. terrestris*) and the presence of plants reduced N₂O cumulative emissions. Furthermore, earthworms, plants, and their interaction modulated SWC fluctuations and jointly affected weekly N₂O and CO₂ emissions. Moreover, we found that GHG emissions were partly explained by increased macropore volume in the first soil layer, resulting from earthworm burrowing activity.

By imposing soil moisture fluctuations and allowing earthworms and plants to modulate these fluctuations, our results illustrate how soil water availability controls N₂O and CO₂
<table><caption>Table 2. Covariables summary statistics table.</caption><thead><tr><th>Plant</th><th>Earthworm</th><th>Total shoot biomass<br>dry g</th><th>Litter cover<br>%</th><th>DEA<br>g N g soil<sup>-1</sup> h<sup>-1</sup></th><th>Cmic<br>g C-CO<sub>2</sub> g soil<sup>-1</sup> h<sup>-1</sup></th><th>BR<br>g C-CO<sub>2</sub> g soil<sup>-1</sup> h<sup>-1</sup></th><th>Met_Q<br>Ratio</th><th>NO<sub>3</sub><sup>-</sup><br>mg kg<sup>-1</sup></th><th>NH<sub>4</sub><sup>+</sup><br>mg kg<sup>-1</sup></th></tr></thead><tbody><tr><td>- plant</td><td>without</td><td></td><td>100 ± 5.6 c</td><td>0.04 ± 0.01 ab</td><td>56.9 ± 2.5 a</td><td>0.8 ± 0 a</td><td>0.5 ± 0 a</td><td>31.7 ± 5.9 abc</td><td>0.4 ± 0.1 a</td></tr><tr><td>+ plant</td><td>without</td><td>3.6 ± 0.8 a</td><td>100 ± 6.5 c</td><td>0.05 ± 0.01 ab</td><td>62.3 ± 2.9 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>12.9 ± 6.8 a</td><td>0.3 ± 0.1 a</td></tr><tr><td>- plant</td><td>anecic</td><td></td><td>5 ± 6 a</td><td>0.07 ± 0.01 ab</td><td>63 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>54.6 ± 6.3 cd</td><td>0.4 ± 0.1 a</td></tr><tr><td>+ plant</td><td>anecic</td><td>5.6 ± 0.8 a</td><td>6.4 ± 6 a</td><td>0.05 ± 0.01 ab</td><td>58.9 ± 2.6 a</td><td>0.9 ± 0.1 a</td><td>0.5 ± 0 a</td><td>32.4 ± 6.3 abc</td><td>0.2 ± 0.1 a</td></tr><tr><td>- plant</td><td>endogeic</td><td></td><td>100 ± 6 c</td><td>0.04 ± 0.01 a</td><td>54.4 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.6 ± 0 a</td><td>49.1 ± 6.3 bcd</td><td>0.2 ± 0.1 a</td></tr><tr><td>+ plant</td><td>endogeic</td><td>4.5 ± 0.8 a</td><td>100 ± 6 c</td><td>0.05 ± 0.01 ab</td><td>59 ± 2.6 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>23.1 ± 6.3 ab</td><td>0.2 ± 0.1 a</td></tr><tr><td>- plant</td><td>both</td><td></td><td>40 ± 6.5 b</td><td>0.06 ± 0.01 ab</td><td>59.7 ± 2.9 a</td><td>0.8 ± 0.1 a</td><td>0.5 ± 0 a</td><td>65.9 ± 6.8 d</td><td>0.2 ± 0.1 a</td></tr><tr><td>+ plant</td><td>both</td><td>5 ± 0.7 a</td><td>15 ± 5.6 ab</td><td>0.07 ± 0.01 b</td><td>61.3 ± 2.5 a</td><td>0.9 ± 0 a</td><td>0.5 ± 0 a</td><td>39.2 ± 5.9 abcd</td><td>0.4 ± 0.1 a</td></tr><tr><td colspan="10">Source: Statistical results (minimal adequate models)</td></tr><tr><td></td><td>Ew</td><td>ns</td><td>74.39***</td><td>2.75<sup>+</sup></td><td>ns</td><td>ns</td><td>ns</td><td>15.95***</td><td>2.22<sup>+</sup></td></tr><tr><td></td><td>Plant</td><td>NA</td><td>0.93</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>46.41***</td><td>ns</td></tr><tr><td></td><td>Ew×Plant</td><td>NA</td><td>3.08*</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td></td><td>m<sup>2</sup>r<sup>2</sup></td><td>0.0</td><td>0.91</td><td>0.14</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.73</td><td>0.11</td></tr></tbody></table>

Summary statistics (mean ± standard error) and minimal adequate models for the covariables considered as predictors for the N₂O and CO₂ emissions from the last week of the experiment in addition to the experimental treatments. Different letters represent significantly different treatments according to Tukey's HSD post hoc test.
F-values are shown with significance levels:

***P < 0.001

**P < 0.01

*P < 0.05

+P < 0.1. m<sup>2</sup>r<sup>2</sup> = marginal coefficient of determination, Ew = the earthworm treatment, DEA = denitrifying enzyme activity, BR = basal respiration, Met_Q = metabolic quotient.

https://doi.org/10.1371/journal.pone.0289859.t002

positive fitted coefficient (Table 4). The second model that included Vpores_L4 had an intermediate amount of explained variation (r²m = 0.72) and only detected one additional marginally significant SWC×Vpores_L4 interaction (p value = 0.054, Table 4). Regarding CO₂ emissions, the inclusion of Vpores_L1 largely increased the amount of variance explained (r²m = 0.83 vs. r²m = 0.18, Table 4). The four-way interaction Ew×Plant×SWC×Vpores_L1 was significant with a positive fitted coefficient, indicating higher CO₂ emissions in the presence of earthworms (all treatment combinations containing earthworms) and plants under high SWC levels and high macropore volume in L1.

Discussion

To further advance our understanding of the effects of earthworms on GHG emissions, our study was designed to simultaneously investigate the effects of earthworms, plants, soil moisture fluctuations, and their interactions, with an experimental setup allowing earthworms and plants to affect soil water status and macroporosity. In line with our first hypothesis, we found not only that the presence of earthworms did not increase the CO₂ and N₂O cumulative emissions over 12 weeks but also that the presence of the endogeic species *A. icterica* (alone or with the anecic *L. terrestris*) and the presence of plants reduced N₂O cumulative emissions. Furthermore, earthworms, plants, and their interaction modulated SWC fluctuations and jointly affected weekly N₂O and CO₂ emissions. Moreover, we found that GHG emissions were partly explained by increased macropore volume in the first soil layer, resulting from earthworm burrowing activity.

By imposing soil moisture fluctuations and allowing earthworms and plants to modulate these fluctuations, our results illustrate how soil water availability controls N₂O and CO₂
<table><caption>Table 3. Covariables summary statistics table.</caption><thead><tr><th>Plant</th><th>Earthworm</th><th>Vpores_L1 cm³</th><th>Vpores_L2 cm³</th><th>Vpores_L3 cm³</th><th>Vpores_L4 cm³</th><th>Vpores_tot cm³</th></tr></thead><tbody><tr><td>- plant</td><td>without</td><td>127.1 ± 10.9 a</td><td>79.7 ± 5.9 bc</td><td>39.4 ± 5.7 ab</td><td>9 ± 5.6 a</td><td>255.2 ± 23 abc</td></tr><tr><td>+ plant</td><td>without</td><td>98.8 ± 12.6 a</td><td>57.4 ± 6.8 ab</td><td>20.3 ± 6.6 a</td><td>2.1 ± 6.5 a</td><td>178.6 ± 26.6 a</td></tr><tr><td>- plant</td><td>anecic</td><td>96.9 ± 11.7 a</td><td>60.3 ± 6.3 ab</td><td>57.6 ± 6.1 bc</td><td>39.3 ± 6 bc</td><td>254.1 ± 24.6 abc</td></tr><tr><td>+ plant</td><td>anecic</td><td>99.3 ± 11.7 a</td><td>45.5 ± 6.3 a</td><td>37.7 ± 6.1 ab</td><td>23.5 ± 6 ab</td><td>206 ± 24.6 ab</td></tr><tr><td>- plant</td><td>endogeic</td><td>118.4 ± 11.7 a</td><td>96.6 ± 6.3 c</td><td>92 ± 6.1 d</td><td>77.3 ± 6 d</td><td>384.2 ± 24.6 d</td></tr><tr><td>+ plant</td><td>endogeic</td><td>92.6 ± 11.7 a</td><td>71.8 ± 6.3 abc</td><td>72.2 ± 6.1 cd</td><td>63.5 ± 6 cd</td><td>300 ± 24.6 bcd</td></tr><tr><td>- plant</td><td>both</td><td>105.6 ± 12.6 a</td><td>77.3 ± 6.8 bc</td><td>78.8 ± 6.6 cd</td><td>65.7 ± 6.5 cd</td><td>327.4 ± 26.6 cd</td></tr><tr><td>+ plant</td><td>both</td><td>112.5 ± 10.9 a</td><td>55.1 ± 5.9 ab</td><td>63.2 ± 5.7 bc</td><td>53.6 ± 5.6 cd</td><td>284.4 ± 23 abcd</td></tr><tr><td colspan="7">Source: Statistical results (minimal adequate models)</td></tr><tr><td></td><td>Ew</td><td>ns</td><td>9.04***</td><td>32.04***</td><td>191.71***</td><td>11.91***</td></tr><tr><td></td><td>Plant</td><td>ns</td><td>26.63***</td><td>21.36***</td><td>13.83***</td><td>15.78***</td></tr><tr><td></td><td>Ew×Plant</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td></td><td>m r²</td><td>0.00</td><td>0.48</td><td>0.68</td><td>0.94</td><td>0.59</td></tr></tbody></table>

Summary statistics (mean ± standard error) and minimal adequate models of the porosity variables considered as predictors for the N₂O and CO₂ emissions from the last week of the experiment in addition to the experimental treatments. Different letters represent significantly different treatments according to Tukey's HSD post hoc test. F-values are shown with significance levels:

***P < 0.001

**P < 0.01

*P < 0.05

+P < 0.1. m r² = marginal coefficient of determination, Ew = the earthworm treatment, Vpores = total macroporosity in the four different soil layers from L1 (0–8.5 cm) to L4 (25.5–34 cm).

https://doi.org/10.1371/journal.pone.0289859.t003

emissions in complex ways [17, 64]. In general, a combination of limited substrate diffusion at very low water content and limited gas diffusion at high water content leads to maximal N₂O emissions (via nitrification and denitrification) and CO₂ emissions (via respiration) at intermediate SWC, approximately 75% of the water-filled pore space [22, 23, 65]. In our experiment, SWC varied considerably, with the lowest values in the presence of plants and anecic earthworms, alone or mixed with endogeic earthworms (Fig 1A and 1B). The significant Ew×SWC interaction (Table 1) observed for cumulative and weekly N₂O and CO₂ emissions illustrates the SWC optimal value phenomenon, which occurred at approximately 85% of field capacity or 59% of water-filled pore space (Fig 4A and 4B). Indeed, anecic earthworms create large vertical burrows, increasing water infiltration, and bury leaf litter, increasing the proportions of bare soil and water evaporation [40]. This led to a SWC value lower than the optimum SWC value for microbial activity, thus explaining the positive relationship between SWC and respiration for this species. Conversely, the treatment combinations with the endogeic species, similar to the control with no earthworms, maintained higher than optimal SWC values for soil respiration on average, thus explaining the negative slopes of CO₂ fluxes with increasing SWC. The presence of both earthworm species led to SWC values that spanned across the optimum, and no clear relationship between SWC and CO₂ could be detected. Simultaneously, the presence of *B. distachyon* grass lowered the average SWC compared to the mesocosms without plants (mean ± se = 82.0 ± 0.6% of WHC), and in this soil moisture range, CO₂ emissions increased with SWC (S5 Fig). In the absence of *B. distachyon*, the average SWC was higher (86.2 ± 0.6% of WHC), but increasing the SWC lowered CO₂ emissions as the range of SWC was beyond the optimum, and presumably, soil respiration was limited by O₂ diffusivity under these conditions (S5 Fig). Regarding N₂O emissions, we observed similar patterns in most cases, except for the control group without earthworms, where emissions still exceeded the
<table><caption>Table 3. Covariables summary statistics table.</caption><thead><tr><th>Plant</th><th>Earthworm</th><th>Vpores_L1 cm³</th><th>Vpores_L2 cm³</th><th>Vpores_L3 cm³</th><th>Vpores_L4 cm³</th><th>Vpores_tot cm³</th></tr></thead><tbody><tr><td>- plant</td><td>without</td><td>127.1 ± 10.9 a</td><td>79.7 ± 5.9 bc</td><td>39.4 ± 5.7 ab</td><td>9 ± 5.6 a</td><td>255.2 ± 23 abc</td></tr><tr><td>+ plant</td><td>without</td><td>98.8 ± 12.6 a</td><td>57.4 ± 6.8 ab</td><td>20.3 ± 6.6 a</td><td>2.1 ± 6.5 a</td><td>178.6 ± 26.6 a</td></tr><tr><td>- plant</td><td>anecic</td><td>96.9 ± 11.7 a</td><td>60.3 ± 6.3 ab</td><td>57.6 ± 6.1 bc</td><td>39.3 ± 6 bc</td><td>254.1 ± 24.6 abc</td></tr><tr><td>+ plant</td><td>anecic</td><td>99.3 ± 11.7 a</td><td>45.5 ± 6.3 a</td><td>37.7 ± 6.1 ab</td><td>23.5 ± 6 ab</td><td>206 ± 24.6 ab</td></tr><tr><td>- plant</td><td>endogeic</td><td>118.4 ± 11.7 a</td><td>96.6 ± 6.3 c</td><td>92 ± 6.1 d</td><td>77.3 ± 6 d</td><td>384.2 ± 24.6 d</td></tr><tr><td>+ plant</td><td>endogeic</td><td>92.6 ± 11.7 a</td><td>71.8 ± 6.3 abc</td><td>72.2 ± 6.1 cd</td><td>63.5 ± 6 cd</td><td>300 ± 24.6 bcd</td></tr><tr><td>- plant</td><td>both</td><td>105.6 ± 12.6 a</td><td>77.3 ± 6.8 bc</td><td>78.8 ± 6.6 cd</td><td>65.7 ± 6.5 cd</td><td>327.4 ± 26.6 cd</td></tr><tr><td>+ plant</td><td>both</td><td>112.5 ± 10.9 a</td><td>55.1 ± 5.9 ab</td><td>63.2 ± 5.7 bc</td><td>53.6 ± 5.6 cd</td><td>284.4 ± 23 abcd</td></tr><tr><td colspan="7">Source: Statistical results (minimal adequate models)</td></tr><tr><td></td><td>Ew</td><td>ns</td><td>9.04***</td><td>32.04***</td><td>191.71***</td><td>11.91***</td></tr><tr><td></td><td>Plant</td><td>ns</td><td>26.63***</td><td>21.36***</td><td>13.83***</td><td>15.78***</td></tr><tr><td></td><td>Ew×Plant</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td><td>ns</td></tr><tr><td></td><td>m r²</td><td>0.00</td><td>0.48</td><td>0.68</td><td>0.94</td><td>0.59</td></tr></tbody></table>

Summary statistics (mean ± standard error) and minimal adequate models of the porosity variables considered as predictors for the N₂O and CO₂ emissions from the last week of the experiment in addition to the experimental treatments. Different letters represent significantly different treatments according to Tukey's HSD post hoc test. F-values are shown with significance levels:

***P < 0.001

**P < 0.01

*P < 0.05

+P < 0.1. m r² = marginal coefficient of determination, Ew = the earthworm treatment, Vpores = total macroporosity in the four different soil layers from L1 (0–8.5 cm) to L4 (25.5–34 cm).

https://doi.org/10.1371/journal.pone.0289859.t003

emissions in complex ways [17, 64]. In general, a combination of limited substrate diffusion at very low water content and limited gas diffusion at high water content leads to maximal N₂O emissions (via nitrification and denitrification) and CO₂ emissions (via respiration) at intermediate SWC, approximately 75% of the water-filled pore space [22, 23, 65]. In our experiment, SWC varied considerably, with the lowest values in the presence of plants and anecic earthworms, alone or mixed with endogeic earthworms (Fig 1A and 1B). The significant Ew×SWC interaction (Table 1) observed for cumulative and weekly N₂O and CO₂ emissions illustrates the SWC optimal value phenomenon, which occurred at approximately 85% of field capacity or 59% of water-filled pore space (Fig 4A and 4B). Indeed, anecic earthworms create large vertical burrows, increasing water infiltration, and bury leaf litter, increasing the proportions of bare soil and water evaporation [40]. This led to a SWC value lower than the optimum SWC value for microbial activity, thus explaining the positive relationship between SWC and respiration for this species. Conversely, the treatment combinations with the endogeic species, similar to the control with no earthworms, maintained higher than optimal SWC values for soil respiration on average, thus explaining the negative slopes of CO₂ fluxes with increasing SWC. The presence of both earthworm species led to SWC values that spanned across the optimum, and no clear relationship between SWC and CO₂ could be detected. Simultaneously, the presence of *B. distachyon* grass lowered the average SWC compared to the mesocosms without plants (mean ± se = 82.0 ± 0.6% of WHC), and in this soil moisture range, CO₂ emissions increased with SWC (S5 Fig). In the absence of *B. distachyon*, the average SWC was higher (86.2 ± 0.6% of WHC), but increasing the SWC lowered CO₂ emissions as the range of SWC was beyond the optimum, and presumably, soil respiration was limited by O₂ diffusivity under these conditions (S5 Fig). Regarding N₂O emissions, we observed similar patterns in most cases, except for the control group without earthworms, where emissions still exceeded the
<table><caption>Table 4. Summary statistic table of last week gas fluxes.</caption><thead><tr><th rowspan="2">Source</th><th colspan="3">N<sub>2</sub>O model</th><th colspan="2">CO<sub>2</sub> model</th></tr><tr><th>without porosity</th><th>With pore volume L1</th><th>with pore volume L4</th><th>without pore volume</th><th>with pore volume L1</th></tr></thead><tbody><tr><td>Ew</td><td>15.46***</td><td>24.51***</td><td>16.15***</td><td>ns</td><td>2.1</td></tr><tr><td>Plant</td><td>ns</td><td>1.77</td><td>ns</td><td>6.95*</td><td>15.9***</td></tr><tr><td>SWC</td><td>11.93**</td><td>28.25***</td><td>13.15***</td><td>4.12*</td><td>14.33**</td></tr><tr><td>Vpore</td><td>NA</td><td>0.63</td><td>ns</td><td>NA</td><td>36.98***</td></tr><tr><td>Ew×Plant</td><td>ns</td><td>1.22</td><td>ns</td><td>ns</td><td>1.27</td></tr><tr><td>Ew×SWC</td><td>12.19***</td><td>18.73***</td><td>12.81***</td><td>ns</td><td>2.86<sup>+</sup></td></tr><tr><td>Plant×SWC</td><td>ns</td><td>2.00</td><td>ns</td><td>ns</td><td>0.32</td></tr><tr><td>Ew×Vpore</td><td>NA</td><td>5.60**</td><td>ns</td><td>NA</td><td>0.14</td></tr><tr><td>Plant×Vpore</td><td>NA</td><td>1.04</td><td>ns</td><td>NA</td><td>0.15</td></tr><tr><td>SWC×Vpore</td><td>NA</td><td>0.66</td><td>3.83<sup>+</sup></td><td>NA</td><td>0.49</td></tr><tr><td>Ew×Plant×SWC</td><td>ns</td><td>1.19</td><td>ns</td><td>ns</td><td>1.81</td></tr><tr><td>Ew×Plant×Vpore</td><td>NA</td><td>2.04</td><td>ns</td><td>NA</td><td>1.14</td></tr><tr><td>Ew×SWC×Vpore</td><td>NA</td><td>11.79***</td><td>ns</td><td>NA</td><td>1.39</td></tr><tr><td>Plant×SWC×Vpore</td><td>NA</td><td>1.95</td><td>ns</td><td>NA</td><td>1.94</td></tr><tr><td>Ew×Plant×SWC×Vpore</td><td>NA</td><td>ns</td><td>ns</td><td>NA</td><td>3.76*</td></tr><tr><td>m<sub>r</sub>r<sup>2</sup></td><td>0.69</td><td>0.86</td><td>0.72</td><td>0.18</td><td>0.83</td></tr></tbody></table>

Minimal adequate models presenting the results explaining the CO₂ and N₂O fluxes from the last sampling (week 12) where the soil porosity-related variables were included in the model as potential predictors (compared with the models without the soil-porosity variables). “NA” stands for non-applicable, “ns” stands for variables that were not significant and were not retained in the minimal adequate models and m<sub>r</sub>r² represents the marginal coefficient of determination. F-values are shown with significance levels:

https://doi.org/10.1371/journal.pone.0289859.t004

optimal SWC values observed in the treatments with earthworms. This finding only partially supports the hypothesis of an optimal SWC mechanism and its effect on N₂O emissions. These results suggest that while the interactions between earthworms, plants, and SWC strongly influence N₂O emissions, other factors likely come into play when earthworms are absent, such as a significantly different soil porosity status.

The inclusion of the soil porosity data revealed that the total volume occupied by soil macro-pores in the upper soil layer was an important predictor of GHG emissions (with a negative coefficient, Fig 4G). This suggests that increasing porosity/aeration in the topsoil layer (0–8.5 cm) can decrease N₂O and CO₂ emissions, presumably by reducing the SWC in the upper and most microbially active soil layers. Interestingly, in line with our third hypothesis, porosity in the bottom layer (25.5–34 cm) was a good predictor of N₂O emissions (Fig 4D) and was the only variable that was influenced by earthworm species in the same way as cumulative N₂O emissions. Indeed, the number of burrows in the deepest layer was higher in the presence of the endogeic *A. icterica* (alone or alongside the anecic species), a species with high affinity for the deepest soil layers [66], and presumably prevented the development of denitrification-stimulating anaerobic sites. The reduction of N₂O emissions via increased soil aeration was also previously suggested by several studies [67, 68] but was not explicitly shown to our knowledge. Our results indicate that this effect was more prevalent in the presence of the endogeic species and seemed to be related to the higher number of burrows produced by this species in contrast to the larger but less numerous semipermanent burrows produced by the anecic species (S3 Fig)

<footnote>***P < 0.001</footnote>

<footnote>**P < 0.01</footnote>

<footnote>*P < 0.05</footnote>

<footnote>+P < 0.1.</footnote>
<table><caption>Table 4. Summary statistic table of last week gas fluxes.</caption><thead><tr><th rowspan="2">Source</th><th colspan="3">N<sub>2</sub>O model</th><th colspan="2">CO<sub>2</sub> model</th></tr><tr><th>without porosity</th><th>With pore volume L1</th><th>with pore volume L4</th><th>without pore volume</th><th>with pore volume L1</th></tr></thead><tbody><tr><td>Ew</td><td>15.46***</td><td>24.51***</td><td>16.15***</td><td>ns</td><td>2.1</td></tr><tr><td>Plant</td><td>ns</td><td>1.77</td><td>ns</td><td>6.95*</td><td>15.9***</td></tr><tr><td>SWC</td><td>11.93**</td><td>28.25***</td><td>13.15***</td><td>4.12*</td><td>14.33**</td></tr><tr><td>Vpore</td><td>NA</td><td>0.63</td><td>ns</td><td>NA</td><td>36.98***</td></tr><tr><td>Ew×Plant</td><td>ns</td><td>1.22</td><td>ns</td><td>ns</td><td>1.27</td></tr><tr><td>Ew×SWC</td><td>12.19***</td><td>18.73***</td><td>12.81***</td><td>ns</td><td>2.86<sup>+</sup></td></tr><tr><td>Plant×SWC</td><td>ns</td><td>2.00</td><td>ns</td><td>ns</td><td>0.32</td></tr><tr><td>Ew×Vpore</td><td>NA</td><td>5.60**</td><td>ns</td><td>NA</td><td>0.14</td></tr><tr><td>Plant×Vpore</td><td>NA</td><td>1.04</td><td>ns</td><td>NA</td><td>0.15</td></tr><tr><td>SWC×Vpore</td><td>NA</td><td>0.66</td><td>3.83<sup>+</sup></td><td>NA</td><td>0.49</td></tr><tr><td>Ew×Plant×SWC</td><td>ns</td><td>1.19</td><td>ns</td><td>ns</td><td>1.81</td></tr><tr><td>Ew×Plant×Vpore</td><td>NA</td><td>2.04</td><td>ns</td><td>NA</td><td>1.14</td></tr><tr><td>Ew×SWC×Vpore</td><td>NA</td><td>11.79***</td><td>ns</td><td>NA</td><td>1.39</td></tr><tr><td>Plant×SWC×Vpore</td><td>NA</td><td>1.95</td><td>ns</td><td>NA</td><td>1.94</td></tr><tr><td>Ew×Plant×SWC×Vpore</td><td>NA</td><td>ns</td><td>ns</td><td>NA</td><td>3.76*</td></tr><tr><td>m<sub>r</sub>r<sup>2</sup></td><td>0.69</td><td>0.86</td><td>0.72</td><td>0.18</td><td>0.83</td></tr></tbody></table>

Minimal adequate models presenting the results explaining the CO₂ and N₂O fluxes from the last sampling (week 12) where the soil porosity-related variables were included in the model as potential predictors (compared with the models without the soil-porosity variables). “NA” stands for non-applicable, “ns” stands for variables that were not significant and were not retained in the minimal adequate models and m<sub>r</sub>r² represents the marginal coefficient of determination. F-values are shown with significance levels:

https://doi.org/10.1371/journal.pone.0289859.t004

optimal SWC values observed in the treatments with earthworms. This finding only partially supports the hypothesis of an optimal SWC mechanism and its effect on N₂O emissions. These results suggest that while the interactions between earthworms, plants, and SWC strongly influence N₂O emissions, other factors likely come into play when earthworms are absent, such as a significantly different soil porosity status.

The inclusion of the soil porosity data revealed that the total volume occupied by soil macro-pores in the upper soil layer was an important predictor of GHG emissions (with a negative coefficient, Fig 4G). This suggests that increasing porosity/aeration in the topsoil layer (0–8.5 cm) can decrease N₂O and CO₂ emissions, presumably by reducing the SWC in the upper and most microbially active soil layers. Interestingly, in line with our third hypothesis, porosity in the bottom layer (25.5–34 cm) was a good predictor of N₂O emissions (Fig 4D) and was the only variable that was influenced by earthworm species in the same way as cumulative N₂O emissions. Indeed, the number of burrows in the deepest layer was higher in the presence of the endogeic *A. icterica* (alone or alongside the anecic species), a species with high affinity for the deepest soil layers [66], and presumably prevented the development of denitrification-stimulating anaerobic sites. The reduction of N₂O emissions via increased soil aeration was also previously suggested by several studies [67, 68] but was not explicitly shown to our knowledge. Our results indicate that this effect was more prevalent in the presence of the endogeic species and seemed to be related to the higher number of burrows produced by this species in contrast to the larger but less numerous semipermanent burrows produced by the anecic species (S3 Fig)

<footnote>***P < 0.001</footnote>

<footnote>**P < 0.01</footnote>

<footnote>*P < 0.05</footnote>

<footnote>+P < 0.1.</footnote>
[41]. In contrast, plants significantly reduced soil macropore volume, likely due to roots improving soil structure and stability against crack formation through the production of exudates acting as binding agents or by root mechanical engagement with soil aggregates [69, 70].

Our experiment also allows us to discuss the importance of nutrient availability for GHG emissions. The observed 19.8% reduction in N₂O emissions in the presence of plants also occurred, likely in part due to plant N uptake, as we found that the amounts of soil NO₃⁻ and NH₄⁺ were 43% and 20% lower, respectively, in the presence of plants (in line with our second hypothesis), independent of earthworm presence. This finding supports our hypothesis that plants can compete with microorganisms for nutrients and therefore limit bulk microbial activity [37], given the importance of nitrogen availability for nitrification and denitrification [71]. The absence of a positive effect of plants on CO₂ fluxes (either weekly fluxes or microbial potential activity at final harvest) is surprising, notably because our experimental design only allowed the combined measurements of CO₂ originating from heterotrophic and root respiration. This could be explained by either nutrient (nitrate and ammonium) limitations or an overall low plant effect due to the relatively low plant biomass production of *B. distachyon* in our experiment. The soil NO₃⁻ concentration increased in the presence of anecic and endogeic earthworms and even more so when both ecological categories were present, leading to 2.1- and 3-fold increases relative to the control in mesocosms with and without plants, respectively. These results can be explained by the combined effect of local vertical litter burial by *L. terres-tris* and the horizontal redistribution in the extensive burrow system of *A. icterica* as well as the higher nitrogen concentration in earthworm casts compared to bulk soil [26, 72]. The accelerated burial of surface litter by the anecic species therefore likely contributed to the higher N₂O emissions via increased N availability after watering events, but this effect faded with time and soil drying. Simultaneously, this higher nutrient availability likely contributed to higher plant growth in the presence of earthworms [6]. Despite this increase in nutrient availability, the decreased soil moisture due to earthworm burrowing and the formation of cracks in the top-soil still reduced microbial activity and GHG production. Overall, our study supports the idea that investigating the effect of earthworms on GHG emissions requires the use of an experimental setup that includes plants over a sufficiently long period (> 3 months) [7] and allows water availability to fluctuate due to biological activity.

Mesocosm experiments are highly valuable tools for global change research, but care must be taken in interpreting and extrapolating results, and potential caveats of our study must be mentioned [73, 74]. Because soil properties can strongly influence GHG emissions as well as earthworm cast properties [75, 76], we cannot be sure of the transferability of our results to other soil types. Second, as gas diffusion will occur at the soil-air interface, with gases that are more concentrated in the soil moving toward the atmospheric air, where the concentrations are lower, the bottoms of mesocosms should be as airtight as possible while still allowing for drainage. In our case, the holes at the bottom represented 0.4% of the surface area and were obstructed by the table, thus presumably limiting this bias. Future studies could address this issue by placing the mesocosms on a layer of the same experimental soil or using active drainage systems consisting of suction pumps and tubing equipped with valves that can be closed after drainage. We also acknowledge that the size of our mesocosms, although larger than those in many other studies, could have interfered with earthworm burrowing behavior, especially for deep-burrowing anecic earthworms [51]. Furthermore, as only one earthworm species per ecological category was used, it is unknown whether our findings are transferable to other species from the same ecological categories or whether these findings are also valid for epigeic earthworm species that have been reported to also increase N₂O emissions [7]. It must be noted that earthworm ecological categories were not conceptualized to describe a functional role but rather ecological and morphological groups, which can also explain the high
[41]. In contrast, plants significantly reduced soil macropore volume, likely due to roots improving soil structure and stability against crack formation through the production of exudates acting as binding agents or by root mechanical engagement with soil aggregates [69, 70].

Our experiment also allows us to discuss the importance of nutrient availability for GHG emissions. The observed 19.8% reduction in N₂O emissions in the presence of plants also occurred, likely in part due to plant N uptake, as we found that the amounts of soil NO₃⁻ and NH₄⁺ were 43% and 20% lower, respectively, in the presence of plants (in line with our second hypothesis), independent of earthworm presence. This finding supports our hypothesis that plants can compete with microorganisms for nutrients and therefore limit bulk microbial activity [37], given the importance of nitrogen availability for nitrification and denitrification [71]. The absence of a positive effect of plants on CO₂ fluxes (either weekly fluxes or microbial potential activity at final harvest) is surprising, notably because our experimental design only allowed the combined measurements of CO₂ originating from heterotrophic and root respiration. This could be explained by either nutrient (nitrate and ammonium) limitations or an overall low plant effect due to the relatively low plant biomass production of *B. distachyon* in our experiment. The soil NO₃⁻ concentration increased in the presence of anecic and endogeic earthworms and even more so when both ecological categories were present, leading to 2.1- and 3-fold increases relative to the control in mesocosms with and without plants, respectively. These results can be explained by the combined effect of local vertical litter burial by *L. terres-tris* and the horizontal redistribution in the extensive burrow system of *A. icterica* as well as the higher nitrogen concentration in earthworm casts compared to bulk soil [26, 72]. The accelerated burial of surface litter by the anecic species therefore likely contributed to the higher N₂O emissions via increased N availability after watering events, but this effect faded with time and soil drying. Simultaneously, this higher nutrient availability likely contributed to higher plant growth in the presence of earthworms [6]. Despite this increase in nutrient availability, the decreased soil moisture due to earthworm burrowing and the formation of cracks in the top-soil still reduced microbial activity and GHG production. Overall, our study supports the idea that investigating the effect of earthworms on GHG emissions requires the use of an experimental setup that includes plants over a sufficiently long period (> 3 months) [7] and allows water availability to fluctuate due to biological activity.

Mesocosm experiments are highly valuable tools for global change research, but care must be taken in interpreting and extrapolating results, and potential caveats of our study must be mentioned [73, 74]. Because soil properties can strongly influence GHG emissions as well as earthworm cast properties [75, 76], we cannot be sure of the transferability of our results to other soil types. Second, as gas diffusion will occur at the soil-air interface, with gases that are more concentrated in the soil moving toward the atmospheric air, where the concentrations are lower, the bottoms of mesocosms should be as airtight as possible while still allowing for drainage. In our case, the holes at the bottom represented 0.4% of the surface area and were obstructed by the table, thus presumably limiting this bias. Future studies could address this issue by placing the mesocosms on a layer of the same experimental soil or using active drainage systems consisting of suction pumps and tubing equipped with valves that can be closed after drainage. We also acknowledge that the size of our mesocosms, although larger than those in many other studies, could have interfered with earthworm burrowing behavior, especially for deep-burrowing anecic earthworms [51]. Furthermore, as only one earthworm species per ecological category was used, it is unknown whether our findings are transferable to other species from the same ecological categories or whether these findings are also valid for epigeic earthworm species that have been reported to also increase N₂O emissions [7]. It must be noted that earthworm ecological categories were not conceptualized to describe a functional role but rather ecological and morphological groups, which can also explain the high
variability of earthworm species effects within the same ecological category [77]. To the best of our knowledge, our study is the first to investigate the link between earthworm-induced macroporosity and greenhouse gas fluxes; however, the size of the mesocosms combined with the imposed drying-rewetting cycles led to the formation of cracks that unfortunately made the analysis more difficult. As cracks are even more unstable than burrows under drying-rewetting cycles [41], they should be avoided if possible or taken into account in analyses in future experiments. Another caveat is that our weekly measurement frequency of CO₂ and N₂O fluxes over 12 weeks may have missed daily variations or higher emission peaks following the watering events. However, our measurements still detected peaks after watering (Fig 2), and while stimulation of emissions under the anecic treatment was detectable, for the endogeic species, the N₂O emission rates were consistently lower than those of the control. Finally, our study duration (3 months) was intermediate, as classified by Lubbers et al. 2013, and earthworm densities in two of our three earthworm treatments were higher than natural densities. Future studies should perform longer experiments with natural densities [7], which, as in our case, can prove to be technically challenging.

In conclusion, our study highlights new mechanisms by which earthworms and plants influence soil GHG emissions in an experimental setup integrating earthworm engineering effects on soil water fluxes and soil porosity, two major mechanisms that have been neglected thus far. The presence of earthworms did not increase CO₂ and N₂O emissions and revealed that the endogeic earthworm *A. icterica*, a common species present in Europe and North America, even has the potential to reduce N₂O emissions. Our study is an additional step toward a better understanding of the interactions between soil biota and soil physicochemical properties underlying GHG emissions. Future research on these mechanisms would be highly valuable, especially in an agricultural context, as agriculture is the first sector of N₂O emissions [78], and many mitigation practices have been proposed (e.g., reduced tillage or cover crops) [69]. At the same time, these practices also affect earthworm communities, soil porosity, soil compaction, and water infiltration, which interact and affect soil functions [79–82]. Future research should address these points in experimental setups where the earthworm engineering effect on soil water status and aeration is allowed to take place in a realistic manner.

## Supporting information

**S1 Fig. Schematic depicting the mesocosms dimensions and the elements (base collar and static chamber) used for measuring the N₂O and CO₂ emissions.** L1 (0–8.5 cm), L2 (8.5–17 cm), L3 (17–25.5 cm) and L4 (25.5–34 cm depth) represent the four different soil layers that were separately analyzed for soil porosity variables.
(TIF)

**S2 Fig.**
**A)** The effects of treatments on macroporosity volume differentiated as earthworm burrows, cracks, and total macroporosity (burrows + cracks). Error bars represent ± 1 SEM (L1 for 0–8.5 cm, L2 for 8.5–17 cm, L3 for 17–25.5 cm and L4 for 25.5–34 cm depth; see Fig 1).
**B)** Examples of 3D reconstruction of the soil macroporosity differentiated as burrows, cracks and total (burrows + cracks) for the two earthworm species alone.
(TIF)

**S3 Fig. Correlation matrix of all predictors measured in the last week of the experiment.**
Ew_bm = Earthworm biomass (g FW mesocosm-1), Ew_no = Number of added earthworms per mesocosm (number), SWC = Soil water content relative to field capacity (% of field capacity), Plant_bm = Aboveground plant biomass (g DW mesocosm⁻¹), Litter = Percentage of soil surface covered by litter (%), EA = Potential denitrification enzymatic activity (µg N g⁻¹ soil
variability of earthworm species effects within the same ecological category [77]. To the best of our knowledge, our study is the first to investigate the link between earthworm-induced macroporosity and greenhouse gas fluxes; however, the size of the mesocosms combined with the imposed drying-rewetting cycles led to the formation of cracks that unfortunately made the analysis more difficult. As cracks are even more unstable than burrows under drying-rewetting cycles [41], they should be avoided if possible or taken into account in analyses in future experiments. Another caveat is that our weekly measurement frequency of CO₂ and N₂O fluxes over 12 weeks may have missed daily variations or higher emission peaks following the watering events. However, our measurements still detected peaks after watering (Fig 2), and while stimulation of emissions under the anecic treatment was detectable, for the endogeic species, the N₂O emission rates were consistently lower than those of the control. Finally, our study duration (3 months) was intermediate, as classified by Lubbers et al. 2013, and earthworm densities in two of our three earthworm treatments were higher than natural densities. Future studies should perform longer experiments with natural densities [7], which, as in our case, can prove to be technically challenging.

In conclusion, our study highlights new mechanisms by which earthworms and plants influence soil GHG emissions in an experimental setup integrating earthworm engineering effects on soil water fluxes and soil porosity, two major mechanisms that have been neglected thus far. The presence of earthworms did not increase CO₂ and N₂O emissions and revealed that the endogeic earthworm *A. icterica*, a common species present in Europe and North America, even has the potential to reduce N₂O emissions. Our study is an additional step toward a better understanding of the interactions between soil biota and soil physicochemical properties underlying GHG emissions. Future research on these mechanisms would be highly valuable, especially in an agricultural context, as agriculture is the first sector of N₂O emissions [78], and many mitigation practices have been proposed (e.g., reduced tillage or cover crops) [69]. At the same time, these practices also affect earthworm communities, soil porosity, soil compaction, and water infiltration, which interact and affect soil functions [79–82]. Future research should address these points in experimental setups where the earthworm engineering effect on soil water status and aeration is allowed to take place in a realistic manner.

**Supporting information**

**S1 Fig. Schematic depicting the mesocosms dimensions and the elements (base collar and static chamber) used for measuring the N₂O and CO₂ emissions.** L1 (0–8.5 cm), L2 (8.5–17 cm), L3 (17–25.5 cm) and L4 (25.5–34 cm depth) represent the four different soil layers that were separately analyzed for soil porosity variables.
(TIF)

**S2 Fig.**
**A)** The effects of treatments on macroporosity volume differentiated as earthworm burrows, cracks, and total macroporosity (burrows + cracks). Error bars represent ± 1 SEM (L1 for 0–8.5 cm, L2 for 8.5–17 cm, L3 for 17–25.5 cm and L4 for 25.5–34 cm depth; see Fig 1).
**B)** Examples of 3D reconstruction of the soil macroporosity differentiated as burrows, cracks and total (burrows + cracks) for the two earthworm species alone.
(TIF)

**S3 Fig. Correlation matrix of all predictors measured in the last week of the experiment.**
Ew_bm = Earthworm biomass (g FW mesocosm-1), Ew_no = Number of added earthworms per mesocosm (number), SWC = Soil water content relative to field capacity (% of field capacity), Plant_bm = Aboveground plant biomass (g DW mesocosm⁻¹), Litter = Percentage of soil surface covered by litter (%), EA = Potential denitrification enzymatic activity (µg N g⁻¹ soil
DW h-1), Cmic = Microbial biomass C (µg Cmic g soil⁻¹ DW), BR = Microbial basal respiration (µg C-CO₂ g⁻¹ soil DW h⁻¹), Met_Q = Microbial metabolic quotient (µg C-CO₂ µg⁻¹ h⁻¹), NH₄⁺ = Ammonium content in soil at the end of the experiment (mg kg⁻¹), NO₃⁻ = Nitrate content in soil at the end of the experiment (mg kg⁻¹), Vpores_L1-4 & tot = Macropore (burrow + cracks) volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³), Vburrows_L1-4 & tot = Burrow volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³), Vcracks_L1-4 & tot = Cracks volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³).

(TIF)

**S4 Fig.** Effects of experimental treatments on several predictors of CO₂ and N₂O fluxes measured at the end of the experiment. Green = in presence of plant, dark grey = in absence of plant. Different letters represent significantly different treatments according to Tukey's HSD post hoc test. Error bars represent ± 1 SEM.

(TIF)

**S5 Fig.** Cumulative N₂O (A) and CO₂ (B) emissions as affected by Plant×SWC interaction, where SWC represents the 3-month average SWC and the emissions the total cumulative gas emissions.

(TIF)

**S1 Table.** Effects earthworms (Ew) and plant treatments on the total macroporosity volume (pores) as well as differentiated as burrows and cracks (see Table 2 for detailed variable description). The “ns” abbreviation stands for variables that were not significant and were not retained in the minimal adequate models whereas $m_r^2$ represents the marginal coefficient of determination. ***P < 0.001; **P < 0.01; *P < 0.05; †P < 0.1.

(DOCX)

## Acknowledgments

This study benefited from the CNRS human and technical resources allocated to the ECO-TRONS Research Infrastructure and AnaEE France. We thank Thierry Morvan for his assistance with the organization of the soil excavation at the EFELE experimental site as well as Thierry Mathieu and David Degueldre from the technical platform Terrain d'Expériences du C.E.F.E. for their support with the greenhouse and production of gas sampling collars. Microbial and soil analyses were done at Plateforme d'Analyses Chimiques en Ecologie (PACE) supported by the LabEx CeMEB, an ANR "Investissements d'avenir" programme (ANR-10-LABX-04-01).

## Author Contributions

**Conceptualization:** Johanne Nahmani, Alexandru Milcu.

**Data curation:** Pierre Ganault, Yvan Capowiez, Alexandru Milcu.

**Formal analysis:** Pierre Ganault, Alexandru Milcu.

**Funding acquisition:** Johanne Nahmani, Alexandru Milcu.

**Investigation:** Pierre Ganault.
DW h-1), Cmic = Microbial biomass C (µg Cmic g soil⁻¹ DW), BR = Microbial basal respiration (µg C-CO₂ g⁻¹ soil DW h⁻¹), Met_Q = Microbial metabolic quotient (µg C-CO₂ µg⁻¹ h⁻¹), NH₄⁺ = Ammonium content in soil at the end of the experiment (mg kg⁻¹), NO₃⁻ = Nitrate content in soil at the end of the experiment (mg kg⁻¹), Vpores_L1-4 & tot = Macropore (burrow + cracks) volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³), Vburrows_L1-4 & tot = Burrow volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³), Vcracks_L1-4 & tot = Cracks volume estimated from CT scan in the 0–8.5 cm layer (L1), 8.5–17 cm layer (L2), 17–25.5 cm layer (L3), 25.5–34 cm layer (L4), and in the whole mesocosm (tot) (cm³).

(TIF)

**S4 Fig. Effects of experimental treatments on several predictors of CO₂ and N₂O fluxes measured at the end of the experiment.** Green = in presence of plant, dark grey = in absence of plant. Different letters represent significantly different treatments according to Tukey's HSD post hoc test. Error bars represent ± 1 SEM.

(TIF)

**S5 Fig.** Cumulative N₂O (A) and CO₂ (B) emissions as affected by Plant×SWC interaction, where SWC represents the 3-month average SWC and the emissions the total cumulative gas emissions.

(TIF)

**S1 Table.** Effects earthworms (Ew) and plant treatments on the total macroporosity volume (pores) as well as differentiated as burrows and cracks (see Table 2 for detailed variable description). The “ns” abbreviation stands for variables that were not significant and were not retained in the minimal adequate models whereas $m_r^2$ represents the marginal coefficient of determination. ***P < 0.001; **P < 0.01; *P < 0.05; †P < 0.1.

(DOCX)

**Acknowledgments**

This study benefited from the CNRS human and technical resources allocated to the ECO-TRONS Research Infrastructure and AnaEE France. We thank Thierry Morvan for his assistance with the organization of the soil excavation at the EFELE experimental site as well as Thierry Mathieu and David Degueldre from the technical platform Terrain d'Expériences du C.E.F.E. for their support with the greenhouse and production of gas sampling collars. Microbial and soil analyses were done at Plateforme d'Analyses Chimiques en Ecologie (PACE) supported by the LabEx CeMEB, an ANR "Investissements d'avenir" programme (ANR-10-LABX-04-01).

**Author Contributions**

**Conceptualization:** Johanne Nahmani, Alexandru Milcu.

**Data curation:** Pierre Ganault, Yvan Capowiez, Alexandru Milcu.

**Formal analysis:** Pierre Ganault, Alexandru Milcu.

**Funding acquisition:** Johanne Nahmani, Alexandru Milcu.

**Investigation:** Pierre Ganault.
**Methodology:** Pierre Ganault, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, Isabelle Bertrand, Bruno Buatois, Alexandru Milcu.

**Validation:** Pierre Ganault.

**Visualization:** Pierre Ganault, Alexandru Milcu.

**Writing – original draft:** Pierre Ganault, Alexandru Milcu.

**Writing – review & editing:** Pierre Ganault, Alexandru Milcu.

## References

1. Schlesinger WH. Carbon Balance in Terrestrial Detritus. Annu Rev Ecol Syst. 1977; 8: 51–81. https://doi.org/10.1146/annurev.es.08.110177.000411

2. Wardle DA. Communities and ecosystems: linking the aboveground and belowground components. Princeton, N.J: Princeton University Press; 2002.

3. Lavelle P, Spain AV. Soil Ecology. New York, Boston, Dordrecht, London, Moscow: Kluwer Academic Publishers; 2001. Available: https://doi.org/10.1007/0-306-48162-6

4. Scheu S. Effects of earthworms on plant growth: patterns and perspectives. Pedobiologia. 2003; 47: 846–856. https://doi.org/10.1078/0031-4056-00270

5. Blouin M, Hodson ME, Delgado EA, Baker G, Brussaard L, Butt KR, et al. A review of earthworm impact on soil function and ecosystem services: Earthworm impact on ecosystem services. Eur J Soil Sci. 2013; 64: 161–182. https://doi.org/10.1111/ejss.12025

6. van Groenigen JW, Lubbers IM, Vos HMJ, Brown GG, De Deyn GB, van Groenigen KJ. Earthworms increase plant production: a meta-analysis. Sci Rep. 2015; 4: 6365. https://doi.org/10.1038/srep06365 PMID: 25219785

7. Lubbers IM, van Groenigen KJ, Fonte SJ, Six J, Brussaard L, van Groenigen JW. Greenhouse-gas emissions from soils increased by earthworms. Nature Clim Change. 2013; 3: 187–194. https://doi.org/10.1038/nclimate1692

8. Stocker T, editor. Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change. New York: Cambridge University Press; 2014.

9. Baggs EM. Soil microbial sources of nitrous oxide: recent advances in knowledge, emerging challenges and future direction. Current Opinion in Environmental Sustainability. 2011; 3: 321–327. https://doi.org/10.1016/j.cosust.2011.08.011

10. Firestone MK, Davidson EA. Microbiological Basis of NO and N2O production and consumption in soil. ohn Wiley&Sons Ltd., Exchange of Trace Gases between Terrestrial Ecosystems and the Atmosphere. ohn Wiley&Sons Ltd., 1989. pp. 7–21.

11. Bardgett RD, Freeman C, Ostle NJ. Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2008; 2: 805–814. https://doi.org/10.1038/ismej.2008.58 PMID: 18615117

12. Angst Š, Mueller CW, Cajthaml T, Angst G, Lhotáková Z, Bartuška M, et al. Stabilization of soil organic matter by earthworms is connected with physical protection rather than with chemical changes of organic matter. Geoderma. 2017; 289: 29–35. https://doi.org/10.1016/j.geoderma.2016.11.017

13. Cotrufo MF, Soong JL, Horton AJ, Campbell EE, Haddix ML, Wall DH, et al. Formation of soil organic matter via biochemical and physical pathways of litter mass loss. Nature Geosci. 2015; 8: 776–779. https://doi.org/10.1038/ngeo2520

14. Schmidt MWI, Torn MS, Abiven S, Dittmar T, Guggenberger G, Janssens IA, et al. Persistence of soil organic matter as an ecosystem property. Nature. 2011; 478: 49–56. https://doi.org/10.1038/nature10386 PMID: 21979045

15. Chen C, Whalen JK, Guo X. Earthworms reduce soil nitrous oxide emissions during drying and rewetting cycles. Soil Biology and Biochemistry. 2014; 68: 117–124. https://doi.org/10.1016/j.soilbio.2013.09.020

16. Kitzler B, Zechmeister-Boltenstern S, Holtermann C, Skiba U, Butterbach-Bahl K. Nitrogen oxides emission from two beech forests subjected to different nitrogen loads. Biogeosciences. 2006; 3: 293–310. https://doi.org/10.5194/bg-3-293-2006

17. Butterbach-Bahl K, Baggs EM, Dannenmann M, Kiese R, Zechmeister-Boltenstern S. Nitrous oxide emissions from soils: how well do we understand the processes and their controls? Philosophical Transactions of the Royal Society B: Biological Sciences. 2013; 368: 20130122–20130122. https://doi.org/10.1098/rstb.2013.0122 PMID: 23713120

18. Oertel C, Matschullat J, Zurba K, Zimmermann F, Erasmi S. Greenhouse gas emissions from soils—A review. Geochemistry. 2016; 76: 327–352. https://doi.org/10.1016/j.chemer.2016.04.002
**Methodology:** Pierre Ganault, Yvan Capowiez, Nathalie Fromin, Ammar Shihan, Isabelle Bertrand, Bruno Buatois, Alexandru Milcu.

**Validation:** Pierre Ganault.

**Visualization:** Pierre Ganault, Alexandru Milcu.

**Writing – original draft:** Pierre Ganault, Alexandru Milcu.

**Writing – review & editing:** Pierre Ganault, Alexandru Milcu.

## References

1. Schlesinger WH. Carbon Balance in Terrestrial Detritus. Annu Rev Ecol Syst. 1977; 8: 51–81. https://doi.org/10.1146/annurev.es.08.110177.000411

2. Wardle DA. Communities and ecosystems: linking the aboveground and belowground components. Princeton, N.J: Princeton University Press; 2002.

3. Lavelle P, Spain AV. Soil Ecology. New York, Boston, Dordrecht, London, Moscow: Kluwer Academic Publishers; 2001. Available: https://doi.org/10.1007/0-306-48162-6

4. Scheu S. Effects of earthworms on plant growth: patterns and perspectives. Pedobiologia. 2003; 47: 846–856. https://doi.org/10.1078/0031-4056-00270

5. Blouin M, Hodson ME, Delgado EA, Baker G, Brussaard L, Butt KR, et al. A review of earthworm impact on soil function and ecosystem services: Earthworm impact on ecosystem services. Eur J Soil Sci. 2013; 64: 161–182. https://doi.org/10.1111/ejss.12025

6. van Groenigen JW, Lubbers IM, Vos HMJ, Brown GG, De Deyn GB, van Groenigen KJ. Earthworms increase plant production: a meta-analysis. Sci Rep. 2015; 4: 6365. https://doi.org/10.1038/srep06365 PMID: 25219785

7. Lubbers IM, van Groenigen KJ, Fonte SJ, Six J, Brussaard L, van Groenigen JW. Greenhouse-gas emissions from soils increased by earthworms. Nature Clim Change. 2013; 3: 187–194. https://doi.org/10.1038/nclimate1692

8. Stocker T, editor. Climate change 2013: the physical science basis: Working Group I contribution to the Fifth assessment report of the Intergovernmental Panel on Climate Change. New York: Cambridge University Press; 2014.

9. Baggs EM. Soil microbial sources of nitrous oxide: recent advances in knowledge, emerging challenges and future direction. Current Opinion in Environmental Sustainability. 2011; 3: 321–327. https://doi.org/10.1016/j.cosust.2011.08.011

10. Firestone MK, Davidson EA. Microbiological Basis of NO and N2O production and consumption in soil. ohn Wiley&Sons Ltd., Exchange of Trace Gases between Terrestrial Ecosystems and the Atmosphere. ohn Wiley&Sons Ltd., 1989. pp. 7–21.

11. Bardgett RD, Freeman C, Ostle NJ. Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2008; 2: 805–814. https://doi.org/10.1038/ismej.2008.58 PMID: 18615117

12. Angst Š, Mueller CW, Cajthaml T, Angst G, Lhotáková Z, Bartuška M, et al. Stabilization of soil organic matter by earthworms is connected with physical protection rather than with chemical changes of organic matter. Geoderma. 2017; 289: 29–35. https://doi.org/10.1016/j.geoderma.2016.11.017

13. Cotrufo MF, Soong JL, Horton AJ, Campbell EE, Haddix ML, Wall DH, et al. Formation of soil organic matter via biochemical and physical pathways of litter mass loss. Nature Geosci. 2015; 8: 776–779. https://doi.org/10.1038/ngeo2520

14. Schmidt MWI, Torn MS, Abiven S, Dittmar T, Guggenberger G, Janssens IA, et al. Persistence of soil organic matter as an ecosystem property. Nature. 2011; 478: 49–56. https://doi.org/10.1038/nature10386 PMID: 21979045

15. Chen C, Whalen JK, Guo X. Earthworms reduce soil nitrous oxide emissions during drying and rewetting cycles. Soil Biology and Biochemistry. 2014; 68: 117–124. https://doi.org/10.1016/j.soilbio.2013.09.020

16. Kitzler B, Zechmeister-Boltenstern S, Holtermann C, Skiba U, Butterbach-Bahl K. Nitrogen oxides emission from two beech forests subjected to different nitrogen loads. Biogeosciences. 2006; 3: 293–310. https://doi.org/10.5194/bg-3-293-2006

17. Butterbach-Bahl K, Baggs EM, Dannenmann M, Kiese R, Zechmeister-Boltenstern S. Nitrous oxide emissions from soils: how well do we understand the processes and their controls? Philosophical Transactions of the Royal Society B: Biological Sciences. 2013; 368: 20130122–20130122. https://doi.org/10.1098/rstb.2013.0122 PMID: 23713120

18. Oertel C, Matschullat J, Zurba K, Zimmermann F, Erasmi S. Greenhouse gas emissions from soils—A review. Geochemistry. 2016; 76: 327–352. https://doi.org/10.1016/j.chemer.2016.04.002
19. Chapuis-Lardy L, Wrage N, Metay A, Chotte J-L, Bernoux M. Soils, a sink for N₂O? A review. Global Change Biology. 2007; 13: 1–17. https://doi.org/10.1111/j.1365-2486.2006.01280.x

20. Bateman EJ, Baggs EM. Contributions of nitrification and denitrification to N2O emissions from soils at different water-filled pore space. Biol Fertil Soils. 2005; 41: 379–388. https://doi.org/10.1007/s00374-005-0858-3

21. Seitzinger S, Harrison JA, Böhlke JK, Bouwman AF, Lowrance R, Peterson B, et al. Denitrification across landscapes and waterscapes: a synthesis. Ecological Applications. 2006; 16: 2064–2090. https://doi.org/10.1890/1051-0761(2006)016[2064:dalawa]2.0.co;2 PMID: 17205890

22. Davidson EA, Samanta S, Caramori SS, Savage K. The Dual Arrhenius and Michaelis-Menten kinetics model for decomposition of soil organic matter at hourly to seasonal time scales. Glob Change Biol. 2012; 18: 371–384. https://doi.org/10.1111/j.1365-2486.2011.02546.x

23. Skopp J, Jawson MD, Doran JW. Steady-State Aerobic Microbial Activity as a Function of Soil Water Content. Soil Science Society of America Journal. 1990; 54: 1619–1625. https://doi.org/10.2136/sssaj1990.03615995005400060018x

24. Weier KL, Doran JW, Power JF, Walters DT. Denitrification and the Dinitrogen/Nitrous Oxide Ratio as Affected by Soil Water, Available Carbon, and Nitrate. Soil Science Society of America Journal. 1993; 57: 66–72. https://doi.org/10.2136/sssaj1993.03615995005700010013x

25. Hénault C, Grossel A, Mary B, Roussel M, LéOnard J. Nitrous oxide emission by agricultural soils: a review of spatial and temporal variability for mitigation. Pedosphere. 2012; 22: 426–433. https://doi.org/10.1016/S1002-0160(12)60029-0

26. Decaëns T, Rangel AF, Asakawa N, Thomas RJ. Carbon and nitrogen dynamics in ageing earthworm casts in grasslands of the eastern plains of Colombia. Biology and Fertility of Soils. 1999; 30: 20–28. https://doi.org/10.1007/s003740050582

27. Drake HL, Horn MA. As the worm turns: the earthworm gut as a transient habitat for soil microbial biomes. Annu Rev Microbiol. 2007; 61: 169–189. https://doi.org/10.1146/annurev.micro.61.080706.093139 PMID: 17506687

28. Angst G, Mueller CW, Prater I, Angst Š, Frouz J, Jílková V, et al. Earthworms act as biochemical reactors to convert labile plant compounds into stabilized soil microbial necromass. Commun Biol. 2019; 2: 441. https://doi.org/10.1038/s42003-019-0684-z PMID: 31815196

29. Angst G, Frouz J, Van Groenigen JW, Scheu S, Kögel-Knabner I, Eisenhauer N. Earthworms as catalysts in the formation and stabilization of soil microbial necromass. Global Change Biology. 2022; 28: 4775–4782. https://doi.org/10.1111/gcb.16208 PMID: 35543252

30. Frouz J. Effects of soil macro- and mesofauna on litter decomposition and soil organic matter stabilization. Geoderma. 2018; 332: 161–172. https://doi.org/10.1016/j.geoderma.2017.08.039

31. Edwards CA, Arancon NQ. Biology and ecology of earthworms. Fourth edition. New York: Springer; 2022.

32. Chapin FS, Matson PA, Vitousek PM. Principles of Terrestrial Ecosystem Ecology. New York, NY: Springer New York; 2011. Available: http://link.springer.com/10.1007/978-1-4419-9504-9

33. Sokol NW, Kuebbing SaraE, Karlsen-Ayala E, Bradford MA. Evidence for the primacy of living root inputs, not root or shoot litter, in forming soil organic carbon. New Phytologist. 2019; 221: 233–246. https://doi.org/10.1111/nph.15361 PMID: 30067293

34. Cotrufo MF, Wallenstein MD, Boot CM, Denef K, Paul E. The Microbial Efficiency-Matrix Stabilization (MEMS) framework integrates plant litter decomposition with soil organic matter stabilization: do labile plant inputs form stable soil organic matter? Glob Change Biol. 2013; 19: 988–995. https://doi.org/10.1111/gcb.12113 PMID: 23504877

35. Hinsinger P, Plassard C, Jaillard B. Rhizosphere: A new frontier for soil biogeochemistry. Journal of Geochemical Exploration. 2006; 88: 210–213. https://doi.org/10.1016/j.gexplo.2005.08.041

36. Angers DA, Caron J. Plant-induced changes in soil structure: processes and feedbacks. Biogeochemistry. 1998; 42: 55–72. https://doi.org/10.1023/A:1005944025343

37. Hodge A, Robinson D, Fitter A. Are microorganisms more effective than plants at competing for nitrogen? Trends in Plant Science. 2000; 5: 304–308. https://doi.org/10.1016/s1360-1385(00)01656-3 PMID: 10871903

38. Moreau D, Bardgett RD, Finlay RD, Jones DL, Philippot L. A plant perspective on nitrogen cycling in the rhizosphere. Power S, editor. Funct Ecol. 2019; 33: 540–552. https://doi.org/10.1111/1365-2435.13303

39. Capowiez Y, Sammartino S, Michel E. Burrow systems of endogeic earthworms: Effects of earthworm abundance and consequences for soil water infiltration. Pedobiologia. 2014; 57: 303–309. https://doi.org/10.1016/j.pedobi.2014.04.001

40. Capowiez Y, Bottinelli N, Sammartino S, Michel E, Jouquet P. Morphological and functional characterization of the burrow systems of six earthworm species (Lumbricidae). Biology and Fertility of Soils. 2015; 51: 869–877. https://doi.org/10.1007/s00374-015-1036-x
19. Chapuis-Lardy L, Wrage N, Metay A, Chotte J-L, Bernoux M. Soils, a sink for N₂O? A review. Global Change Biology. 2007; 13: 1–17. https://doi.org/10.1111/j.1365-2486.2006.01280.x

20. Bateman EJ, Baggs EM. Contributions of nitrification and denitrification to N2O emissions from soils at different water-filled pore space. Biol Fertil Soils. 2005; 41: 379–388. https://doi.org/10.1007/s00374-005-0858-3

21. Seitzinger S, Harrison JA, Böhlke JK, Bouwman AF, Lowrance R, Peterson B, et al. Denitrification across landscapes and waterscapes: a synthesis. Ecological Applications. 2006; 16: 2064–2090. https://doi.org/10.1890/1051-0761(2006)016[2064:dalawa]2.0.co;2 PMID: 17205890

22. Davidson EA, Samanta S, Caramori SS, Savage K. The Dual Arrhenius and Michaelis-Menten kinetics model for decomposition of soil organic matter at hourly to seasonal time scales. Glob Change Biol. 2012; 18: 371–384. https://doi.org/10.1111/j.1365-2486.2011.02546.x

23. Skopp J, Jawson MD, Doran JW. Steady-State Aerobic Microbial Activity as a Function of Soil Water Content. Soil Science Society of America Journal. 1990; 54: 1619–1625. https://doi.org/10.2136/sssaj1990.03615995005400060018x

24. Weier KL, Doran JW, Power JF, Walters DT. Denitrification and the Dinitrogen/Nitrous Oxide Ratio as Affected by Soil Water, Available Carbon, and Nitrate. Soil Science Society of America Journal. 1993; 57: 66–72. https://doi.org/10.2136/sssaj1993.03615995005700010013x

25. Hénault C, Grossel A, Mary B, Roussel M, LéOnard J. Nitrous oxide emission by agricultural soils: a review of spatial and temporal variability for mitigation. Pedosphere. 2012; 22: 426–433. https://doi.org/10.1016/S1002-0160(12)60029-0

26. Decaëns T, Rangel AF, Asakawa N, Thomas RJ. Carbon and nitrogen dynamics in ageing earthworm casts in grasslands of the eastern plains of Colombia. Biology and Fertility of Soils. 1999; 30: 20–28. https://doi.org/10.1007/s003740050582

27. Drake HL, Horn MA. As the worm turns: the earthworm gut as a transient habitat for soil microbial biomes. Annu Rev Microbiol. 2007; 61: 169–189. https://doi.org/10.1146/annurev.micro.61.080706.093139 PMID: 17506687

28. Angst G, Mueller CW, Prater I, Angst Š, Frouz J, Jílková V, et al. Earthworms act as biochemical reactors to convert labile plant compounds into stabilized soil microbial necromass. Commun Biol. 2019; 2: 441. https://doi.org/10.1038/s42003-019-0684-z PMID: 31815196

29. Angst G, Frouz J, Van Groenigen JW, Scheu S, Kögel-Knabner I, Eisenhauer N. Earthworms as catalysts in the formation and stabilization of soil microbial necromass. Global Change Biology. 2022; 28: 4775–4782. https://doi.org/10.1111/gcb.16208 PMID: 35543252

30. Frouz J. Effects of soil macro- and mesofauna on litter decomposition and soil organic matter stabilization. Geoderma. 2018; 332: 161–172. https://doi.org/10.1016/j.geoderma.2017.08.039

31. Edwards CA, Arancon NQ. Biology and ecology of earthworms. Fourth edition. New York: Springer; 2022.

32. Chapin FS, Matson PA, Vitousek PM. Principles of Terrestrial Ecosystem Ecology. New York, NY: Springer New York; 2011. Available: http://link.springer.com/10.1007/978-1-4419-9504-9

33. Sokol NW, Kuebbing SaraE, Karlsen-Ayala E, Bradford MA. Evidence for the primacy of living root inputs, not root or shoot litter, in forming soil organic carbon. New Phytologist. 2019; 221: 233–246. https://doi.org/10.1111/nph.15361 PMID: 30067293

34. Cotrufo MF, Wallenstein MD, Boot CM, Denef K, Paul E. The Microbial Efficiency-Matrix Stabilization (MEMS) framework integrates plant litter decomposition with soil organic matter stabilization: do labile plant inputs form stable soil organic matter? Glob Change Biol. 2013; 19: 988–995. https://doi.org/10.1111/gcb.12113 PMID: 23504877

35. Hinsinger P, Plassard C, Jaillard B. Rhizosphere: A new frontier for soil biogeochemistry. Journal of Geochemical Exploration. 2006; 88: 210–213. https://doi.org/10.1016/j.gexplo.2005.08.041

36. Angers DA, Caron J. Plant-induced changes in soil structure: processes and feedbacks. Biogeochemistry. 1998; 42: 55–72. https://doi.org/10.1023/A:1005944025343

37. Hodge A, Robinson D, Fitter A. Are microorganisms more effective than plants at competing for nitrogen? Trends in Plant Science. 2000; 5: 304–308. https://doi.org/10.1016/s1360-1385(00)01656-3 PMID: 10871903

38. Moreau D, Bardgett RD, Finlay RD, Jones DL, Philippot L. A plant perspective on nitrogen cycling in the rhizosphere. Power S, editor. Funct Ecol. 2019; 33: 540–552. https://doi.org/10.1111/1365-2435.13303

39. Capowiez Y, Sammartino S, Michel E. Burrow systems of endogeic earthworms: Effects of earthworm abundance and consequences for soil water infiltration. Pedobiologia. 2014; 57: 303–309. https://doi.org/10.1016/j.pedobi.2014.04.001

40. Capowiez Y, Bottinelli N, Sammartino S, Michel E, Jouquet P. Morphological and functional characterization of the burrow systems of six earthworm species (Lumbricidae). Biology and Fertility of Soils. 2015; 51: 869–877. https://doi.org/10.1007/s00374-015-1036-x
41. Cheik S, Jouquet P, Maeght J, Capowiez Y, Tran TM, Bottinelli N. X-ray tomography analysis of soil bio-pores structure under wetting and drying cycles. Eur J Soil Sci. 2021; 72: 2128-2132. https://doi.org/10.1111/ejss.13119

42. Ernst G, Felten D, Vohland M, Emmerling C. Impact of ecologically different earthworm species on soil water characteristics. European Journal of Soil Biology. 2009; 45: 207-213. https://doi.org/10.1016/j.ejsobi.2009.01.001

43. Law MMS, Lai DYF. Impacts of wetting-drying cycles on short-term carbon and nitrogen dynamics in Amynthas earthworm casts. Pedosphere. 2021; 31: 423-432. https://doi.org/10.1016/S1002-0160(20)60081-9

44. Bottinelli N, Hedde M, Jouquet P, Capowiez Y. An explicit definition of earthworm ecological categories-Marcel Bouché's triangle revisited. Geoderma. 2020; 372: 114361. https://doi.org/10.1016/j.geoderma.2020.114361

45. Bouché M. Lombriciens de France. Ecologie et systématique. Paris: INRA; 1972.

46. Hallam J, Hodson M. Earthworms and soil water regulation: A review. African and Mediterranean Agricultural Journal - Al Awamia. 2021; 130: 89-135. https://doi.org/10.34874/IMIST.PRSM/afrimed-i130.31381

47. Le Mer G, Bottinelli N, Dignac M-F, Capowiez Y, Jouquet P, Mazurier A, et al. Exploring the control of earthworm cast macro- and micro-scale features on soil organic carbon mineralization across species and ecological categories. Geoderma. 2022; 427: 116151. https://doi.org/10.1016/j.geoderma.2022.116151

48. Hong S-Y, Park J-H, Cho S-H, Yang M-S, Park C-M. Phenological growth stages of Brachypodium distachyon: codification and description: BBCH Brachypodium distachyon growth stages. Weed Research. 2011; 51: 612-620. https://doi.org/10.1111/j.1365-3180.2011.00877.x

49. Brkljacic J, Grotewold E, Scholl R, Mockler T, Garvin DF, Vain P, et al. Brachypodium as a Model for the Grasses: Today and the Future. Plant Physiology. 2011; 157: 3-13. https://doi.org/10.1104/pp.111.179531 PMID: 21771916

50. Rochette P. Towards a standard non-steady-state chamber methodology for measuring soil N2O emissions. Animal Feed Science and Technology. 2011;166-167: 141-146. https://doi.org/10.1016/j.anifeedsci.2011.04.063

51. Capowiez Y, Bastardie F, Costagliola G. Sublethal effects of imidacloprid on the burrowing behaviour of two earthworm species: Modifications of the 3D burrow systems in artificial cores and consequences on gas diffusion in soil. Soil Biology and Biochemistry. 2006; 38: 285-293. https://doi.org/10.1016/j.soilbio.2005.05.014

52. Capowiez Y, Pierret A, Daniel O, Monestiez P, Kretzschmar A. 3D skeleton reconstructions of natural earthworm burrow systems using CAT scan images of soil cores. Biology and Fertility of Soils. 1998; 27: 51-59. https://doi.org/10.1007/s003740050399

53. Anderson EW. A Guide for Estimating Cover. Rangelands. 1986; 8: 236-238.

54. Daubenmire RF. Canopy coverage method of vegetation analysis. Northwest Sci. 1959; 33: 39-64.

55. Stohlgren TJ, Bull KA, Otsuki Y. Comparison of rangeland vegetation sampling techniques in the Central Grasslands. 51. 1998 [cited 8 Nov 2023]. https://doi.org/10.2307/4003202

56. Lensi R, Gourbiere F, Josserand A. Measurement of small amounts of nitrate in an acid soil by N2O production. Soil Biology and Biochemistry. 1985; 17: 733-734. https://doi.org/10.1016/0038-0717(85)90056-2

57. Campbell CD, Chapman SJ, Cameron CM, Davidson MS, Potts JM. A Rapid Microtiter Plate Method To Measure Carbon Dioxide Evolved from Carbon Substrate Amendments so as To Determine the Physiological Profiles of Soil Microbial Communities by Using Whole Soil. Appl Environ Microbiol. 2003; 69: 3593-3599. https://doi.org/10.1128/AEM.69.6.3593-3599.2003 PMID: 12788767

58. Anderson JPE, Domsch KH. A physiological method for the quantitative measurement of microbial biomass in soils. Soil Biology and Biochemistry. 1978; 10: 215-221. https://doi.org/10.1016/0038-0717(78)90099-8

59. Pinheiro JC, Bates DM, DebRoy SS, Sarkar D. Nlme: Linear and Nonlinear Mixed Effects Models. CRAN. 2013.

60. Zuur AF, Ieno EN, Walker N, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R. New York, NY: Springer New York; 2009. https://doi.org/10.1007/978-0-387-87458-6

61. Barton K. Package "MuMin." CRAN. 2020.

62. Breheny P, Huang J. Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann Appl Stat. 2011; 5. https://doi.org/10.1214/10-AOAS388 PMID: 22081779
41. Cheik S, Jouquet P, Maeght J, Capowiez Y, Tran TM, Bottinelli N. X-ray tomography analysis of soil bio-pores structure under wetting and drying cycles. Eur J Soil Sci. 2021; 72: 2128–2132. https://doi.org/10.1111/ejss.13119

42. Ernst G, Felten D, Vohland M, Emmerling C. Impact of ecologically different earthworm species on soil water characteristics. European Journal of Soil Biology. 2009; 45: 207–213. https://doi.org/10.1016/j.ejsobi.2009.01.001

43. Law MMS, Lai DYF. Impacts of wetting-drying cycles on short-term carbon and nitrogen dynamics in Amynthas earthworm casts. Pedosphere. 2021; 31: 423–432. https://doi.org/10.1016/S1002-0160(20)60081-9

44. Bottinelli N, Hedde M, Jouquet P, Capowiez Y. An explicit definition of earthworm ecological categories–Marcel Bouché’s triangle revisited. Geoderma. 2020; 372: 114361. https://doi.org/10.1016/j.geoderma.2020.114361

45. Bouché M. Lombriciens de France. Ecologie et systématique. Paris: INRA; 1972.

46. Hallam J, Hodson M. Earthworms and soil water regulation: A review. African and Mediterranean Agricultural Journal - Al Awamia. 2021; 130: 89–135. https://doi.org/10.34874/IMIST.PRSM/afrimed-i130.31381

47. Le Mer G, Bottinelli N, Dignac M- F, Capowiez Y, Jouquet P, Mazurier A, et al. Exploring the control of earthworm cast macro- and micro-scale features on soil organic carbon mineralization across species and ecological categories. Geoderma. 2022; 427: 116151. https://doi.org/10.1016/j.geoderma.2022.116151

48. Hong S-Y, Park J-H, Cho S-H, Yang M-S, Park C-M. Phenological growth stages of Brachypodium distachyon: codification and description: BBCH Brachypodium distachyon growth stages. Weed Research. 2011; 51: 612–620. https://doi.org/10.1111/j.1365-3180.2011.00877.x

49. Brkljacic J, Grotewold E, Scholl R, Mockler T, Garvin DF, Vain P, et al. Brachypodium as a Model for the Grasses: Today and the Future. Plant Physiology. 2011; 157: 3–13. https://doi.org/10.1104/pp.111.179531 PMID: 21771916

50. Rochette P. Towards a standard non-steady-state chamber methodology for measuring soil N2O emissions. Animal Feed Science and Technology. 2011;166–167: 141–146. https://doi.org/10.1016/j.anifeedsci.2011.04.063

51. Capowiez Y, Bastardie F, Costagliola G. Sublethal effects of imidacloprid on the burrowing behaviour of two earthworm species: Modifications of the 3D burrow systems in artificial cores and consequences on gas diffusion in soil. Soil Biology and Biochemistry. 2006; 38: 285–293. https://doi.org/10.1016/j.soilbio.2005.05.014

52. Capowiez Y, Pierret A, Daniel O, Monestiez P, Kretzschmar A. 3D skeleton reconstructions of natural earthworm burrow systems using CAT scan images of soil cores. Biology and Fertility of Soils. 1998; 27: 51–59. https://doi.org/10.1007/s003740050399

53. Anderson EW. A Guide for Estimating Cover. Rangelands. 1986; 8: 236–238.

54. Daubenmire RF. Canopy coverage method of vegetation analysis. Northwest Sci. 1959; 33: 39–64.

55. Stohlgren TJ, Bull KA, Otsuki Y. Comparison of rangeland vegetation sampling techniques in the Central Grasslands. 51. 1998 [cited 8 Nov 2023]. https://doi.org/10.2307/4003202

56. Lensi R, Gourbiere F, Josserand A. Measurement of small amounts of nitrate in an acid soil by N2O production. Soil Biology and Biochemistry. 1985; 17: 733–734. https://doi.org/10.1016/0038-0717(85)90056-2

57. Campbell CD, Chapman SJ, Cameron CM, Davidson MS, Potts JM. A Rapid Microtiter Plate Method To Measure Carbon Dioxide Evolved from Carbon Substrate Amendments so as To Determine the Physiological Profiles of Soil Microbial Communities by Using Whole Soil. Appl Environ Microbiol. 2003; 69: 3593–3599. https://doi.org/10.1128/AEM.69.6.3593-3599.2003 PMID: 12788767

58. Anderson JPE, Domsch KH. A physiological method for the quantitative measurement of microbial biomass in soils. Soil Biology and Biochemistry. 1978; 10: 215–221. https://doi.org/10.1016/0038-0717(78)90099-8

59. Pinheiro JC, Bates DM, DebRoy SS, Sarkar D. Nlme: Linear and Nonlinear Mixed Effects Models. CRAN. 2013.

60. Zuur AF, Ieno EN, Walker N, Saveliev AA, Smith GM. Mixed effects models and extensions in ecology with R. New York, NY: Springer New York; 2009. https://doi.org/10.1007/978-0-387-87458-6

61. Barton K. Package "MuMin." CRAN. 2020.

62. Breheny P, Huang J. Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann Appl Stat. 2011; 5. https://doi.org/10.1214/10-AOAS388 PMID: 22081779
63. Breheny P. Package 'ncvreg': regularization paths for SCAD and MCP penalized regression models. 2020.

64. Robertson GP, Tiedje JM. Nitrous oxide sources in aerobic soils: Nitrification, denitrification and other biological processes. Soil Biology and Biochemistry. 1987; 19: 187-193. https://doi.org/10.1016/0038-0717(87)90080-0

65. Davidson EA, Keller M, Erickson HE, Verchot LV, Veldkamp E. Testing a Conceptual Model of Soil Emissions of Nitrous and Nitric Oxides. BioScience. 2000; 50: 667. https://doi.org/10.1641/0006-3568 (2000)050[0667:TACMOS]2.0.CO;2

66. Bernier N. Earthworm feeding activity and development of the humus profile. Biology and Fertility of Soils. 1998; 26: 215-223. https://doi.org/10.1007/s003740050370

67. Lubbers IM, Berg MP, De Deyn GB, Putten WH, Groenigen JW. Soil fauna diversity increases CO₂ but suppresses N₂O emissions from soil. Glob Change Biol. 2020; 26: 1886-1898. https://doi.org/10.1111/gcb.14860 PMID: 31587448

68. Kuiper I, de Deyn GB, Thakur MP, van Groenigen JW. Soil invertebrate fauna affect N₂ O emissions from soil. Glob Change Biol. 2013; 19: 2814-2825. https://doi.org/10.1111/gcb.12232 PMID: 23625707

69. Six J, Bossuyt H, Degryze S, Denef K. A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics. Soil and Tillage Research. 2004; 79: 7-31. https://doi.org/10.1016/j.still.2004.03.008

70. Erktan A, Cécillon L, Graf F, Roumet C, Legout C, Rey F. Increase in soil aggregate stability along a Mediterranean successional gradient in severely eroded gully bed ecosystems: combined effects of soil, root traits and plant community characteristics. Plant Soil. 2016; 398: 121-137. https://doi.org/10.1007/s11104-015-2647-6

71. Barnard R, Leadley PW, Hungate BA. Global change, nitrification, and denitrification: A review. Global Biogeochem Cycles. 2005; 19. https://doi.org/10.1029/2004GB002282

72. Van Groenigen JW, Van Groenigen KJ, Koopmans GF, Stokkermans L, Vos HMJ, Lubbers IM. How fer- tile are earthworm casts? A meta-analysis. Geoderma. 2019; 338: 525-535. https://doi.org/10.1016/j.geoderma.2018.11.001

73. De Boeck HJ, Vicca S, Roy J, Nijs I, Milcu A, Kreyling J, et al. Global Change Experiments: Challenges and Opportunities. BioScience. 2015; 65: 922-931. https://doi.org/10.1093/biosci/biv099

74. Stewart RIA, Dossena M, Bohan DA, Jeppesen E, Kordas RL, Ledger ME, et al. Mesocosm Experi- ments as a Tool for Ecological Climate-Change Research. Advances in Ecological Research. Elsevier; 2013. pp. 71-181. https://doi.org/10.1016/B978-0-12-417199-2.00002-1

75. Schindlbacher A. Effects of soil moisture and temperature on NO, NO₂, and N₂ O emissions from Euro- pean forest soils. J Geophys Res. 2004; 109: D17302. https://doi.org/10.1029/2004JD004590

76. Clause J, Barot S, Richard B, Decaëns T, Forey E. The interactions between soil type and earthworm species determinethe properties of earthworm casts. Applied Soil Ecology. 2014; 83: 149-158. https:// doi.org/10.1016/j.apsoil.2013.12.006

77. Bottinelli N, Capowiez Y. Earthworm ecological categories are not functional groups. Biol Fertil Soils. 2021; 57: 329-331. https://doi.org/10.1007/s00374-020-01517-1

78. Smith P, Martino D, Cai Z, Gwary D, Janzen H, Kumar P, et al. Greenhouse gas mitigation in agricul- ture. Philosophical Transactions of the Royal Society B: Biological Sciences. 2007; 363: 789-813. https://doi.org/10.1098/rstb.2007.2184 PMID: 17827109

79. Briones MJI, Schmidt O. Conventional tillage decreases the abundance and biomass of earthworms and alters their community structure in a global meta-analysis. Glob Change Biol. 2017; 23: 4396-4419. https://doi.org/10.1111/gcb.13744 PMID: 28464547

80. Capowiez Y, Sammartino S, Keller T, Bottinelli N. Decreased burrowing activity of endogeic earthworms and effects on water infiltration in response to an increase in soil bulk density. Pedobiologia. 2021;85- 86: 150728. https://doi.org/10.1016/j.pedobi.2021.150728

81. Lipiec J, Kuś J, Słowińska-Jurkiewicz A, Nosalewicz A. Soil porosity and water infiltration as influenced by tillage methods. Soil and Tillage Research. 2006; 89: 210-220. https://doi.org/10.1016/j.still.2005.07.012

82. Pérès G, Bellido A, Curmi P, Marmonier P, Cluzeau D. Relationships between earthworm communities and burrow numbers under different land use systems. Pedobiologia. 2010; 54: 37-44. https://doi.org/10.1016/j.pedobi.2010.08.006
63. Breheny P. Package 'ncvreg': regularization paths for SCAD and MCP penalized regression models. 2020.

64. Robertson GP, Tiedje JM. Nitrous oxide sources in aerobic soils: Nitrification, denitrification and other biological processes. Soil Biology and Biochemistry. 1987; 19: 187-193. https://doi.org/10.1016/0038-0717(87)90080-0

65. Davidson EA, Keller M, Erickson HE, Verchot LV, Veldkamp E. Testing a Conceptual Model of Soil Emissions of Nitrous and Nitric Oxides. BioScience. 2000; 50: 667. https://doi.org/10.1641/0006-3568 (2000)050[0667:TACMOS]2.0.CO;2

66. Bernier N. Earthworm feeding activity and development of the humus profile. Biology and Fertility of Soils. 1998; 26: 215-223. https://doi.org/10.1007/s003740050370

67. Lubbers IM, Berg MP, De Deyn GB, Putten WH, Groenigen JW. Soil fauna diversity increases CO₂ but suppresses N₂O emissions from soil. Glob Change Biol. 2020; 26: 1886-1898. https://doi.org/10.1111/gcb.14860 PMID: 31587448

68. Kuiper I, de Deyn GB, Thakur MP, van Groenigen JW. Soil invertebrate fauna affect N₂ O emissions from soil. Glob Change Biol. 2013; 19: 2814-2825. https://doi.org/10.1111/gcb.12232 PMID: 23625707

69. Six J, Bossuyt H, Degryze S, Denef K. A history of research on the link between (micro)aggregates, soil biota, and soil organic matter dynamics. Soil and Tillage Research. 2004; 79: 7-31. https://doi.org/10.1016/j.still.2004.03.008

70. Erktan A, Cécillon L, Graf F, Roumet C, Legout C, Rey F. Increase in soil aggregate stability along a Mediterranean successional gradient in severely eroded gully bed ecosystems: combined effects of soil, root traits and plant community characteristics. Plant Soil. 2016; 398: 121-137. https://doi.org/10.1007/s11104-015-2647-6

71. Barnard R, Leadley PW, Hungate BA. Global change, nitrification, and denitrification: A review. Global Biogeochem Cycles. 2005; 19. https://doi.org/10.1029/2004GB002282

72. Van Groenigen JW, Van Groenigen KJ, Koopmans GF, Stokkermans L, Vos HMJ, Lubbers IM. How fer- tile are earthworm casts? A meta-analysis. Geoderma. 2019; 338: 525-535. https://doi.org/10.1016/j.geoderma.2018.11.001

73. De Boeck HJ, Vicca S, Roy J, Nijs I, Milcu A, Kreyling J, et al. Global Change Experiments: Challenges and Opportunities. BioScience. 2015; 65: 922-931. https://doi.org/10.1093/biosci/biv099

74. Stewart RIA, Dossena M, Bohan DA, Jeppesen E, Kordas RL, Ledger ME, et al. Mesocosm Experi- ments as a Tool for Ecological Climate-Change Research. Advances in Ecological Research. Elsevier; 2013. pp. 71-181. https://doi.org/10.1016/B978-0-12-417199-2.00002-1

75. Schindlbacher A. Effects of soil moisture and temperature on NO, NO₂, and N₂ O emissions from Euro- pean forest soils. J Geophys Res. 2004; 109: D17302. https://doi.org/10.1029/2004JD004590

76. Clause J, Barot S, Richard B, Decaëns T, Forey E. The interactions between soil type and earthworm species determinethe properties of earthworm casts. Applied Soil Ecology. 2014; 83: 149-158. https:// doi.org/10.1016/j.apsoil.2013.12.006

77. Bottinelli N, Capowiez Y. Earthworm ecological categories are not functional groups. Biol Fertil Soils. 2021; 57: 329-331. https://doi.org/10.1007/s00374-020-01517-1

78. Smith P, Martino D, Cai Z, Gwary D, Janzen H, Kumar P, et al. Greenhouse gas mitigation in agricul- ture. Philosophical Transactions of the Royal Society B: Biological Sciences. 2007; 363: 789-813. https://doi.org/10.1098/rstb.2007.2184 PMID: 17827109

79. Briones MJI, Schmidt O. Conventional tillage decreases the abundance and biomass of earthworms and alters their community structure in a global meta-analysis. Glob Change Biol. 2017; 23: 4396-4419. https://doi.org/10.1111/gcb.13744 PMID: 28464547

80. Capowiez Y, Sammartino S, Keller T, Bottinelli N. Decreased burrowing activity of endogeic earthworms and effects on water infiltration in response to an increase in soil bulk density. Pedobiologia. 2021;85- 86: 150728. https://doi.org/10.1016/j.pedobi.2021.150728

81. Lipiec J, Kuś J, Słowińska-Jurkiewicz A, Nosalewicz A. Soil porosity and water infiltration as influenced by tillage methods. Soil and Tillage Research. 2006; 89: 210-220. https://doi.org/10.1016/j.still.2005.07.012

82. Pérès G, Bellido A, Curmi P, Marmonier P, Cluzeau D. Relationships between earthworm communities and burrow numbers under different land use systems. Pedobiologia. 2010; 54: 37-44. https://doi.org/10.1016/j.pedobi.2010.08.006