The experiment was conducted at the Montpellier European Ecotron (CNRS, Montpellier, France) in the Mesocosms experimental platform,¹ whose capabilities were briefly described by Roy et al. (2021). The experiment included a treatment of functional diversity of soil macrofauna communities (LowFD vs. HighFD), factorially crossed with a drought treatment (Control with typical summer drought vs. Drought including a more severe and prolonged summer drought), with four replicated mesocosms per treatment combination (2 levels of functional diversity of fauna × 2 levels of summer drought × 4 replicates = 16 mesocosms; Supplementary Table 1). Each mesocosm unit consisted of a belowground compartment containing a 1 m³ stainless steel lysimeter with a 1 m² surface and a 3.9 m³ aboveground compartment enclosed by a highly transparent material permeable to light and UV radiation (250 μm thick Teflon-FEP film, DuPont, USA).

The soil was excavated from a Mediterranean old-field (previously a Mediterranean forest) located at Monferrier-sur-Lez (43°40′52.8″N, 3°52′34.3″E) in January 2018, after removing the herbaceous plants. Three soil layers were excavated separately: the first 10 cm (L1), 10–30 cm (L2) and 30–80 cm (L3). Larger rocks were manually removed from each soil layer. The first two layers, L1 and L2 were then defaunated by gamma irradiation at Synergy Health Marseille SAS facilities, Chusclan, France. Before and after sterilization, sub-samples of soil were extracted to measure soil chemical properties and microbial communities. Compared to other defaunation methods such as autoclaving, or freezing and thawing cycles, the gamma irradiation offers two advantages: (1) a lesser effect on soil physicochemical properties (Berns et al., 2008) and (2) the possibility of selectively target different organism groups with varying doses. We chose a relatively low dose of irradiation (2.5 kGy) to preferentially eliminate macrofauna species while leaving at least part of the microbial community undamaged (McNamara et al., 2003). For sterilization, the soil was bagged in batches of 1 ton each to reach 2.5 kGy at the core of each batch; the peripheral areas received higher doses up to 10 kGy.

The soil was then added by layers to each of the 16 lysimeters, reconstructing the natural layers and bulk density. Suction cups were installed at the bottom of each lysimeter to control water potential. Before adding the soil, we covered the suction cups with a 10 cm layer of gravel, which also facilitated drainage. Layer L3 filled each lysimeter between −90 cm and −30 cm, layer L2 between −30 cm and −10 cm, and finally, layer L1 to the top 10 cm. During the filling of each lysimeter, we added TDR sensors for measurements of volumetric soil moisture content and temperature at three different depths (−60 cm, −30 cm and −10 cm; two sensors per depth).

Each lysimeter was then planted with two individuals of two-year-old saplings of each of four different tree species in March 2018: two evergreen species, Quercus ilex L. 1753 and Arbutus unedo L. 1753 and two deciduous species, Quercus pubescens Willd. 1796 and Acer monspessulanum L. 1753. These tree species all co-occur in Mediterranean forests of southern France and differ with regard to their leaf litter quality as well as their mycorrhizal fungi (ectomycorrhiza: the two-oak species, arbuscular mycorrhiza: A. monspessulanum, and A. unedo). Before planting, the root systems of each sapling were carefully washed to remove soil. The saplings were planted equidistantly spaced to allow similar light exposure. Due to high mortality of Q. ilex before the beginning of the experiment, new saplings grown from acorns collected at the Puéchabon forest were planted in February 2019.

We then reconstructed the litter layer in each of the 16 mesocosms in July 2018. The litter layer was composed of equal amounts of all four tree species with a total quantity of 62 ± 0.1 g m⁻². The leaf litter was collected in surrounding natural forests in fall 2017 (for the two deciduous species) and in June 2018 (for the two evergreen species) and dried at 60°C for 72 h to kill all invertebrates. We later added a layer of manually fragmented leaf litter from all four tree species between the soil surface and the litter layer consisting of intact leaves in June 2019, to improve habitat conditions for the treatment-specific soil fauna (see below).

To ensure that the soil and litter layer contained a microbial and microfauna community representative of the targeted model system (a Mediterranean forest ecosystem), we added an inoculum to each mesocosm created from soils of two forest sites: Montarnaud-France (March 2018) and Puéchabon, France (December 2018). The inoculum was extracted from the top 10 cm of natural forest soil including three of the target species Q. ilex, Q. pubescens and A. unedo. We mixed one volume of soil to five volumes of deionized water, manually blended for 30 min. The solution was passed through a 68 μm sieve, and 5 L per mesocosm were homogeneously added. One week after the second inoculation, sub-samples of soil were extracted to measure microbial communities in the “initial conditions.”

The simulated climate conditions, with the exception of the drought treatment, were identical across all mesocosms during the whole experiment and were based on climatic records from 2013 at the nearby experimental site of Puéchabon, a well-studied Mediterranean forest (Lempereur et al., 2015). The year 2013 was selected as it represented an average year for temperature and precipitation during the 2003–2017 period.

During the summer drought period, the mesocosm watering regime deviated from the 2013 Puéchabon precipitation scenario due to the need to determine the appropriate reduction level for inducing realistic water stress under controlled conditions. A tensiometer-guided approach was employed, regulating drought cycles based on soil water potential at 30 cm depth. These cycles were repeated three times annually from early summer to autumn. Soil water potential was monitored in a reference mesocosm per drought level using a full-range tensiometer (FR T15D, UGT, Germany) to define drought onset and duration.

In 2019, a milder Drought treatment was applied to accommodate sapling root establishment. The Control treatment simulated Mediterranean summer drought, with drought cycles spanning the decline of soil water potential from −0.2 MPa to −1 MPa (moderate to severe water stress). Upon reaching −1 MPa, all eight Control mesocosms were rewatered with the cumulative water lost via evapotranspiration, measured by lysimeter weight changes. In the Drought treatment, drought duration was extended by 30% in the first two cycles and doubled in the third cycle, causing a temporary asynchrony between treatments, with the Control mesocosms receiving water earlier than the Drought mesocosms after each cycle. Consequently, Drought mesocosms received 23% less water than Control between June and October 2019.

In 2020, a similar approach was applied, but water stress in the Drought treatment was intensified by extending drought duration by 30% and reducing precipitation input by 30% per cycle relative to Control. The threshold for defining the drought period was lowered, with soil water potential reaching −1.5 MPa (permanent wilting point) before rewatering. These adjustments resulted in a 38% reduction in precipitation for the Drought treatment between June and October 2020. In both years, following three drought cycles, Control and Drought mesocosms were rewatered according to the Puéchabon precipitation scenario (Supplementary Table 2).

We took soil samples before the drought cycles (“pre-drought” in June of each year, 2019 and 2020) and at the end of the three drought cycles (“post-drought” in September/October of each year) (Figure 1). For each sampling date, three soil cores were randomly collected from the L1 layer (i.e., top soil down to 10 cm) of each mesocosm. The soil was homogenized by sieving through a 4 mm mesh. Three subsamples of 10 g for each mesocosms were taken from the sieved bulk soil and frozen at −80°C for DNA analysis.

We manipulated soil macrofauna species composition to create two levels of functional diversity, low functional diversity (LowFD) vs. high functional diversity (HighFD), while keeping species richness constant at three species per mesocosm. Twelve sapro-geophageous species common in Mediterranean forest were used: three isopod species, Armadillidium vulgare Latreille, 1804 (Armadillidae), Armadillo officinalis Duméril, 1816 (Armadillidae), Porcellio laevis Latreille, 1804 (Porcellionidae); three millipedes (Myriapoda: Diplopoda) Glomeris marginata Villers, 1789 (Glomeridae), Cylindroiulus caeruleocinctus Wood, 1864 (Julidae), and Ommatoiulus sabulosus (Linnaeus, 1758) (Julidae); three endogeic earthworms (Annelida: Lumbricidae) Aporrectodea caliginosa (Savigny, 1826), Aporrectodea icterica (Savigny, 1826), and Allolobophora chlorotica (Savigny, 1826); and three anecic or epianecic earthworms Aporrectodea nocturna (Evans, 1946), Lumbricus terrestris (Linnaeus. 1758), and Scherotheca gigas (Dugès, 1828) (Supplementary Table 1). Although the functional trait-based approaches have received a growing interest in soil fauna ecology (Hedde et al., 2022), the choice of relevant traits explaining the effects of soil fauna species on ecosystem carbon (C) fluxes is still challenging, especially when species belong to different clades. Instead, we have used a hybrid classification mixing phylogenetic (woodlice, millipeds, earthworms) and morphoecological groups (anecic or epianecic vs. endogeic). The LowFD communities were composed of three species of the same group, either woodlice, millipede, endogeic or anecic whereas HighFD communities were composed of three species belonging to three different groups. The first fauna addition took place in December 2018, with 4 g fresh weight (FW) per arthropod species and 13 g FW per earthworm species. Subsequent additions of identical biomass were made each year to maintain soil fauna biomass.

The molecular procedure was in accordance with Prada-Salcedo et al. (2022). The DNA extraction and sequencing process were conducted on each soil sample collected from every mesocosm unit. Briefly, the microbial DNA was isolated from 400 mg of each soil sample using a Power Soil DNA Isolation Kit (Qiagen Laboratories Inc., Solana Beach, USA) following the manufacturer’s instructions. After collecting DNA, it was quantified using NanoDrop equipment (Thermo Fisher Scientific, Germany). The V4 region of the 16S rRNA gene was amplified using the primers P5_8N_515F and P5_7N_515F together with P7_2N_806R and P7_1N_806R (Caporaso et al., 2011; Moll et al., 2018). PCR amplification was performed using 7.50 μL of KAPA HiFi HotStart ReadyMix DNA Polymerase (Kapa, Roche or Sigma) and 0.3 μL of each primer with 2 μL of template in 15 μL reaction. The amplification program was initiated at 95°C for 3 min, followed by 30 cycles of 95°C for 50 s, 55°C for 50 s and 72°C for 60 s, with a final extension of 72°C for 7 min. PCR products were purified using AMPure XP beads (Beckman Coulter, Krefeld, Germany). Indexing and sequencing were performed using Nextera XT Illumina Index Kit (Illumina) according to the manufacturer’s instructions. Library quantification was done following Prada-Salcedo et al. (2022) Illumina MiSeq sequencing was performed at the Department of Soil Ecology, UFZ-Helmholtz Centre for Environmental Research in Halle (Saale), Germany.

For subsequent bacterial and archaeal data analyses, the “dadasnake” pipeline was used (Weißbecker et al., 2020), which is based on the DADA2 algorithm (Callahan et al., 2016). This pipeline produces amplicon sequence variants (ASVs). The following parameters were used: minimum read length 170 bp (forward and reverse) for read sequences; maximum expected error after truncation was 0.2. Taxonomic assignment was done using the SILVA database v138.1 (Quast et al., 2013). The statistical analyses were performed with R v.4.2.1. “Phyloseq” (McMurdie and Holmes, 2013) and “vegan” (Oksanen et al., 2014) were the main R-packages used to analyze the retrieved 16S amplicon data. Samples were rarefied at the cutoff of 32,000 reads. The resulting 16,712 ASVs were used to evaluate the microbial communities in this experiment and to predict microbial functional abundances using PICRUSt2 (Douglas et al., 2020) with the “per_sequence_contrib” option. This option predicted the functional gene abundances from 16S rRNA gene sequences and linked them to enzyme groups and pathways from the MetaCyc database.² Linear mixed-effects models (LMM) were used to evaluate differences in relative abundances at the phylum level and microbial diversity indices across treatments and time. The effects of factors on these variables were assessed, and comparisons between treatments and time points were estimated using the “emmeans” R package. The microbial community composition was evaluated with the “vegdist” function using Bray–Curtis dissimilarity and visualized by non-metric multidimensional scaling (NMDS). The impact of soil macrofauna functional diversity on soil microbial communities was tested by permutational multivariate analysis of variance (PERMANOVA). The “indicspecies” (package version 1.7.14) was used to identify microbial indicator taxa associated with drought and macrofauna functional diversity for time points: Prolonged drought and Prolonged Intense drought, because these were the time points to investigate the microbial communities under post-drought conditions. The selected species indicator taxa were tested for differential abundances and represented with heat trees with the “metacoder” (Foster et al., 2017) applying a Wilcoxon Rank Sum test and represented with a log of ratio of median abundances (p value = 0.05).

The results of the predicted functional gene abundances allowed the selection of relevant enzymes associated with abiotic stress for each ASV. This process was determined under two criteria: (1) genes had a relative abundance function between 0.01 and 100, and (2) genes were associated with enzymatic groups related to three types of responses to drought abiotic stress, such as oxidative stress, synthesis of the cell wall, and gene regulation (Table 1). Enzyme abundances were associated with their corresponding taxa at the phylum level and statistically differences analyzed using the Wilcoxon Rank Sum test (p value = 0.05).

Our results indicate that the low dose of gamma irradiation applied after the soil was excavated in 2018 had an effect on microbial species but allowed the majority of the microbial community to persist. Moreover, later in the experiment, after the mesocosms had been established, planted with trees, and supplemented with a litter layer, soil macrofauna, and two inoculations of soil slurry from natural forests, the number of reads and species diversity increased steadily during 2019 and 2020 (Supplementary Figure S1). The continuously changing microbial community over time may also suggest a relatively long establishment phase for microbial communities. We used ASVs from the experimental years 2019 and 2020, to evaluate microbial communities and predict functional gene abundances as a function of the combined treatments of macrofauna functional diversity and drought. The relative abundances at the phylum level throughout the two experimental years for high functional diversity (HighFD) and low functional diversity (LowFD) of macrofauna diversity levels are presented in Figure 2. The results indicate a significant effect of sampling periods, i.e., pre-drought vs. drought periods, on microbial relative abundances (F = 3.29, p < 0.01). We observed a general decrease in relative abundances from 2019 to 2020 for some phyla, such as Gemmatimonadota, Verrucomicrobiota and Planctomycetota, for both macrofauna functional diversity levels. Significant differences across sampling periods, specifically between pre-drought and drought periods, were observed in Proteobacteria and Chloroflexi, which also showed the highest relative abundances. Proteobacteria decreased in abundance during the drought period for both Control and Drought treatments. Conversely, Chloroflexi increased in abundance during drought periods (Figure 2 and Supplementary Table 3).

Overall relative abundances at the phylum level did not differ between high and low macrofauna functional diversity (F = 2.11, p = 0.14). However, the results indicate a significant effect of the drought treatments on the relative abundances (F = 4.99, p = 0.02) and an interaction effect between drought treatments and macrofauna functional diversity levels (F = 3.95, p = 0.04). For example, Gram-positive bacteria, such as Chloroflexi and Actinobacteriota, exhibited higher relative abundances under HighFD conditions during drought (in both years 2019 and 2020) and even during the pre-drought period between the 2019 and 2020 summer droughts compared to the LowFD treatment (Figure 2 and Supplementary Table 4).

Our species indicator analysis for drought periods revealed that from a total of 16,712 ASVs, 130 and 158 ASVs were exclusive to the control treatment for LowFD and HighFD, respectively. In the drought treatments, a total of 173 and 167 ASVs were exclusive for LowFD and HighFD, respectively (across both years 2019 and 2020, Supplementary Figure S2-Venn diagram). Overall, we found differential abundances among taxa during the drought periods of both years (2019 and 2020) between LowFD and HighFD (Figure 3). A consistent trend for both fauna functional diversity levels was that Actinobacteriota was related to drought treatments (2019 and 2020), and at the class level within this phylum, the Thermoleophilia and Rubrobacteria were more abundant under drought than in the control treatments. A similar case was observed for the Planctomycetes, where these taxa were abundant under drought treatments. In an opposite direction, and also consistent between macrofauna functional diversity levels, the Desulfobacterota phylum showed higher abundances in drought than in control treatments (Figure 3). Other class such as Acidimicrobiia, Saccharimonadia, or Myxoccocia showed opposing responses in abundance according to the functional diversity of soil fauna. The differential abundance analysis also revealed differential taxa responding to the control and drought treatments for HighFD and LowFD. In the HighFD, 78 ASVs were detected with significantly higher differential abundances under control conditions, while 87 ASVs exhibited significantly higher differential abundances under drought conditions. For the LowFD samples, 88 ASVs were associated with significantly higher differential abundances in the control treatment, whereas 65 ASVs displayed significantly higher differential abundances under drought conditions (Supplementary Tables 5, 6).

The microbial Shannon diversity was comparable between HighFD and LowFD treatments (Figure 4). Likewise, the Drought treatment had no impact on microbial diversity. However, Shannon diversity generally increased over time for both levels of macrofauna functional diversity (F = 16.44, p < 0.01). The differences over time appeared more pronounced for LowFD than for HighFD treatments. For instance, the Shannon diversity of microbial communities from the HighFD treatment exposed to drought was comparable between 2019 and 2020. Specifically, microbial diversity in the Drought treatment in 2019 was similar to that in the pre-drought period of 2020, and this pre-drought diversity in 2020 was comparable to the Drought treatment diversity in the same year (Figure 4). In contrast, the microbial Shannon diversity of communities from the LowFD treatment differed considerably among measurements over time for the same years, and Drought treatments were less comparable over time (Figure 4 and Supplementary Tables 7, 8).

Our analysis revealed significant differences in microbial community composition. The PERMANOVA results showed a slight but significant effect of macrofauna functional diversity (R² = 0.006, F = 1.51, p = 0.03), and stronger effects were detected among sampling time periods (R² = 0.16, F = 12.65, p < 0.01). Moreover, the results indicate an overall Drought treatment significant effect (R² = 0.008, F = 2.00, p < 0.01). This Drought treatment effect was more pronounced under prolonged and more severe summer drought, showing that HighFD exhibited a weaker difference of the Control Vs. Drought treatment (R² = 0.034, F = 1.64, p < 0.009), for the years 2019 and 2020, compared to the same treatments with LowFD (R² = 0.037, F = 1.77, p < 0.01). NMDS analysis showed similar differences among sampling dates regardless of the macrofauna functional diversity treatment (Figure 5). The first axis primarily separated the data according to time of sampling, with 2019 samples to the left and 2020 samples to the right. Likewise, microbial communities were also mainly separated with ongoing experimental duration, and thus, drought impact along the second axis of the NMDS plot. This resulted in microbial communities shifting from the bottom left to the top right of the NMDS plot in their composition over time, with some differences between Control and Drought.

The estimated abundance of genes coding for microbial enzymes involved in drought response is shown in Figure 6. Overall, the results suggest little impact of prolonged and more severe summer drought on the abundance of drought stress-related enzymes. There were two exceptions, both in the HighFD treatment, with Crenarchaeota and Proteobacteria showing higher enzyme abundances under the Drought condition, compared to the Control (Figure 6). This effect was significant for the histidine kinase and serine–threonine kinase (Supplementary Table 9). Interestingly, there were eleven phyla associated with the expression of drought stress-related enzymes under the LowFD, compared to only nine phyla associated with the HighFD. The Patescibacteria and Thermoplasmatota (Archaea) were phyla present in the LowFD, but absent in the HighFD. These findings suggest that the enzyme abundances were supported by a greater number of taxa under LowFD than under HighFD.