Pesticide effects on soil fauna communities—A meta‐analysis

The wide use of pesticides across the globe to improve crop yield raises concerns regarding their impacts on biodiversity (Bernhardt et al., 2017; Pelosi et al., 2021; Wang et al., 2020). Despite recent initiatives to reduce the use of pesticides in many countries (e.g. European Commission, 2022), it continues to increase globally (Sharma et al., 2019), with a broad range of substances still being applied, and new substances continuously put on the market (Wang et al., 2020). It is, thus, crucial to synthesize the effects of pesticides on biodiversity and identify the most detrimental scenarios to gain a general understanding of their consequences for biodiversity. So far, most syntheses of the ecological consequences of pesticide use have focused on individual and population levels, on specific taxonomic groups such as pollinators and natural enemies (Desneux et al., 2007; Douglas & Tooker, 2016), or on specific pesticide types such as neonicotinoids (Pisa et al., 2015). To date, no quantitative synthesis has evaluated pesticide impacts on soil communities. Soil organisms represent about a quarter of global biodiversity (Bardgett & van der Putten, 2014) and play important roles in multiple ecosystem services (Schuldt et al., 2018; Soliveres et al., 2016; Wagg et al., 2014). Soil biodiversity, quality and fertility are, thus, crucial for humanity (Wall et al., 2015). However, the response of soil communities to environmental change remains understudied compared to above-ground communities (FAO, ITPS, GSBI, SCBD, and EC, 2020; Phillips et al., 2017). This is particularly the case for soil invertebrates (such as earthworms or springtails) that not only constitute a significant fraction of global biodiversity (Eisenhauer et al., 2019; FAO, ITPS, GSBI, SCBD, and EC, 2020) but are also an essential part of terrestrial food webs, with many above-ground animals, including vertebrates, depending on this resource (Barnes et al., 2023; Scherber et al., 2010; Wardle, 2002). Soil invertebrates are exposed to pesticides that often reach the soil through spraying, seed coating, fumigation, leaching and tillage practices (FAO, 2018). Several pesticides are persistent in soils (Hladik et al., 2018; Navarro et al., 2007), sometimes over the course of decades (chlordecone: Cabidoche et al., 2009). For instance, neonicotinoids can be found up to several years after application (Pisa et al., 2015; Wood & Goulson, 2017). Pesticides can induce direct toxic and sublethal effects on soil invertebrates, even at concentrations commonly found in the field (Pelosi et al., 2014; Pisa et al., 2015), and many soil invertebrates can accumulate persistent chemicals such as pesticides in their tissues (Beaumelle et al., 2017; Pelosi et al., 2021; Römbke et al., 2017). Direct toxic effects include a range of responses from the individual to the population level. For example, pesticides can cause DNA damage, alter enzyme activities, growth rates, foraging activity and reproduction of soil invertebrates (Gunstone et al., 2021; Pelosi et al., 2014). On top of these direct effects, pesticides can further affect soil invertebrates through different indirect effects mediated by shifts in resources and biotic interactions with soil microbes, plants and above-ground arthropods (Chen et al., 2013; Clements & Rohr, 2009; Desneux et al., 2007; Morgado et al., 2018; Zaller et al., 2014). A recent review summarized the results of 394 studies and found that a large proportion reports significant negative effects of pesticides on soil invertebrates (70.5% of the individual endpoints reported, Gunstone et al., 2021). Negative effects were far more frequent at the scale of individuals and populations, while soil invertebrate community responses were more variable. However, such qualitative, vote-counting approaches, cannot quantify the magnitude of pesticide effects. Quantitative synthesis is, thus, an important next step to generalize the consequences of pesticides for soil invertebrate communities across a range of environmental contexts (soil types, crops, geographic regions, etc.). Although detrimental effects of pesticides are generally expected, several factors related to the type of exposure and organisms involved can alter the consequences of pesticides for soil invertebrate communities. As a result, it remains unclear whether pesticides generally decrease soil biodiversity and if certain scenarios are more detrimental than others. Indeed, soil invertebrates respond differently to different types of pesticides, according to their mode of action and chemical structure (Biondi et al., 2012; Frampton et al., 2006; Jänsch et al., 2006; Pelosi et al., 2014; Pisa et al., 2015). Pesticides that target invertebrates (insecticides) or other organisms (fungicides and herbicides), and broad-spectrum substances that simultaneously target several pests and diseases, could, thus, elicit contrasted effects on soil communities. Additionally, real-world scenarios often involve combinations of multiple substances with different targets throughout the crop season (Pelosi et al., 2021), and the impacts of pesticides can differ depending on whether they are applied alone or in conjunction with other substances (Wood & Goulson, 2017). Furthermore, exposure conditions in terms of dose and temporal extent (short vs. long-term exposures) can greatly influence soil fauna response to pesticides, and detrimental effects on soil communities may increase with application rates and decrease over time (Amossé et al., 2018; Atwood et al., 2018). Another important moderator of the exposure and response of soil fauna communities to pesticides relates to their functional traits. Indeed, different functional groups of soil animals can elicit contrasting responses to pesticides due to different sensitivities (De Silva et al., 2010; Douglas & Tooker, 2016; Jänsch et al., 2006; Pisa et al., 2015). For example, environmental change is generally expected to have stronger effects on larger organisms (Liess & von der Ohe, 2005; Yvon-Durocher et al., 2011). In addition, soil invertebrates that live in close contact with the soil surface have been found to be more sensitive to pesticides than litter-dwelling invertebrates (Pelosi et al., 2013; Römbke et al., 2004). Such group-specific sensitivities could lead to a wide range of responses at the scale of community diversity and abundance, including neutral responses due to shifts in community composition and indirect effects (Gunstone et al., 2021; Pelosi et al., 2014). Only by addressing these potential context-dependencies, we can generate a comprehensive understanding of pesticide impacts on soil biodiversity able to inform land managers and policymakers. Overarching hypothesis (H1): pesticides generally elicit negative effects on soil fauna diversity and abundance across different types of pesticides and studies, with: (H1.2) the magnitude of pesticide effects depends on the type of pesticide applied, with: Hypotheses related to exposure conditions (H2): The dose and temporal extent of pesticide applications determine the responses of soil fauna, with: Hypotheses related to the sensitivity of different soil fauna groups (H3): Functional traits influence the sensitivity of soil fauna communities to pesticides, with: We conducted a systematic literature search to collect studies on pesticide impacts on soil fauna communities. We used a previously compiled database of studies on global change effects on soil fauna communities (Beaumelle, Thouvenot, et al., 2021; Phillips et al., 2019). The search terms are fully reported in (Phillips et al., 2019). Specifically, they included keywords to retrieve pesticide studies (such as ‘pesticide’, ‘fungicide’, ‘insecticide’ and ‘agrochemical’) and soil fauna community studies (combination of keywords of the form ‘soil’ AND (‘fauna’ OR ‘arthropod’ OR ‘invertebrate’) AND (‘community’ OR ‘biodiversity’)). From the full systematic search, by screening the abstracts and full texts, we identified 104 studies focusing on pesticide impacts (Beaumelle, Thouvenot, et al., 2021). We further screened the full texts of those 104 studies against our inclusion criteria to identify our final set of 54 included studies (PRISMA plot: Figure S2). To verify the efficacy of the search, we checked that the initial literature search retrieved a set of preidentified papers relevant to the current meta-analysis scope (Burrows & Edwards, 2002; Fountain et al., 2007; Knacker et al., 2004; Scholz-Starke et al., 2011, 2013; Velcheva et al., 2012). For each study, we collected soil fauna abundance, biomass, richness, and diversity indices along with several moderators. First, we defined the reference and treatment levels for each study based on the levels or intensities of pesticide use. For studies investigating multiple sites located along a gradient of pesticide intensities, the reference level was set as the condition with the lowest pesticide use intensity or concentration, and the treatment level was the treatment with the highest pesticide use intensity or concentration. We collected community data at the lowest taxonomic resolution reported by the study but higher than the family level (except for Enchytraeidae). When the study presented the response of several taxa or several pesticide treatments (different products, or different mixtures of pesticides), we collected them all and created separate observations for each taxon and pesticide treatment. If studies measured the same plot over multiple time points, we extracted control and treatment data from the latest time point to reflect long-term effects and control for temporal pseudoreplication (Ferlian et al., 2020; Hedges et al., 1999; Yue et al., 2018). We extracted means and variances for the reference and treatment levels from tables, text or figures (Rohatgi, 2018). To test the effects of different pesticides (H1.2), we recorded pesticide identity (active ingredients) and assigned them to a pesticide type (Table S1). Pesticide types were herbicides, fungicides, insecticides, or broad-spectrum substances (i.e. a single substance targeting several pests across taxonomic groups, e.g. fumigants). As a single study addressed the effect of a nematicide, we included that study in the insecticides category. When pesticide treatments involved multiple products combined into single, or several applications, we categorized the case study as ‘multiple substances’. We acknowledge that this categorization does not allow us to address the effect of mixtures of pesticides in its strict sense. From a soil fauna community perspective, realistic pesticide exposure can involve tank mixtures (i.e. pesticides applied together, from the same sprayer), as well as the combination of several pesticides applied over the crop season. Both scenarios can lead to significant detrimental effects on soil biodiversity, and our approach enabled us to address them together in the present analysis. We extracted pesticide application rates and the temporal extent of studies to test our hypotheses (H2.1) and (H2.2). We recorded if pesticides were applied at their recommended application rates or not based on information provided in the primary studies. When that information was missing, we compared the applied rates with the recommended application rates retrieved from the database ANSES E-Phy (ANSES, 2023) for similar crops and target pests. We also recorded if studies addressed the effects of a single pesticide application, or applied pesticides repeatedly (e.g. throughout the crop season). The temporal extent of the study reflected the time of the last sampling event after the start of the experiment (that is, the first pesticide application), with three categories: (1) ‘short-term’ studies measured soil fauna response up to 3 months after pesticide application; (2) ‘intermediate-term’ studies spanned 4 months to a year; and (3) ‘long-term’ studies from over a year to several years. It is important to note that our approach does not allow testing the temporal dynamics of pesticide effects on soil communities, as this was outside the scope of our study. Quantifying temporal effects would have required to use a different modelling framework to reduce temporal autocorrelation within each study. It would also further reduce the amount of available data, as fewer studies measure biodiversity changes over time. Instead, our approach, commonly used in other meta-analyses (Ferlian et al., 2020; Hedges et al., 1999; Yue et al., 2018), enables us to account for any variability in our mean effect sizes driven by different temporal extents in the primary studies. To quantify whether functional traits moderated the response of soil fauna communities to pesticides (H3), we first assigned soil fauna taxa to taxonomic groups defined by the Global Soil Biodiversity Atlas (Orgiazzi et al., 2016) before assigning two functional traits. Taxonomic groups were categorized into (1) two functional groups based on the presence (arthropods) or absence (enchytraeids, nematodes and earthworms) of an exoskeleton and (2) three body size categories (microfauna, mesofauna and macrofauna, Table S2). Although soil animals operate along a continuum of body size, body size categories are often used in soil ecology, as they reflect different food web compartments (Decaëns, 2010; Potapov et al., 2021; Thakur et al., 2020) and are useful for comparing the response of different soil fauna groups to global change drivers in meta-analytical approaches (Phillips et al., 2019). We chose Hedge's g as the metric for our effect sizes because we were interested in the magnitude difference between control and pesticide treatments (Koricheva et al., 2013). Hedge's g represents a standardized mean difference (i.e. the raw mean difference is divided by the pooled standard deviation of the two groups) corrected for positive bias (Viechtbauer, 2010). Hedge's g is widely used in ecology, making our results comparable to other studies. Furthermore, this metric is more suitable than the log-response ratio to incorporate the zero or close to zero values in our dataset (such as very low abundance or richness values due to pesticide treatments; Lajeunesse, 2013). For Shannon data, 11 observations from six studies lacked standard deviations (SDs) because studies did not report them. Based on our dataset that included 50 means and SDs of Shannon index responses (from nine other studies), we regressed the SD against the mean and used the slope to impute missing SDs (Lajeunesse, 2013). We note that a limitation of this approach was that the lack of data prevented the possibility of including the effects of moderators, such as pesticide type or community metric, on the relationship between means and SDs. For other diversity metrics, such as evenness, the number of studies was insufficient to use a similar approach, and observations that lacked SDs were not considered. We fitted meta-analytical models to test our hypotheses. All models included a random effect of the study identity, as data from the same study are not independent. In addition, an observation-level random effect was also used following Viechtbauer (2022), to ensure that the models did not assume identical effects for different observations within the same study. Effect sizes were weighted by the inverse of their variance in all models, giving more weight to well-replicated studies (Koricheva et al., 2013). We used Q-Q plots of the standardized deleted residuals and verified that they did not strongly depart from a normal error distribution (Viechtbauer, 2020). We also used Cook's distance plots to identify influential observations and checked that excluding them from the models yielded similar results. We tested interactive effects between moderators using likelihood ratio tests (LRTs). When interactions were not significant, we refitted the models without the interaction to test the main effect of each moderator using Wald-type Chi-square tests (QM). Main and interactive effects from the models were considered significant when p < 0.05. We identified the most detrimental scenarios based on the predicted mean effect sizes and their confidence intervals. A significant mean effect size is shown when the confidence interval does not overlap with zero. We used R software and the R package metafor (R Core Team, 2019; Viechtbauer, 2010). First, we fitted a main model of the form: Y ~ PesticideType x CommunityMetric, with Y being the effect size, PesticideType, a categorical moderator with five levels (herbicides, fungicides, insecticides, broad-spectrum and multiple substances), and CommunityMetric, a categorical moderator with two levels: abundance (densities and biomass data) or diversity (richness and diversity indices) because there were too few studies to test separately abundance, biomass, diversity indices and richness. This main model allowed us to test our overarching hypothesis (H1) by quantifying the marginal grand mean effect size across all studies, along with the mean effect sizes of different types of pesticides and community metrics (H1.1). According to (H1.2) we expected stronger negative effects of multiple substances, broad spectrum substances targeting a wide range of organisms, and insecticides that have direct toxic effects on invertebrates. To formally test (H1.2), we used a post-hoc analysis and fitted a model similar to the main model, but with four categories of pesticides: (1) multiple substances, (2) broad spectrum substances, (3) insecticides and (4) fungicides and herbicides as a group of substances that do not target soil invertebrates (Table 1). Publication bias was assessed based on the main model that included pesticide type and community metric as moderators. The funnel plot indicated no clear sign of publication bias (Figure S3). We further fitted Egger's regression between standard normal deviates of the effect sizes and precision (i.e. the inverse of their standard errors) with pesticide type and community metric as covariates interacting with precision (Koricheva et al., 2013). The results did not indicate publication bias, as none of the intercepts from Egger's regression differed significantly from zero showing no significant funnel plot asymmetry (Table S3). We further tested the influence of four moderators related to exposure conditions: ‘recommended rates’ (H2.1) and ‘temporal extent’ (H2.2), and to soil fauna functional traits: ‘body size’ (H3.1), and ‘exoskeleton’ (H3.2). Those subsequent models tested each of the four moderators in interaction with pesticide type because we expected that the influence of the moderators would depend on the type of pesticide (H1.2). Due to the lack of studies, we could not include interactions with the community metric (few diversity observations). The second-level models, therefore, addressed the combined responses of diversity and abundance. As the main model highlighted community metric as a significant moderator, we included a random effect of community metric in the second-level models to account for differences in diversity and abundance response across studies, in addition to the random effect of study identity. These second-level models were of the form: Y ~ PesticideType x Moderator, with a random effect structure: (1|CommunityMetric) + (1|StudyID/ObservationID), and Moderator being either the temporal extent, body size or exoskeleton. For recommended rates, we tested Y ~ PesticideType with the same random effect structure. Each model was based on different subsets of the data suitable to test our hypotheses (Table 1). For recommended rates, the data subset included all studies and observations of pesticides used at their recommended rates (n = 31 studies; 158 observations), and all pesticide types were represented (although broad-spectrum and fungicides had few diversity observations). For temporal extent, we excluded broad-spectrum substances that did not have long-term observations (Table 1), which removed four studies and 37 observations for that analysis. For body size, three studies and 12 observations were removed. Primarily, this was due to the observations of soil invertebrates being at a taxonomic level that inherently spanned multiple body size categories (e.g. measuring all soil invertebrates). As a result of this subsetting, broad-spectrum pesticides were removed as only one study remained (addressing their impacts on macrofauna). Similarly, for the exoskeleton model, five studies and 14 observations were removed, as the observations were not at a taxonomic level that could be identified as with or without an exoskeleton (e.g. measuring all macrofauna). However, this model did include all five pesticide categories (broad spectrum, fungicides, herbicides, insecticides and multiple substances: Table 1). Finally, we performed supplementary analyses focusing on the impact of multiple pesticides. We compared the mean effect sizes of multiple substances combining different types of pesticides versus multiple substances belonging to the same pesticide type. We also quantified the mean effect sizes of two well-represented pesticides in our dataset (glyphosate and neonicotinoids) and compared their effects when used alone or in combination with other substances. Mean effect sizes were estimated with a similar modelling approach as described above (random effect structure, Wald-type Chi-square tests and data subsets of multiple substances, glyphosate or neonicotinoid studies). Our study encompasses the results from 54 publications, with a total of 294 individual observations. Most were field experiments (n = 51 studies; 277 observations), covering a wide range of environmental contexts, although mostly located in the Northern Hemisphere (Figure S1). Studies addressed a total of 86 active substances (Table S1) belonging to different pesticide types (Figure 1): insecticides (n = 19; 116), herbicides (n = 10; 37), fungicides (n = 8; 42) and broad-spectrum pesticides (such as fumigants, n = 7; 37). Twenty studies addressed the impact of multiple substances, either applied together or as separate applications throughout the crop season. Multiple substances included combinations of several pesticide types (n = 7; 14), or several products belonging to the same pesticide type (n = 14; 48). Most studies applied pesticides at recommended application rates (n = 31, 158), but their temporal extent varied greatly, with fauna sampled from 7 days up to 22 years after pesticide application. Studies reported soil fauna community abundance (n = 47; 190), diversity indices (n = 16; 41), richness (n = 12; 28), biomass (n = 9; 27) and evenness indices (n = 3; 8). They covered 15 taxonomic groups (Figure 1; Table S2). Nematodes were the most represented (n = 21 studies; 80 observations), followed by Acari (n = 13; 50), Collembola (n = 13; 42) and earthworms (n = 12; 55). Several studies also reported the response of multiple taxa together (macro-arthropods, macrofauna, micro-arthropods or soil fauna). All pesticide types had a similar balance of observations among the most represented taxa (Figure 1). However, studies investigated insecticide effects on a much wider range of soil fauna groups (Figure 1, 16 groups) than the other pesticide types (9 for multiple substances, 7 for broad-spectrum and fungicides and 6 for herbicides). The results revealed a significant decrease in soil fauna abundance and diversity in response to pesticide use across studies. The grand mean effect size (Hedge's g) was −0.30 (CI: −0.47; −0.14; Figure 2a). Richness and diversity indices were more negatively affected than abundance and biomass of soil fauna (significant effect of community metric: QM(df = 1) = 3.94, p = 0.0471), an effect that was consistent across different pesticides (Figure 2a: non-significant pesticide type x community metric interaction, LRT(df = 4) = 3.43; p = 0.4886). Multiple substances and broad-spectrum substances significantly decreased soil fauna abundance and diversity, and insecticides significantly decreased soil fauna diversity (Figure 2a). Although the model indicated that all pesticides generally had a negative effect (non-significant effect of pesticide type QM(df = 4) = 8.35, p = 0.0797), mean effect sizes of herbicides and fungicides had confidence intervals that overlapped with zero indicating non-significant (or neutral) effects when averaging across studies. A post-hoc analysis revealed that this was due to similarities between the effects of fungicides and herbicides. Grouping together these two types of pesticides, that do not target invertebrates, revealed a significant effect of pesticide type (QM(df = 3) = 8.20, p = 0.0420), showing that broad-spectrum, multiple substances and insecticides elicited stronger negative effects compared with substances not targeting invertebrates (fungicides and herbicides). Supplementary analyses suggested stronger negative effect sizes of multiple substances combining different types of pesticides (herbicides and fungicides, herbicides and insecticides, insecticides and fungicides) than of multiple substances combining several pesticides of the same type (broad-spectrum substances, herbicides, fungicides or insecticides; Figure S4a). Although that trend, based on a subset of the data (only for recommended application rates to avoid confounding effects of application rates and type of multiple substances) was not statistically significant (QM(df = 1) = 1.04, p = 0.3086, n = 15; 30), the magnitude of the difference between mean effects was 0.43 (CI: −1.25; 0.39). We further compared the impacts of two types of widely used pesticides, glyphosate and neonicotinoids, when they were applied alone or combined with other pesticides. The results indicated a neutral mean effect size of these pesticides when applied alone (Figure S4). When glyphosate was combined with other pesticides, the mean effect size was also not significant (Figure S4b). However, we found a strong and significant negative effect of neonicotinoids applied in combination with other substances (Figure S4c). When analysing studies that applied pesticides at their recommended rates (Figure 2b, n = 31; 158), the mean effect size on soil fauna abundance and diversity had a magnitude similar to that of the main analysis (−0.30) but with confidence intervals that overlapped zero (−0.66; 0.06). Responses did not differ according to pesticide types (QM(df = 5) = 9.65, p = 0.0858), and only broad-spectrum substances caused significant declines when applied at recommended rates (Table S4). However, in contrast with the main meta-analysis (Figure 2a), this model combined diversity and abundance responses due to the lack of diversity observations (Figure 2b). Focusing on pesticide types that had enough observations to quantify the separate mean effect for diversity and abundance (multiple substances, insecticides and herbicides), we found that recommended rates of multiple substances and insecticides had significant negative effects on diversity, and neutral mean effects on abundance, while herbicides had neutral effects on both diversity and abundance (Figure 2b). The temporal extent of the study (the time between the start of the experiment and soil fauna sampling) did not significantly affect soil fauna communities (Figure 3a; QM(df = 2) = 1.26, p = 0.5320, n = 50; 257). Across pesticide types, short- and long-term studies had the lowest mean effect sizes, but only the long-term mean effect size was significantly negative (Figure 3a). Long-term studies were dominated by insecticides (n = 5; 23) and multiple substances (n = 9; 18), while long-term effects of fungicides and herbicides were only addressed in n = 2; 4 and n = 2; 3 studies and observations, respectively. Furthermore, most long-term studies focused on repeated pesticide applications (n = 14 out of 17), while short- and intermediate-term studies mostly addressed the impacts of a single pesticide application (n = 15 out of 22, and 14 out of 19 studies, in short- and intermediate term, respectively). Testing the effect of the temporal extent for single versus repeated applications separately across all pesticide types would have needed more data. For insecticides, we observed similar short-, intermediate- and long-term effects of single applications on the abundance of soil fauna (QM(df = 2) = 3.01, p = 0.2219). Single insecticide applications had significant negative effects in the short term (−0.22; CI: −0.44; −0.01; n = 7; 63) and intermediate term (−0.45; CI: −0.85; −0.04; n = 3; 11). In the long-term category, we had only two studies addressing the effects of the same active ingredient (diflubenzuron) and the effect size was neutral (0.03; CI: −0.33; 0.39; n = 2; 15). There was no evidence that body size and the presence/absence of an exoskeleton significantly modulated the response of soil fauna communities to pesticides (Figure 3). In general, we found similar responses for the three body size categories (QM(df = 2) = 2.41, p = 0.2994, n = 47; 247). Although the presence of an exoskeleton was associated with significant negative effects (effect size non-overlapping with zero, Figure 3), communities of organisms with an exoskeleton or a soft body had broadly overlapping effect sizes that were not significantly different from each other (QM(df = 1) = 1.11, p = 0.2917, n = 49; 280). None of those traits interacted with pesticide types (Body size: LRT(df = 6) = 5.45, p = 0.4877; Exoskeleton: LRT(df = 4) = 3.29, p = 0.5102). Despite raising concerns about the ecological consequences of the widespread and intensive use of pesticides in agriculture, our study is the first meta-analysis to investigate their impacts on soil fauna communities. The results confirm the overarching hypothesis (H1) t