CC-BY 4.0
Introduction
Soil is a non-renewable resource at human life span that carries out several ecosystem services supporting the life of many organisms: it sustains food and fibers production, stores CO2 and water, and helps in disease suppression (Dominati et al., 2010; Aislabie and Deslippe, 2013; Moebius-Clune et al., 2016). Soil microbial communities are major drivers of nutrient cycles that sustain plant growth and productivity. The nitrogen cycle (N-cycle), which regulates the availability of N in soil, is carried out by specific microbial guilds (Dominati et al., 2010; Aislabie and Deslippe, 2013; Whitman et al., 1998; Singh, 2015). Indeed, the N2 in the atmosphere is not available for direct use by plants, while symbiotic nitrogen-fixing bacteria can fix it into plants as ammonia that is then used to build amino acids and proteins (Moebius-Clune et al., 2016). The ammonia oxidizing bacteria and archaea take part to the nitrification converting the NH4+ derived from organic nitrogen into NO2- (Lehtovirta-Morley, 2018), while the newly discovered comammox can transform ammonia directly into NO3- (Koch et al., 2019). The inverse process from NO2- to N2 is performed by the denitrifiers or anammox bacteria (Kuypers et al., 2018). Overall, the N-cycle involves a wide diversity of microbial taxa that are important to ensure its functioning all along seasonal changes of soil temperature and humidity.
Climate change and extreme climatic events can induce shifts in microbial community composition and taxa abundance (Barnard et al., 2013; Jurburg et al., 2017a; Mooshammer et al., 2017; Calderón et al., 2018; Meisner et al., 2018), with effects on nutrient cycles (Philippot et al., 2013) and ecosystem resilience (Jurburg et al., 2017a). Frey et al. (2008) in their long-term experiment in mixed forest concluded that 12 years of soil heating reduced total microbial biomass, and induced a shift in community composition, with a negative impact on substrate use. The N-cycle is also deeply influenced by temperature (Calderón et al., 2018; Mooshammer et al., 2017; Sahrawat, 2008; Szukics et al., 2010), and the perturbation of one of the steps, e.g. nitrification, poses challenges for nitrogen supply in the soil, leading to the disruption of nitrogen availability.
Pesticides are commonly applied in conventional agriculture and in integrated pest management to ensure an adequate production of food and fibres. According to the Food and Agriculture Organization (FAO, 2023), 354’082 tonnes of pesticides were sprayed for agricultural use in the EU in 2021, of which only a part of the applied compounds gets to the target. The rest (30-50%) can reach the soil environment (Rodríguez Eugenio et al., 2018) where it may impact non-target organisms and transfer to other environmental compartments, including water resources. Several studies have already documented that pesticide residues can impact soil microbial community leading to i) the emergence of pesticide-degrading strains able to use pesticide as carbon source for their growth, or ii) the alteration of microbial abundance and diversity and possibly affecting ecosystems stability and their resilience (Karas et al., 2018; Kumar et al., 2019; Puglisi, 2012; Storck et al., 2018; Zabaloy et al., 2010). During the entire crop cycle, often more than one pesticide is applied at the same field. Persistent residues of active ingredients as well as their transformation products accumulate to a cocktail of chemicals to which soil living organisms are exposed (Geissen et al., 2021; Tang and Maggi, 2021). The overall ecotoxicological impact resulting from the exposure of non-target organisms to chemical cocktails is the result of the combination of additive, antagonistic or synergistic effects, that until now have been rarely studied with regard to soil microorganisms. Cedergreen (2014) reviewed the available literature on chemical mixtures, concluding that additive interactions generally dominate whereas evidences for synergistic interactions were limited to cholinesterase inhibitors (insecticide) and cytochrome P450 inhibitors (fungicide). Toxicity endpoints of aquatic organisms dominated this review (Cedergreen 2014), while microbial endpoints were underrepresented. Hence, there is a strong need for a better understanding of the impact of accumulated pesticide residues particularly on the soil microbial community.
Ecosystems are usually exposed to more than one disturbance, that act independently or influence each other, with different effects on the ecological recovery of the ecosystem functions. Very often, a legacy effect of disturbances (Paine et al., 1998) can be observed and might decrease the ability of species to resist and then recover to subsequent perturbations. Jurburg et al. (2017b) demonstrated that a microbial community needed at least 29 days to recover from a heat shock. In addition, the disturbance history (heat, and heat or cold shock) affects the microbial community response to the subsequent perturbations (Jurburg et al., 2017a). Extensive works have been done to assess the effect of factors mimicking extreme climate events (Bardgett and Caruso, 2020; Calderón et al., 2018; Meisner et al., 2018; Mooshammer et al., 2017; Philippot et al., 2021), but up to now, only a few studies have considered the potential combined effect of environmental disturbances and pesticides on soil microbial communities (Franco-Andreu et al., 2016; Pesaro et al., 2004).
Within this context, the objectives of the present work are, therefore, to compare effects of global change-related environmental disturbances and impacts of pesticidal active ingredients on soil microbial community structure and functioning and to evaluate their compounded effects. We hypothesize that, as compared to undisturbed soil i) extreme climatic events, such as elevated temperature and heavy rainfall will change microbial community’s structure and functioning; ii) subsequent exposure to different pesticidal active ingredients applied alone or in combination may lead to compounded effects on the structure and functioning of some key microbial community members, namely guilds involved in N-cycle.
Material and methods
Soil sampling and microcosms set up
Soil from the surface layer (0-20 cm depth) was sampled from Sayens field site located in the vicinity of the INRAE experimental farm in Bretenière, France (47° 14' 03.1" N, 5° 06' 37.9" E) during August 2021. This soil is silty clay loam, characterized by 34.6% clay, 59.0% silt, and 6.3% sand, with an organic carbon content of 27.3 g.kg-1 and a nitrogen content of 2,7 g.kg-1. The water holding capacity (WHC) was 64. The soil was homogenized and sieved to 4 mm and stored at 4°C until further use. The soil humidity after sieving was 26% which corresponds to ~ 40% of the WHC. The fresh soil was divided into 250 microcosms, each containing 63 g fresh soil (corresponding to 50 g dry soil), and incubated at room temperature for 6 days prior to exposure, or not, to environmental disturbance.
Environmental disturbances and active ingredients application
We studied two different environmental disturbance factors: elevated temperature (heat) and heavy rainfall (high humidity). Respectively, 80 soil microcosms were subjected to the heat disturbance, 80 to high humidity and 80 received no environmental disturbance (control). Heat and high humidity disturbances were applied over a window of 7 days. The heat disturbance consisted in two periods of 30 h at 42°C, separated by a period of 40 h at room temperature. At the end of the disturbance window, the humidity was adjusted to 22%. The high humidity disturbance consisted in watering soil microcosms to reach 35% of humidity, followed by an evaporation period of 7 days. Three days after the end of the disturbance window, i.e. at day 19 of the experiment, five microcosms per each disturbance condition tested (heat, high humidity and control) were treated with one of the three pesticides (technical material of active ingredients, a.i.) or with the mixture of the three a.i at 1x or 10x the agronomical dose. The pesticides were: 3,6-dichloropyridine-2-carboxylic acid (i.e. clopyralid, pesticide analytical standard, CAS 1702-17-6, supplier: Sigma-Aldrich, purity ≥ 98.0%, agronomical dose: 125 g/ha corresponding to 4.5 µg/50 g dry soil), [cyano-(3-phenoxyphenyl)methyl] 3-(2,2-dichloroethenyl)-2,2-dimethylcyclopropane-1-carboxylate (i.e. cypermethrin, pesticide analytical standard, CAS 52315-07-8, supplier: Sigma-Aldrich, purity ≥ 90.0%, agronomical dose: 100 g/ha corresponding to 3.4 µg/50 g dry soil), and methyl N-[2-[[1-(4-chlorophenyl)pyrazol-3-yl]oxymethyl]phenyl]-N-methoxycarbamate (i.e. pyraclostrobin, pesticide analytical standard, CAS 175013-18-0, supplier: Sigma-Aldrich, purity ≥ 98.0%, agronomical dose: 200 g /ha corresponding to 7.1 µg/50 g dry soil). These three a.i are all constituents of commercial formulated products commonly used in agriculture, and can be all used in corn cropping during a growth cycle. Clopyralid was dissolved in water, and cypermethrin and pyraclostrobin were dissolved in acetone, and equal quantities of water (1 mL) and acetone (10 µL) a.i. solutions were applied to all microcosms with the respective treatment (single a.i or a.i. mixture). The solutions were applied onto the soil surface, without mixing the soil after application. In a similar manner, water and acetone were applied to 10 untreated control microcosms. Soil moisture was adjusted and kept to 26% in all microcosms. All the treated microcosms and the untreated controls were done in five independent replicates, for a total of 250 microcosms (Figure 1).
From each treatment combination and from the untreated controls, five microcosms were destructively sampled at day 29 (T30: 10 days after a.i application), and the remaining five at day 61 (T60: 42 days after a.i application). The soil samples were processed fresh or stored at -20°C for further analyses.
Figure 1 - Graphic representation of the experimental timeline with a description of the experimental design
Nitrogen pools quantification
The nitrogen species NO3- and NH4+ were extracted from 10 g (fresh weight) of the sampled soil using 50 ml of K2SO4 0.5 M. The soil suspension was shaken for one hour, allowed to settle for two-three hours at room temperature, and then decanted. The supernatant was filtered and kept frozen until quantification. The quantification was performed by colorimetry according to ISO 14256-2:2005. Two blanks per series were included. The results of the NO3- and NH4+ quantifications were expressed as mg of N/kg of dry soil.
DNA extraction
DNA was extracted from 250 mg of dry weight soil using the DNeasy PowerSoil DNA extraction Kit (Qiagen, FR) following the kit’s instructions. The DNA was quantified with the Quant-iT™ dsDNA Assay Kit (Thermo Fischer Scientific, FR), and stored at -20°C until further use.
Assessment of microbial community composition and diversity
Bacterial community composition and diversity were monitored via sequencing the 16S rDNA V3-V4 hypervariable region. We use a two-step PCR amplification method (Berry et al., 2011). In the first step, 25 amplification cycles were performed using the fusion primers U341F (5’-CCTACGGGRSGCAGCAG-3’) and 805R (5’-GACTACCAGGGTATCTAAT-3’) as described in Takahashi et al. (2014), with overhang adapters (forward: 5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3’, adapter: 5’-GTCTCGTG GGCTCGGAGATGTGTATAAGAGACAG-3’). The thermal cycling conditions of the first PCR were: 3 min at 94°C, 25*(30 sec at 94 °C, 30 sec at 55 °C, and 30 sec at 72 °C), and 10 min at 72°C. Duplicate of the first PCR were pooled and used in a second PCR, were the amplification linked a unique combination of multiplex primer pair to the overhang adapters for each sample. The thermal cycling conditions of the second PCR were: 3 min at 98°C, 8*(15 sec at 98°C, 30 sec at 55°C, 30 sec at 72°C), and 10 min at 72°C.
The amplicons from the second PCR were pooled and visualized in a 2% agarose gel for size and amplification check. They were then pooled and purified with the SequalPrep Normalization Plate kit (Invitrogen, Frederick, MD, USA). Sequencing was performed by GenoScreen (Lille, FR) on MiSeq (Illumina 2*250 bp) using the MiSeq reagent kit v2 (500 cycles).
Quantification of microbial communities
The total abundances of bacteria (16S), fungi (ITS), and protists (18S), as well as the N-cycle microbial guilds were estimated through real-time quantitative PCR (qPCR). The 16S rDNA and ITS2 primers described by Muyzer et al. (1993) and White et al. (1990) were used to estimate the total 16S and ITS communities abundances, respectively. The 18S community abundance was estimated from the 18S rDNA using the primers EK-565F (5′-GCAGTTAAAAAGCTCGTAGT-3′) and 18S-EUK-1134-R–UnonMet (5′-TTTAAGTTTCAGCCTTGCG-3′) with cycling condition of 3 min at 95°C, 35*(10 sec at 95 °C, 30 sec at 58 °C, 30 sec at 72 °C and 20 sec at 80°C). The ammonia oxidizing archaea (AOA) and bacteria (AOB) were quantified through gene amoA (Bru et al., 2011), comammox A (ComaA) and B (ComaB) were targeted with the clade A and B amoA genes (Pjevac et al., 2017), respectively.
The quantifications were carried out in a ViiA7TM (Life Technologies, Carlsbad, CA, USA) with a reaction volume of 15 µL containing 7.5 µL of Takyon Master Mix (Eurogentec, Liège, BE), 1 µM of each primer, 250 ng of T4 gene 32 (MP Biomedicals, Santa Ana, CA, USA), 3 ng of DNA. Two independent replicates were used for each qPCR assay. Prior to qPCR analysis, an inhibition test was performed by mixing the soil DNA extracts with a control plasmid DNA (pGEM-T Easy Vector, Promega, Madison, WI, USA) or water. No inhibition was detected. The results from the quantification analysis were expressed as number of gene copies/ng of DNA.
Bioinformatics analysis
Sequencing data were analyzed using an in-house developed Jupyter Notebooks (Kluyver et al., 2016) piping together different bioinformatics tools. Briefly, for 16S, R1 and R2 sequences were assembled using PEAR (Zhang et al., 2014) with default settings. Further quality checks were conducted using the QIIME pipeline (Caporaso et al., 2010) and short sequences were removed (< 400 bp for 16S). Reference based and de novo chimera detection, as well as clustering in OTUs were performed using VSEARCH (Rognes et al., 2016) and the SILVA v132 reference database. The identity thresholds were set at 94% for 16S rDNA data based on replicate sequencings of a bacterial mock community containing 40 bacterial species. Representative sequences for each OTU were aligned using infernal (Nawrocki and Eddy, 2013) and a 16S phylogenetic tree was constructed using FastTree (Price et al., 2009). Taxonomy was assigned using BLAST (Altschul et al., 1990) and the SILVA reference database v132 (Quast et al., 2013).
Diversity metrics, that is, Faith’s Phylogenetic Diversity (Faith, 1992), richness (observed species) and evenness (Simpson’s reciprocal index), describing the structure of microbial communities were calculated based on rarefied OTU tables (8300 sequences per sample). Weighted UniFrac distance matrices (Lozupone and Knight, 2005) were also computed to detect global variations in the composition of microbial communities.
Statistical analysis of the microbial endpoints and of the nitrogen pools
The statistical analysis was carried out using RStudio statistical software version 4.2.1. The data were checked for homogeneity of variance with Levene’s test (p < 0.05). Data were log-transformed or cos-transformed when necessary and outliers were removed when detected (a maximum of 2 values were removed for any given variable).
The values from 16S, ITS, 18S, AOA, AOB, ComaA, ComaB abundances, NO3- and NH4+ concentration, and α-diversity indices which are observed species (OS), Faith’s Phylogenetic Diversity (PD) and Simpson’s reciprocal index (SR) were analyzed with linear ANOVAs (n=5 per treatment at each sampling date) to determine any significant variance in the measured variables that can be attributed to the time of sampling, environmental disturbance (heat, high humidity, control), pesticide treatment (three different a.i and their mixture) and dose (1x, 10x) combinations on microbial parameters. Post-hoc analyses were performed with Tukey’s honestly significant difference test (p value < 0.05) using the TukeyHSD() function from the stats package in RStudio.
The data from β-diversity weighted UniFrac were analyzed with PermANOVA (n=5) using the function adonis from the R package vegan to detect significant differences between time, environmental disturbance, pesticide and dose combinations on bacterial community composition. Post-hoc analyses with pairwise comparisons were performed when significant differences were observed (p value < 0.05).
To estimate which OTUs significantly differ between our experimental conditions, we used a generalized linear mixed model. This model accounts for the non-normal distribution of abundance data that typically follow a Poisson distribution, and contains both fixed effects (experimental conditions: time of sampling, environmental disturbances, pesticide treatments and doses) and random effects (sampling effects). Multiple pairwise comparisons were performed with the emmeans() function of the emmeans package with the Bonferroni correction for p-value adjustment.
Results
We assessed here the impact of environmental disturbances (heat or high humidity) and a.i exposure (three different a.i and their mixture) on several endpoints: i) abundance of 16S, ITS and 18S communities, ii) structure and composition of the soil bacterial community (α- and β-diversity) and iii) N-cycle functional guilds, by targeting the nitrification step by different means (NH4+ and NO3- pools, AOA, AOB, ComaA and ComaB abundances).
Effect of environmental disturbances on microbial endpoints
We performed a Principal Component Analysis (PCA) to decipher the contribution of each measured endpoint to the observed response of soil microcosms to environmental disturbances as compared to the control at each time of sampling (T30 and T60, Figure 2 panel A and panel B, respectively). Most of the variance was explained by the two first axes (respectively ~75 % and ~65 % at T30 and T60) and separated the heat-disturbed microcosms from the high humidity and the control ones (Figure 2). At both timepoints, this was mostly explained by α-diversity indices and NO3- pool (Table S1).
Figure 2 - Principal Component Analysis (PCA) of the microbial endpoints (abundance of microbial communities 16S, 18S and ITS; bacterial diversity indices OS, PD and SR; N-cycle functional guilds: abundance of AOA, AOB, ComaA and comaB and abundance of nitrate (NO3-) and ammonium (NH4+) pools) measured at T30 (panel A) and T60 (panel B). PC1 and PC2 axes represent the two main principal components capturing the maximum variability in the dataset. Arrows indicate the direction and magnitude of each endpoint’s contribution to the first and second principal components. Each soil microcosm is depicted by a symbol whose shape, color and fill respectively indicate the disturbance, the active ingredient (a.i) and the a.i. doses that were applied.
Focusing on disturbed and control microcosms irrespectively of the applied a.i, we observed significant differences between heat-disturbed microcosms and the others (Table 1). The abundances of 16S, ITS and 18S communities, as well as the diversity of the bacterial community (OS, PD and SR) were significantly lower in heat-disturbed microcosms (ANOVAs followed by Tukey HSD tests; p < 0.05) as compared to the control. When focusing on N-cycle, we also observed a significant decrease in the abundance of the ComaA and ComaB guilds and a significant increase in NH4+ (only at T30) and NO3- pools, as well as in the abundance of AOB in heat-disturbed microcosms. For most of the endpoints, the high humidity disturbance did not induce changes as compared to control microcosms, except for the diversity level of the bacterial community at T30, and the abundance of the protist community (18S).
We also performed a Principal Coordinates Analysis (PCoA) of the weighted UniFrac distance matrix to examine bacterial compositional changes in response to environmental disturbances applied to soil microcosms (Figure 3). In line with results presented above, heat-disturbed microcosms were significantly different from the other two disturbance treatments (high humidity and control) at both timepoints (PermANOVAs; p < 0.05). In addition, the high-humidity-disturbed microcosms were significantly slightly different from control microcosms at T30. Overall, among the 336 dominant OTUs (defined as those having abundance > 0.1 %), ~80% showed at both T30 and T60, significant differences in their abundance between heat-disturbed and control microcosms, while less than 5% displayed significant differences between high-humidity-disturbed and control ones (Figure S1).
Table 1 - Microbial endpoints (abundance of microbial communities 16S, 18S and ITS; bacterial diversity indices OS, PD and SR; N-cycle functional guilds: abundance of AOA, AOB, ComaA and comaB and abundance of nitrate (NO3-) and ammonium (NH4+) pools) measured at T30 and T60 on disturbed (HEAT, HIGH HUMIDITY) and control (CTRL) soil microcosms. For each endpoint, the mean ± the standard error across all pesticide treated microcosms (n= 40 for each disturbance treatment and control) measured at each time point and for each soil disturbance is indicated. For a given time, an asterisk (*) indicates that the endpoint is significantly different from the one measured in the undisturbed soil (ANOVA followed by a Tukey HSD test; p < 0.05). SR: Simpson Reciprocal; PD: PD Whole Tree; OS: Observed Species.
T30 | T60 | ||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
CTRL | HIGH HUMIDITY | HEAT | CTRL | HIGH HUMIDITY | HEAT | ||||||||||
Microbial abundance | ITS (x 103 per ng) | 10.9 | ± 0.5 | 11.9 | ± 0.7 | 5.2 | ± 0.3* | 13.5 | ± 0.8 | 15.5 | ± 0.8 | 7.2 | ± 0.4* | ||
18S (x 103 per ng) | 52.6 | ± 1.9 | 56.7 | ± 3.2* | 23.3 | ± 0.9* | 65.8 | ± 3.5 | 76.6 | ± 3.5* | 31.7 | ± 1.3* | |||
16S (x 103 per ng) | 128 | ± 4 | 140 | ± 7 | 110 | ± 4* | 134 | ± 6 | 145 | ± 5 | 123 | ± 4* | |||
Bacterial α-diversity | SR | 240 | ± 3 | 247 | ± 3 | 140 | ± 2* | 262 | ± 3 | 265 | ± 2 | 179 | ± 4* | ||
PD | 139 | ± 1 | 142 | ± 1* | 122 | ± 1* | 143 | ± 1 | 142 | ± 0 | 128 | ± 1* | |||
OS | 1799 | ± 8 | 1832 | ± 7* | 1542 | ± 9* | 1840 | ± 8 | 1836 | ± 6 | 1626 | ± 9* | |||
N-cycle functional guild | ComaA (per 104 16S) | 71.3 | ± 1.6 | 76.2 | ± 3.1 | 52.5 | ± 1.6* | 80.0 | ± 1.7 | 83.3 | ± 2.6 | 67.0 | ± 1.6* | ||
ComaB (per 104 16S) | 14.3 | ± 0.4 | 15.4 | ± 0.5 | 10.5 | ± 0.3* | 21.5 | ± 0.5 | 22.5 | ± 0.8 | 18.6 | ± 0.5* | |||
AOA (per 104 16S) | 276 | ± 12 | 301 | ± 12 | 270 | ± 12 | 325 | ± 10 | 326 | ± 12 | 358 | ± 12 | |||
AOB (per 104 16S) | 27 | ± 1.0 | 28 | ± 0.9 | 57 | ± 3.0* | 30 | ± 1.0 | 33 | ± 1.3 | 45 | ± 2.8* | |||
NO3- | 1.8 | ± 0.3 | 1.3 | ± 0.2 | 46.2 | ± 0.6* | 1.6 | ± 0.2 | 1.6 | ± 0.2 | 46.8 | ± 1.1* | |||
NH4+ | 1.1 | ± 0.1 | 0.9 | ± 0.1 | 2.1 | ± 0.3* | 1.4 | ± 0.1 | 1.7 | ± 0.2 | 1.2 | ± 0.1 | |||
Figure 3 - Principal Coordinate Analysis (PCoA) of the bacterial community structure in the different soil microcosms based on the weighted UniFrac distance metric. The two main principal coordinates PCo1 and PCo2 captured respectively 54.3 % and 7.5 % of the observed variability between bacterial communities. For T30 (panel A) and T60 (panel B), each sample is represented by a symbol whose shape, color and fill respectively indicate the disturbance, the active ingredient and the a.i. dose that were applied.
Evaluating the impact of active ingredients on microbial endpoints
We then focused on the results obtained for the disturbance control microcosms exposed to a.i applied alone (clopyralid, i.e. CLO; cypermethrin i.e. CYP; pyraclostrobin, i.e. PYR) or as mixture (i.e. MIX) to evaluate their ecotoxicological impacts on measured endpoints. When considering all the endpoints together in a PCA (Figure 4, panels A and B), we did not observe a strong clustering of the samples, neither according to the a.i or mixture of the three a.i, nor to the dose of applications (1x or 10x) at both timepoints (T30 and T60). However, when considering the individual endpoints, we detected a significant effect of a.i treatments on NH4+ and NO3- pools, as well as on the abundance of AOA, ComaA and ComaB (even though no significant differences were detected for those groups in pairwise comparisons with the untreated control microcosms) and to a smaller extent on the abundance of 16S, 18S and ITS (ANOVAs, Pesticide_Dose and/or Pesticide_Dose-by-Time effects p < 0.05). At T30, the NO3- pool was the most sensitive endpoint, with significantly higher levels detected in the CLO1x, CLO10x, CYP1x, PYR1x and MIX1x treated soil microcosms as compared to the untreated controls (Table 2). The NH4+ pool was significantly higher only in the CLO1x and MIX1x treated microcosms as compared to the untreated control. At T60, the 16S abundance was the sole endpoint that was significantly decreased in response to exposure to CYP1x and PYR1x.
Table 2 - Microbial endpoints (abundance of microbial communities 16S, 18S and ITS; bacterial diversity indices OS, PD and SR; N-cycle functional guilds: abundance of AOA, AOB, ComaA and comaB and abundance of nitrate (NO3-) and ammonium (NH4+) pools) measured at T30 and T60 on the undisturbed treated and untreated control (u.c.) soil microcosms. For each endpoint, the mean ± the standard error across all control treated microcosms and u.c. (n=5 for each active ingredient-dose and u.c.) measured at each time point and for each active ingredient is indicated. For a given time, an asterisk (*) indicates that the endpoint is significantly different from the one measured in the u.c. soil (ANOVA followed by a Tukey HSD test; p < 0.05). SR: Simpson Reciprocal; PD: PD Whole Tree; OS: Observed Species.
| Part 1 | T30 | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| U.C. | CLO1x | CLO10x | CYP1x | CYP10x | PYR1x | PYR10x | MIX1x | MIX10x | ||
| Microbial abundance | ITS (x 103 per ng) | 10.0 ± 0.5 | 11.9 ± 1.3 | 10.5 ± 1.2 | 13.3 ± 1.8 | 11.0 ± 1.3 | 11.8 ± 2.1 | 8.7 ± 0.6 | 9.6 ± 0.7 | 10.4 ± 1.1 |
| 18S (x 103 per ng) | 46.9 ± 3.1 | 55.3 ± 3.6 | 56.3 ± 6.3 | 61.9 ± 4.8 | 49.7 ± 4.0 | 61.3 ± 8.5 | 44.4 ± 1.8 | 42.2 ± 3.0 | 49.5 ± 5.0 | |
| 16S (x 103 per ng) | 102 ± 6 | 132 ± 7 | 139 ± 12 | 137 ± 9 | 131 ± 13 | 142 ± 21 | 115 ± 7 | 113 ± 4 | 115 ± 5 | |
| Bacterial α-diversity | SR | 244 ± 6 | 244 ± 12 | 229 ± 6 | 232 ± 8 | 251 ± 6 | 240 ± 6 | 234 ± 12 | 240 ± 6 | 254 ± 8 |
| PD | 142 ±1 | 139 ± 2 | 138 ± 1 | 138 ± 1 | 140 ± 3 | 139 ± 0 | 137 ± 2 | 138 ± 2 | 144 ± 2 | |
| OS | 1827 ±12 | 1792 ± 22 | 1777 ± 7 | 1780 ± 11 | 1790 ± 35 | 1793 ± 12 | 1771 ± 23 | 1787 ± 26 | 1872 ± 20 | |
| N-cycle functional guild | ComaA (per 104 16S) | 73.7 ± 5.6 | 64.9 ± 2.0 | 62.4 ± 2.7 | 65.0 ± 2.7 | 80.4 ± 4.1 | 62.7 ± 2.1 | 71.0 ± 2.5 | 81.4 ± 2.9 | 83.1 ± 4.3 |
| ComaB (per 104 16S) | 20.1 ± 1.6 | 12.5 ± 0.3 | 13.6 ± 0.8 | 13.5 ± 0.7 | 15.5 ± 1.5 | 12.1 ± 0.5 | 14.8 ± 0.5 | 15.2 ± 0.6 | 16.9 ± 1.1 | |
| AOA (per 104 16S) | 265 ± 21 | 233 ± 26 | 235 ± 7 | 274 ± 41 | 287 ± 31 | 249 ± 34 | 284 ± 26 | 313 ± 37 | 330 ± 49 | |
| AOB (per 104 16S) | 26.7 ± 0.9 | 30.9 ± 3.9 | 24.9 ± 2 | 27.9 ± 2.3 | 31.8 ± 3.0 | 29.2 ± 3.0 | 22.0 ± 0.6 | 25.9 ± 3.4 | 27.2 ± 2.9 | |
| NO3- | 0.1 ± 0.1 | 1.4 ± 0.5* | 3.1 ± 0.7* | 3.8 ± 1.0* | 0.1 ± 0.0 | 4.4 ± 1.0* | 0.1 ± 0.0 | 1.3 ± 0.6* | 0.6 ± 0.4 | |
| NH4+ | 0.8 ± 0.1 | 1.7 ± 0.3* | 0.8 ± 0.1 | 0.9 ± 0.2 | 0.6 ± 0.1 | 1.4 ± 0.2 | 0.8 ± 0.1 | 1.9 ± 0.4* | 1.1 ± 0.1 | |
| Part 2 | T60 | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| U.C. | CLO1x | CLO10x | CYP1x | CYP10x | PYR1x | PYR10x | MIX1x | MIX10x | ||
| Microbial abundance | ITS (x 103 per ng) | 18.7 ± 2.1 | 11.4 ± 2.6 | 15.1 ± 1.2 | 11.4 ± 2.7 | 13.8 ± 1.9 | 10.9 ± 1.4 | 15.9 ± 2.4 | 13.1 ± 2.0 | 18.5 ± 2.5 |
| 18S (x 103 per ng) | 79.4 ± 6.2 | 60.7 ± 12.3 | 70.2 ± 4.2 | 54.0 ± 12.0 | 72.2 ± 9.2 | 54.9 ± 7.1 | 83.4 ± 10.5 | 67.4 ± 10.8 | 80.0 ± 10.1 | |
| 16S (x 103 per ng) | 200 ± 16 | 133 ± 25 | 129 ± 7 | 107 ± 4* | 128 ± 18 | 119 ± 10* | 148 ± 13 | 136 ± 20 | 163 ± 20 | |
| Bacterial α-diversity | SR | 259 ± 2 | 264 ± 9 | 264 ± 6 | 266 ± 7 | 265 ± 10 | 258 ± 8 | 257 ± 9 | 268 ± 8 | 251 ± 5 |
| PD | 142 ± 1 | 142 ± 1 | 141 ± 2 | 144 ± 1 | 143 ± 2 | 142 ± 2 | 143 ± 2 | 144 ± 2 | 142 ± 2 | |
| OS | 1852 ± 7 | 1825 ± 20 | 1814 ± 24 | 1868 ± 16 | 1839 ± 21 | 1833 ± 21 | 1848 ± 30 | 1861 ± 25 | 1830 ± 23 | |
| N-cycle functional guild | ComaA (per 104 16S) | 77.3 ± 3.3 | 69.5 ± 1.6 | 86.4 ± 4.8 | 79.0 ± 3.3 | 96.3 ± 3.5 | 80.7 ± 6.5 | 85.0 ± 12.3 | 73.8 ± 3.4 | 76.7 ± 2.3 |
| ComaB (per 104 16S) | 20.5 ± 1 | 19.2 ± 1.5 | 22.7 ± 1.2 | 21.4 ± 1.0 | 24.5 ± 1.7 | 22.4 ± 1.5 | 26.8 ± 3.9 | 19.5 ± 1.6 | 21.6 ± 1.2 | |
| AOA (per 104 16S) | 357 ± 35 | 253 ± 15 | 329 ± 33 | 310 ± 29 | 395 ± 23 | 320 ± 23 | 344 ± 23 | 318 ± 18 | 315 ± 26 | |
| AOB (per 104 16S) | 33.7 ± 3.2 | 30.7 ± 4 | 29.7 ± 1.4 | 26.7 ± 2.1 | 31.1 ± 3.7 | 27.1 ± 1.4 | 26.5 ± 1.1 | 30.9 ± 3.9 | 27.7 ± 1.9 | |
| NO3- | 1.1 ± 0.1 | 1.8 ± 0.8 | 1.4 ± 0.2 | 3.4 ± 0.2 | 2.6 ± 0.8 | 0.8 ± 0.2 | 0.7 ± 0.1 | 1.4 ± 0.2 | 0.7 ± 0.1 | |
| NH4+ | 0.8 ± 0.1 | 2.4 ± 0.5 | 0.8 ± 0.1 | 1.8 ± 0.1 | 1.3 ± 0.2 | 1.5 ± 0.2 | 1.4 ± 0.1 | 1.4 ± 0.2 | 1.5 ± 0.2 | |
Figure 4 - Principal Component Analysis (PCA) performed on the microbial endpoints (abundance of microbial communities 16S, 18S and ITS; bacterial diversity indices OS, PD and SR; N-cycle functional guilds: abundance of AOA, AOB, ComaA and comaB and abundance of nitrate (NO3-) and ammonium (NH4+) pools) measured in the undisturbed microcosms at T30 (panel A) and T60 (panel B). The PC1 and PC2 axes represent the two main principal components capturing the maximum variability in the dataset. Arrows indicate the direction and magnitude of each endpoint’s contribution to the first and second principal components. Each soil microcosm is depicted by a circle whose color denotes the active ingredient that it received. Empty circles represent samples having received the a.i. at the agronomical dose while the full circle represent the one who received 10 times the agronomical dose.
Figure 5. Principal Coordinate Analysis (PCoA) of the bacterial community structure in the undisturbed microcosms based on weighted UniFrac distance metric. The two main principal coordinates PCo1 and PCo2 captured respectively 27.2% and 11.5% of the variability observed between the bacterial communities. For T30 (panel A) and T60 (panel B), each sample is represented by a symbol whose color and fill respectively indicate the active ingredient and its dose applied in soil microcosms.
Table 3 - Analysis of variance (ANOVA) results assessing the impact of time of sampling, the disturbance, and the applied pesticide (nature/dose) on various microbial endpoints (16S: bacterial abundance; OS: Observed Species; AOB: Ammonia Oxydizing Bacteria abundance) and nitrate (NO3-) and ammonium pools (NH4+). Only the endpoints for which a significant 2-way ‘Disturbance-by-Pesticide_Dose’ or 3-way ‘Disturbance-by-Pesticide_Dose-by-TimeOfSampling’ interactions are presented (p < 0.05).
| 16S | OS | AOB | NO3- | NH4+ | ||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Df | Mean Sq | F value | Pr(>F) | Mean Sq x10 3 | F value | Pr(>F) | Mean Sq | F value | Pr(>F) | Mean Sq | F value | Pr(>F) | Mean Sq | F value | Pr(>F) | |||||
| TimeOfSampling | 1 | 0.04 | 4.3 | 0.040 | 114 | 60 | < 0.001 | 0.00 | 0 | 0.985 | 1.90 | 4.70 | 0.031 | 0.537 | 13.4 | < 0.001 | ||||
| Disturbance | 2 | 0.14 | 13.5 | < 0.001 | 1569 | 821 | < 0.001 | 1.31 | 114 | < 0.001 | 1.35 | 3.32 | 0.038 | 0.130 | 3.2 | 0.041 | ||||
| Pesticide_Dose | 7 | 0.02 | 2.2 | 0.037 | 10 | 5 | < 0.001 | 0.01 | 1 | 0.554 | 0.55 | 1.36 | 0.224 | 0.151 | 3.8 | 0.001 | ||||
| Disturbance-by-TimeOfSampling | 2 | 0.01 | 0.8 | 0.466 | 32 | 17 | < 0.001 | 0.15 | 13 | < 0.001 | 0.06 | 0.16 | 0.852 | 0.964 | 24.1 | < 0.001 | ||||
| Pesticide_Dose-by-TimeOfSampling | 7 | 0.01 | 1.1 | 0.373 | 3 | 2 | 0.089 | 0.02 | 2 | 0.053 | 0.81 | 1.99 | 0.059 | 0.145 | 3.6 | 0.001 | ||||
| Disturbance-by-Pesticide_Dose | 14 | 0.02 | 1.9 | 0.030 | 3 | 1 | 0.192 | 0.01 | 1 | 0.263 | 1.37 | 3.40 | < 0.001 | 0.133 | 3.3 | < 0.001 | ||||
| Disturbance-by-Pesticide_Dose-by-TimeOfSampling | 14 | 0.01 | 0.7 | 0.781 | 4 | 2 | 0.013 | 0.03 | 2 | 0.007 | 0.68 | 1.69 | 0.061 | 0.064 | 1.6 | 0.084 | ||||
PCoA on the Weighted UniFrac distance matrix (Figure 5) showed that in the absence of heat or humidity disturbance the bacterial community structure was not changed by a.i exposure as compared to the untreated control. We noticed that the structure of the bacterial communities of the microcosms exposed to CLO10x were significantly different from those exposed to CYP10x (Pairwise PermANOVAs, p < 0.05). When looking at the dominant OTUs, at both sampling times, we found only 12 OTUs that showed significant differences in relative abundance in a.i treated microcosms versus the untreated control ones (4 Acidobacterias, 3 Actinobacterias, 2 Verrucomicrobias, 2 Proteobacteria and 1 Chloroflexi). Among those, only four were consistently and significantly impacted at both timepoints: they were observed only in CLO-treated microcosms (three in CLO10x and one in CLO1x) and three of them are related to uncultured members of the Acidobacteria phylum.
Estimating compounded effects of environmental disturbances and pesticides inputs on microbial endpoints
A compounded effect in our ANOVA model would manifest with a significant two-way ‘Disturbance-by-Pesticide_Dose’ or three-way ‘Disturbance-by-Pesticide_Dose-by-TimeOfSampling’ effect, indicating that the effect of a given ‘Pesticide_Dose’ at a given ‘TimeOfSampling’ is conditioned by the ‘Disturbance’. The AOB abundance and the OS of the bacterial community were the only two endpoints showing a significant three-way interaction (Table 3). However, none of the pairwise comparisons were significant for OS, while the only significant, but slight, difference observed for AOB was between T30-Heat-PYR10x (38 ± 8 copies per 104 copies 16S rDNA) and T30-Heat-CLO10x (75 ± 2 copies per 104 copies 16S rDNA) samples, indicating that the heat disturbance conditioned the relative effects of PYR and CLO at the 10x dose. The two-way interaction was only significant for the abundance of 16S, and NH4+ and NO3- pools. More specifically, CYP1x-treated microcosms showed significantly higher levels of NO3- compared to PYR1x, PYR10x and CYP10x-treated microcosms, but only in the treated controls, while PYR10x-treated microcosms displayed significantly higher levels of NH4+ than PYR1x, CYP1x, CLO1x and MIX10x-treated ones, but only for the heat-disturbed microcosms. Overall, compounded effects are relatively sparse, idiosyncratic and of small scale which is in line with the relatively small ecotoxicological impact detected for the chosen pesticides in undisturbed environmental conditions.
Discussion
The present study aimed to evaluate the effects of elevated temperature (heat) and heavy rainfall (high humidity), two global change related environmental disturbances, and of three commonly used a.i., CLO (herbicide), CYP (insecticide), PYR (fungicide) and their mixture on soil microbial community abundance, structure, composition and functioning. We used qPCR techniques to detect any effect on the abundance of some important microbial guilds, and Illumina next generation sequencing of the 16S rDNA amplicons to study the effects on the bacterial community composition.
Most of the measured microbial endpoints were affected by the heat stress, with significant effects at T30 (short term). The α-diversity indices and important N-cycle-related microbial guilds were the endpoints mostly impacted by the heat disturbance. We found decreasing α-diversity indices and accumulation of NO3-, which is in line with former studies in which the disruption of the N-cycle was described as a consequence of a lower functional redundancy due to community shifts (Calderón et al., 2018). In another research, the heat-disturbed microbial communities never recovered, and were always different from the control, independently from the soil disturbance history (Jurburg et al., 2017a). This response is in line with our observations, where the effect of heat stress persisted all along the experiment. We also measured a higher AOB abundance in the heat disturbed samples compared to the control, which suggests a higher tolerance to heat of the AOB community compared to other members of N microbial guilds. This finding is supported by other studies, in which AOB were more responsive to soil warming than AOA (Szukics et al., 2010; Wang et al., 2021). We observed that ComaA and ComaB abundances were negatively impacted by the heat stress: not much is known about the effect of heat on these two microbial guilds, but a recent research suggests that commamox are favoured at low temperature in wastewater treatment plants (Gonzalez-Martinez et al., 2016). We can therefore conclude that the N-cycle was impacted by heat stress at several levels, and exhibited a limited recovery, i.e. low resilience, at the time scale of the experiment.
The effect of pesticides on soil microbial communities has been studied for decades (Bollen, 1961). The microbial response to pesticide application is variable, and depends on many factors such as the pesticidal mode of action, the experimental set up and the studied endpoints (Jacobsen and Hjelmsø, 2014; Puglisi, 2012). Indeed, the few responses observed in our experiment are not consistent across the various endpoints analysed. Most are related to the N-cycle and exhibit a transient effect, i.e. recovery in the longer term. Some functional guilds, e.g. AOA and AOB, have been described to be sensitive to pesticide application (Karas et al., 2018). To our best knowledge, no studies have explored the influence of pesticides on comammox. Even though we could not detect a strong effect, our findings indicate that the application of technical a.i might affect these N-cycle microbial guilds. The exposure of soil to the a.i, and particularly to the herbicide CLO, led to a significant increase of the NO3- and NH4+ pools which is supported by previous observations reporting the disruption of the N-cycle in response to pesticide application (Brochado et al., 2023; Damin and Trivelin, 2011; Hernández et al., 2011). Interestingly, one could observe that exposure to a.i had significant effects on the above mentioned endpoints already at agronomical dose, contrary to what is reported in the existing literature where effects are mostly detected at higher doses (Crouzet et al., 2010; Puglisi, 2012; Romdhane et al., 2019). The application of a.i alone or in a mixture had no strong impact on microbial community composition. This might be due to the large variability observed between biological replicates within given experimental conditions and to the limited effect of the pesticide. Among all the OTUs, just three affiliated to the Acidobacteria were responsive to one pesticide (CLO). This phylum is very abundant and ubiquitous in many ecosystems, and takes part to various metabolic pathways like carbon and nitrogen cycle (Kalam et al., 2020). Previous studies described this phyla to be a good biological indicator of land-use change from forest to farmland because of its sensitivity to various toxic metals, potentially deriving from metal containing pesticides (Kim et al., 2021).
The three a.i used in our study are considered from EFSA as showing low toxicity towards soil living organisms (European Food Safety Authority (EFSA) et al., 2018a, 2018b, European Food Safety Authority (EFSA) 2018), but nothing is known about the effect of these three compounds applied in combination. Contrary to our expectation, the exposure of soil microcosms to the mixture of these three a.i at the agronomic application rate and at 10-fold this rate did not affect the measured endpoints differently than exposure to the respective single a.i. Studies on the effect of pesticide mixtures on the soil microbial community are scarce (Baćmaga et al., 2015; Joly et al., 2012, 2015; Schuster and Schröder, 1990), and the responses are variable. Based on knowledge for other endpoints (Cedergreen 2014), the probability of detecting synergistic effects is low. Hence, we can conclude that the present experiment confirms the findings of Cedergreen (2014) as no exceptionally stronger effect was found in the mixture treatments (pointing at a potential synergistic interaction) compared to the single a.i. treatments.
A compounded disturbance is the phenomenon in which ecosystems are subjected to multiple disturbances occurring over time. The effect of the initial disturbance may have legacy effect and consequently alter the ability of the microbial community to cope with the following ones. The responses to the following stresses are variable (Vinebrooke et al., 2004): i.e. the exposure to multiple perturbations of the same nature might lead to microbial selection and high community resilience to that specific stress (Calderón et al., 2018). However, if the subsequent perturbation is of a different nature, the disturbed community might either become more resistant to the following stress because of community priming (Rillig et al., 2015), or more sensitive (Calderón et al., 2018) because of negative species co-tolerance (Vinebrooke et al., 2004). In our study, the application of different a.i to a previously environmentally disturbed community had no apparent further impact on the microbial communities. Responses to a.i application were few and transient, but it is worth mentioning that the heat stress seemed to condition the response to the different tested a.i. This was the case of AOB or NH4+, especially to PYR treatment. According to the literature, abiotic disturbances, such as extreme temperatures show a great potential into shaping and modifying the microbial community (Bardgett and Caruso, 2020; Castro et al., 2010; Islam et al., 2020) as compared to pesticides where the effects are often variable and confined to certain endpoints.
Overall, the temperature disturbance had a major and persisting impact on microbial communities. On the opposite, the transient increase in soil humidity was without effects. Similarly, the effects of the tested a.i and their mixture on the control community were only slight and transient, and seem to impact mainly some N-cycle endpoints. We did not find any strong and persisting effect on microbial endpoints of the compounded disturbances. Future studies should be conducted in different soils with various physico-chemical properties to generalize our conclusions. Also, consideration of denitrifiers and N2O emitters would be of interest to give a more comprehensive understanding of the compounded effect of environmental disturbances and pesticides on N-cycle in soils.
Acknowledgements
Preprint version 4 of this article has been peer-reviewed and recommended by Peer Community in Ecotoxicology and Environmental Chemistry (https://doi.org/10.24072/pci.ecotoxenvchem.100177; Hanano & Mougin, 2025).
Fundings
The project leading to this publication, as well as Camilla Drocco’s PhD, have received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 956496.
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.
Data, scripts, code, and supplementary information availability
Data and scripts are available online: https://doi.org/10.5281/zenodo.10556208 (Drocco et al., 2024).
Appendices
Figure S1. Venn Diagram illustrating the number of Operational Taxonomic Units (OTUs) impacted by heat or high humidity disturbances at T30 and T60 across all pesticide treated microcosms. OTUs whose relative abundance was significantly impacted by disturbances were identified by performing a GLM model on the most abundant 336 OTUs (abundance > 0.1% and present at least in 5/5 of the samples of one Time-A.I.-Dose-Disturbance condition). Each circle represents the set of OTUs impacted by one disturbance (high humidity or heat) at one time of sampling (T30 or T60).
Table S1. Contribution of the microbial endpoints to the principal components of the PCA at T30 and T60. For each axis, values represent the percentage of variance explained by each endpoint in the dataset.
T30 | T60 | ||||
|---|---|---|---|---|---|
Axis 1 | Axis 2 | Axis 1 | Axis 2 | ||
OS | 13.5 | 0.0 | 15.0 | 0.7 | |
PD | 13.3 | 0.1 | 14.5 | 1.1 | |
SR | 13.0 | 0.1 | 15.2 | 0.0 | |
16S | 2.1 | 28.7 | 2.9 | 22.7 | |
18S | 8.0 | 17.6 | 10.6 | 12.4 | |
ITS | 8.1 | 11.2 | 9.8 | 9.5 | |
AOA | 1.6 | 10.3 | 0.0 | 8.5 | |
AOB | 9.6 | 0.3 | 5.2 | 1.7 | |
ComaA | 7.1 | 14.7 | 6.1 | 18.8 | |
ComaB | 7.1 | 15.9 | 4.7 | 23.7 | |
NH4+ | 4.0 | 0.7 | 0.9 | 1.0 | |
NO3- | 12.6 | 0.3 | 15.0 | 0.0 | |
PD: Faith’s Phylogenetic Diversity; OS: observed species; SR: Simpson’s reciprocal index.