Neonicotinoid Seed Treatments Influence Soil Nematode Taxonomic Composition and the Soil Microbial Cooccurrence Networks

In modern agriculture, pesticides are broadly used to prevent crop losses due to diseases, weeds, or insect pests. Neonicotinoids are a class of systemic, neuroactive insecticides discovered in the late 1980s (Kagabu 1996; Tomizawa and Casida 2005). These pesticides have chemical structures similar to that of nicotine and act on nicotinic acetylcholine receptors. They disrupt the central nervous system's neurotransmission, particularly in invertebrates (e.g., arthropods, annelids, and nematodes) that are selectively more sensitive to these insecticides (Tomizawa and Casida 2005). In field crops, neonicotinoids are usually used prophylactically as seed coverings to control various early-season insect pests (Goulson 2013; Labrie et al. 2020; Samson-Robert et al. 2017). Previous studies have reported their nontarget negative impacts on agriculturally beneficial organisms; for instance, insect pollinators such as honeybees (Bonmatin et al. 2015; Iwasa et al. 2004; Samson-Robert et al. 2014) and butterflies (Pisa et al. 2015; Whitehorn et al. 2018), and terrestrial invertebrates such as earthworms (Pisa et al. 2015; Qi et al. 2018). However, there are only a few studies on the nontarget effects of neonicotinoids on the composition and function of beneficial microbial communities such as bacteria (Cycoń et al. 2013; Li et al. 2018; Parizadeh et al. 2021; Yu et al. 2020), fungi (Li et al. 2018), and free-living nematodes (Bradford et al. 2020; Hopewell et al. 2017). This study evaluates the effects of neonicotinoid seed treatments on soil nematodes and nematode-bacteria cooccurrence networks.

In sustainable agriculture, which refers to food production using techniques that protect environmental and human health, it is crucial to maintain soil quality and health. Soil quality is the capacity of soil to function, in order to sustain biological productivity and maintain or improve environmental quality (air and water) and health (plant, animal and human) (Doran and Safley 1997). It is a combination of physical, chemical and biological (including microbial) indicators (Sharma et al. 2010; Trivedi et al. 2016). One of the most important factors in sustainable agriculture is the living soil system, including soil microbiota. Soil quality is directly tied to the diversity and dynamics of microbial communities inhabiting it (Bender et al. 2016; de Vries and Wallenstein 2017; Sharma et al. 2010). These communities often vary quickly in response to small environmental changes due to stress and perturbations such as pesticide application (Paterson et al. 2011), making them potential biological indicators to monitor and evaluate soil quality and ecological functions (Kennedy and Stubbs 2006; Sharma et al. 2010). Among soil microbes, eukaryotic microbes and mesofauna, especially free-living nematodes, are excellent biological indicators of ecological changes in soil (Bongers 1990; Bongers and Ferris 1999) because of their temporal stability compared with more dynamic microbial communities such as bacteria (Nannipieri et al. 1990; Neher 2001). Nematodes are one of the most abundant soil invertebrates (van den Hoogen et al. 2019). Free-living nematodes are among the main consumers of soil microbiota (Nannipieri et al. 1990; Neher 2001). They play a key role in shaping soil microbial communities, influencing soil functions and soil food web stability (Ferris 2010a; Neilson et al. 2020). However, our knowledge about the effects of neonicotinoids on these ecologically essential invertebrates is still limited to mostly laboratory-based studies on single species. A previous study observed some nontarget adverse neurological effects of neonicotinoids on the model organism nematode Caenorhabditis elegans (e.g., their impacts on the chemosensory ability of nematodes) (Hopewell et al. 2017). Another recent article reported a decrease in the overall growth and egg laying of C. elegans exposed to neonicotinoids (Bradford et al. 2020).

Studying community diversity and compositional variation is vital to understanding the effects of pesticide application on microbial community structure. However, these communities are constantly in interaction with each other. Thus, assessing the changes in their cooccurrence patterns may offer a more in-depth insight into community stability and dynamics (Allesina and Tang 2012; García-Callejas et al. 2018). Our previous study demonstrated that neonicotinoids alter soil bacterial taxonomic composition and community structure (Parizadeh et al. 2021). Because neonicotinoids are selectively more toxic to invertebrates, we hypothesized that neonicotinoids might affect bacterial community responses by changing the relationships between these bacteria and the higher trophic levels that feed on them (e.g., microarthropods and free-living nematodes). Cooccurrence networks of microbial communities help us identify the shifts in microbial cooccurrence patterns caused by neonicotinoid application, determine the keystone species in the ecosystem, and better understand the ecological processes that drive the complexity, stability, and resilience of the networks (Williams et al. 2014). Previous studies have shown that some pesticides reduce the complexity of soil bacterial networks (Fournier et al. 2020; Zhang et al. 2021a, b) but, in general, there is a lack of empirical data on the impacts of neonicotinoids and other pesticides on the cooccurrence patterns of soil microbial networks, especially at different trophic levels.

In this study, we use amplicon sequencing to (i) characterize communities of soil free-living nematodes, (ii) quantify variation in the composition and trophic functions of nematodes, and (iii) assess the changes in nematode-bacteria cooccurrence patterns in response to neonicotinoid seed treatment in a 3-year soybean-corn crop rotation. We hypothesized that neonicotinoid seed treatment (i) affects the community composition of soil free-living nematodes, (ii) decreases soil nematode diversity, (iii) modulates soil nematode taxonomic composition and trophic functions, and (iv) influences the cooccurrence patterns of bacteria and the free-living nematodes that feed on them. To address these objectives, we used 16S (Parizadeh et al. 2021) and 18S ribosomal RNA (rRNA) gene amplicon sequencings to quantify soil bacterial and nematode communities and evaluate nematode-bacteria cooccurrence networks in a 3-year soybean-corn crop rotation in Quebec, Canada.

MATERIALS AND METHODS

Study site.

In a 3-year crop rotation, we planted soybean and corn on an experimental farm at Agriculture and Agri-Food Canada in L'Acadie, Quebec, Canada (45°17′38.0″N; 73°20′58.0″W). This site has a temperate climate and clay loam soil and is located in the Canadian hardiness zone 5a (McKenney et al. 2014). Prior to the experiment, the field was a meadow and had not been treated by any type of neonicotinoids during the past 3 years. Each year in mid-May, we cultivated four replicates of non-neonicotinoid-treated (control) and neonicotinoid-treated soybean (2016 and 2018) or corn (2017) seeds. Each plot was 100 by 3 m and contained four rows. Two extra neonicotinoid-treated plots surrounded the experimental field. All seeds were covered by three fungicides (difenoconazole, metalaxyl-M, and sedaxane). The neonicotinoid-treated seeds were also coated with thiamethoxam at 0.25 mg/seed. The fields were under a no-tilling regime, and weeds were controlled by glyphosate (before and 1 month after seeding). On the corn field, we also applied NPK fertilizer (15-15-15) at 400 kg/ha before seeding and nitrogen fertilizer at 222 kg/ha (27.5%) 1 month after seeding the crops. Soil physicochemical properties (e.g., pH, electrical conductivity, and organic matter) did not differ across the field or between the growing seasons. For more details, see our previous study (Parizadeh et al. 2021).

Soil sample collection.

We retrieved 48 soil samples during each growing season: 2 samples per plot at three sampling times (in July, August, and September), for a total of 144 samples. We used a 2-cm-diameter corer to take soil samples from the upper 12-to-15-cm layer of the bulk soil (soil that does not adhere to plant roots). For each sample, we collected soil from six different spots in a zigzag pattern, at 10 cm around 6 to 10 close plants, then pooled them into one sample (400 to 500 g). We then transferred the samples to the laboratory in a cooler and stored them at −80°C for DNA extraction.

DNA extraction.

We used the MoBio/QIAGEN PowerSoil DNA isolation kit to extract DNA from the soil samples. Rather than using 0.25 g of soil for DNA extraction, as recommended by the manufacturer, we extracted DNA from 0.5 g of soil from each pooled soil sample of 400 to 500 g. We also extracted DNA twice from each sample, then pooled the samples. We performed the rest of the extraction procedure according to the user manual. We then used Qubit (Thermo Fisher Scientific) and Nanodrop (Thermo Fisher Scientific) to quantify and qualify the extracted DNA before storing them at −80°C.

DNA amplification, library preparation, normalization, and purification.

To target the soil nematodes and bacterial communities, we amplified the 18S and 16S rRNA genes, respectively, using an Agilent SureCycler 8,800. Each PCR assay (25 μl) contained 1 μl of 1:10 diluted genomic DNA, 5 μl of 5× HF buffer (Thermo Scientific), 0.5 μl of dNTPs (10 mM each), 0.75 μl of dimethyl sulfoxide, 0.25 μl of Phusion Hot Start II polymerase (Thermo Scientific), 1.0 μl of each of the two specific primers of the target gene amplification (5 mM), and 15.50 μl of double-distilled water. To amplify the 18S rRNA genes, we performed a two-step seminested PCR following Sapkota's protocol (Sapkota and Nicolaisen 2015). We first preamplified the DNA templates with the primers suggested by the protocol, NemF and 18Sr2b. NemF is a primer designed to target nematodes and other metazoans and to exclude plant and fungal DNA regarding the mismatches at the 3′ end (Sapkota and Nicolaisen 2015). The amplification started with an initial DNA denaturation at 94°C for 5 min; followed by 28 cycles of denaturation at 94°C for 30 s, annealing at 53°C for 30 s, and elongation at 72°C for 1 min; followed by 10 min of final elongation at 72°C. To generate smaller amplicons that were appropriate for sequencing, we diluted the PCR products (1:10) and used 1.0 μl of each of them as a template in the second step of the PCR, which was performed using NF1 (forward) and 18Sr2b (reverse) primers (Creer et al. 2010; Porazinska et al. 2009, 2010). The second amplification included 5 min of initial DNA denaturation at 94°C; followed by 22 cycles of denaturation at 94°C for 30 s, annealing at 58°C for 30 s, and elongation at 72°C for 1 min; and ended with a final elongation at 72°C for 10 min. We randomly distributed all samples into several 96-well PCR plates. Each PCR plate contained one positive and one negative control. The positive controls contained Ditylenchus dipsaci, Globodera rostochiensis, Meloidogyne hapla, Pratylenchus crenatus, and P. penetrans. We used nuclease-free, diethyl pyrocarbonate-treated, and autoclaved water as the PCR plate negative control. We also included negative extraction controls for the sampling bag, 50-ml tube, and extraction kit. After the DNA amplification, we electrophoresed all PCR products on a 2% agarose gel in 1× Tris-acetate-EDTA buffer, stained them with AMRESCO's EZ-Vision dye as a loading buffer (VWR Life Science), then visualized them by G:BOX gel doc (Syngene). For bacteria, we used 16S rRNA gene sequences of soil samples collected from the same field and under the same conditions from our previous study on the effects of neonicotinoid seed treatments on soil bacterial communities (Parizadeh et al. 2021).

Normalization, library preparation, and sequencing.

We used the SequalPrep PCR Normalization kit (Thermo Fisher Scientific) to normalize the PCR products. We pooled all of the amplified, normalized DNA and prepared one library per PCR plate. We then determined the concentration of each library using Qubit. For each sequencing run, we pooled and purified an equimolar concentration of each library using Ampure XP (Beckman Coulter by Thermo Fisher Scientific) based on the manufacturer's protocol. We measured the final concentration and quality of the purified DNA using Qubit and Bioanalyzer DNA analysis kit (Agilent). We denatured the 4 nM DNA using 0.2 N NaOH, then diluted it to a 14 pM library, according to MiSeq Illumina guidance. Finally, we sequenced the final DNA samples on Illumina MiSeq at Agriculture and Agri-Food Canada using a 600-cycle (2 × 300 bp paired-end) MiSeq reagent kit v3.

Bioinformatic analyses.

Taxonomic annotation.

We used BBDuk v35.85 (http://jgi.doe.gov/data-and-tools/bb-tools/) to filter the Illumina adapters. Then, we removed barcodes and primers and demultiplexed the Illumina reads. We applied the DADA2 pipeline (Callahan et al. 2016) to denoise and remove low-quality sequences, join the paired-end sequences, remove the chimeric sequences, and identify amplicon sequence variants (ASVs). To analyze the 18S and 16S rRNA gene data, we used the denoise-paired function of DADA2 wrapped by the QIIME2 v2019.7.0 plugin (Bolyen et al. 2019) and the DADA2 package v1.12.1 (Callahan et al. 2016) in R v4.0.0 (R Core Team 2020), respectively. To annotate the taxonomic identity of the ASVs, we assigned them to species-level taxonomy using a custom-trained 18S rRNA naive Bayesian classifier (Bokulich et al. 2018) trained on the 99% SILVA 132 database (Glöckner et al. 2017; Quast et al. 2013; Yilmaz et al. 2014) in QIIME2 for the nematode ASVs, and the RDP naive Bayesian classifier method implemented in DADA2 package with the SILVA 132 rRNA database for the bacterial ASVs. We used default parameter settings for all functions in these analyses, except for some modifications for the 16S rRNA gene data analysis; for more details, see our previous study (Parizadeh et al. 2021).

Quality filtering and rarefaction.

After evaluating the mock community composition in the positive controls of each dataset, we removed these samples. We then performed the following steps of sequence quality filtering to minimize the sequence artifacts derived from PCR and sequencing procedures (Acinas et al. 2005). (i) Our 18S rRNA primers were designed to capture nematodes and other metazoan sequences and exclude those of fungi and plants. To study soil nematode communities, we subset the Nematoda phylum from the 18S rRNA dataset (30.7% of total relative abundance), which eliminated 12 samples, including four negative controls. We also removed the ASVs that were not annotated to any bacterial phylum from the 16S rRNA dataset to keep only characterized bacteria phyla (> 99.9% of total relative abundance). (ii) Then, we excluded the samples with fewer than 1,000 sequences from both datasets (18 nematode samples, including the 4 remaining negative controls, and 17 bacteria samples, including 13 negative controls). (iii) From the bacteria dataset, we removed the bacterial DNA contaminants characterized by the prevalence method (26 ASVs, probability threshold = 0.5) of the decontam package v1.1.2 (Davis et al. 2018) in R, and the chlroplast and mitochondria ASVs (0.07% of ASVs). (iv) We also filtered the rare ASVs with fewer than 10 reads from both datasets (excluding 25% of remaining nematode ASVs and 40% of remaining bacterial ASVs). (v) We cleaned the nematode dataset from a few ASVs that were likely misannotated as a nematode (eight ASVs assigned as the Bilateria genus) or were animal parasitic nematodes (one ASV assigned as the Eucoleus genus). (vi) Finally, we eliminated outlier samples (three nematode and four bacteria samples, including the last remaining negative control) that had a highly different composition, based on the nonmetric multidimensional scaling on Bray-Curtis dissimilarities (Bray and Curtis 1957) and Shannon diversity from the other samples.

We then investigated the ASV rarefaction curves (Supplementary Fig. S1) to choose rarefaction cutoffs that reached saturation in the rarefaction curves while keeping the most possible samples. To study the soil microbial community composition, we first rarefied the nematode sequences to 1,000 reads per sample (Supplementary Fig. S1A). We also rarefied the bacteria dataset to 10,000 reads per sample (Supplementary Fig. S1B). These cutoffs excluded no nematode samples and five bacteria samples that did not have enough sequences, as well as 20 nematode ASVs and 185 bacterial ASVs. Our final denoised and rarefied datasets consisted of 119 nematode samples (including 58 control and 61 neonicotinoid-treated samples), with an average of 25.1 ± 0.6 observed ASV richness (mean ± standard error [SE]) per sample, and 132 bacteria samples (including 67 control and 65 neonicotinoid-treated samples), with an average of 2,255 ± 29.8 observed ASV richness (mean ± SE) per sample. We then used R to analyze these datasets.

Statistical analyses.

Effects of neonicotinoid seed treatment on free-living nematode community composition and diversity.

To determine the variation in the soil nematode community structure in response to neonicotinoid seed treatment, we analyzed the rarefied nematode dataset. We first performed a permutational multivariate analysis of variance (PERMANOVA) with 999 permutations on the community matrix (model: . ∼ neonicotinoid seed treatment × host species × year × month) using the adonis2 function of the vegan package v2.5.7 (Oksanen et al. 2020) in R v4.0.3 (R Core Team 2020) to assess the relationships between nematode communities and neonicotinoid seed treatment, host species, time (year and month), and their interactions. To display the composition variation in the nematode communities, we performed principal coordinate analysis ordinations (based on Bray-Curtis dissimilarities). We used the Shannon index to estimate the nematode α diversity and applied the nonparametric Wilcoxon rank-sum test (Wilcoxon 1945) to determine whether the nematode α diversity statistically significantly varied among the control and neonicotinoid-treated samples.

Effects of neonicotinoid seed treatment on nematode taxonomic composition and functions.

To evaluate soil nematode taxonomic composition and functional variation, we first assigned the nematode families to trophic functional groups and classified them into the following categories: bacterivores, fungivores, herbivores, omnivores, and predators, according to previous studies (Cesarz et al. 2015; Ferris 1999; Yeates et al. 1993). For the families that were associated with more than one trophic group based on the references, we chose the group that was assigned to their identified genera. We used the rarefied nematode dataset (119 nematode samples, including 58 control and 61 neonicotinoid-treated samples), agglomerated at the family level, to calculate the relative abundance of each nematode trophic group and family in control and neonicotinoid-treated samples. We then compared the total relative abundance of each nematode trophic group, as well as each family, between control and neonicotinoid-treated samples using the nonparametric Wilcoxon rank-sum test and adjusted the P values using Holm's method to control the bias caused by false-positive results. After that, to assess nematode community and trophic group composition and functional variation in response to neonicotinoid seed treatment, we performed the nematode faunal analysis using the R-based nematode indicator joint analysis (NINJA) tool (Sieriebriennikov et al. 2014). NINJA uses the relative weighted abundance of nematode families and their trophic groups to calculate the nematode community-level indices of the soil food web that serve as ecological indicators of environmental disturbances. For this purpose, we first ordered the nematode families based on their life strategies on a colonizer-persister (cp) scale, ranging from cp1 (early colonizers with short life cycles, resistant to perturbations) to cp5 (persisters in undisturbed environments with long life cycles, sensitive to perturbations) (Bongers 1990; Ferris et al. 2001). We then computed the following typically used soil food web indices: (i) maturity index, representing the proportion of various cp groups in nematode communities; a stable, undisturbed soil food web shows a high number of persisters and, therefore, a high maturity index value; (ii) enrichment index, representing the responsiveness of the primary decomposers and bacterivorous and fungivorous nematodes; (iii) structure index, indicating the sensitivity of the nematode trophic groups to perturbations; and (iv) metabolic footprint of each trophic group, indicating the nematode's contribution in ecosystem services and functions based on its carbon utilization and energy flow (Bongers and Bongers 1998; Ferris 2010b; Ferris et al. 2001). Finally, we applied analysis of variance (ANOVA) tests to monitor the significant overall and temporal changes between control and neonicotinoid-treated samples.

Effects of neonicotinoid seed treatment on the cooccurrence of nematodes and bacterial families.

We used network analysis to quantify the cooccurrence of microbial communities and identify the impacts of neonicotinoid application on the interkingdom cooccurrence of soil nematodes and bacterial families. To achieve this, we used the rarefied datasets agglomerated at the family level and kept only the samples that passed the quality filtering and rarefaction steps in both datasets. Thus, out of 119 nematode and 132 bacteria samples, we ended up analyzing 109 samples, which included 28 nematode and 199 bacterial families in total. We then subset control and neonicotinoid-treated samples (54 and 55 samples, respectively) and filtered families that were present in fewer than five samples. We produced one correlation matrix per treatment based on the cooccurrence and relative abundance of nematode and bacterial families by applying the sparCC algorithm (Friedman and Alm 2012) using the SpiecEasi package v1.1.0 (Kurtz et al. 2015) in R v4.0.3 (R Core Team 2020) as a measure of cooccurrence, where positive and negative correlations indicate positive and negative cooccurrence relationships, respectively, and examined the statistical significance of the correlations using bootstrapped estimates with 999 permutations. Significant correlations were defined as the pairs of families with an absolute correlation coefficient threshold > 0.3 and adjusted P value (Benjamini and Hochberg 1995) < 0.01 that had occurred together at least once (109 families in control and 106 families in neonicotinoid-treated samples). For each treatment, we constructed one ecological network of significant correlations, with nodes representing microbial families and edges (links between nodes) representing the correlations between them (Albert and Barabási 2002), using the igraph package v1.2.5 (Csardi and Nepusz 2006) in R v4.0.3 (R Core Team 2020). To determine whether the observed networks were significantly different from a network structure expected by chance, we created null models of random networks by randomly generating communities with the same number of nodes (microbial families) but randomizing the links among nodes. Thus, for each treatment, we simulated null models of 999 random networks by randomizing the cooccurrence of pairs of microbial families while keeping the same number of families and their cooccurrences (i.e., nodes and edges of the network). To assess whether neonicotinoid application induced changes in the nematode-bacteria cooccurrence networks, we calculated the following topological network properties for the whole network and for individual network nodes (i.e., families) using the functions provided by the igraph package v1.2.5 (Csardi and Nepusz 2006), designated as network-level metrics and network node metrics. Network-level metrics consisted of (i) motifs: subgraphs representing patterns of interactions between families (Milo et al. 2002) and (ii) modularity: the degree of division of a network into modules containing cooccurring taxa (Newman 2006). Network node metrics consisted of (i) betweenness centrality: the number of times a node serves as a bridge within the shortest path of the other nodes in a network (Freeman 1977); the nodes (microbial families) with the highest betweenness centrality might be identified as keystone microbial families (Agler et al. 2016); and (ii) coreness: a measure to determine whether each node belongs to a more or less densely connected part of the network (Batagelj and Zaversnik 2003). We used a Z test that compared the observed value of each metric to the distribution of the randomly generated metrics from null models to investigate whether the topological network properties varied significantly among the observed and random networks. We assessed the Z score (distance of the metric value from null model mean divided by null model standard deviation), rank (of the observed metric compared with 999 null models), and a two-tailed P value (the observed metric rank divided by the total number of random networks) for each network property against the null models.

RESULTS

Impacts of neonicotinoid seed treatment on nematode community composition variation.

Neonicotinoid seed treatment significantly affected nematode community composition (2.0% of compositional variation explained, PERMANOVA P < 0.001) (Table 1; Fig. 1). Host species and year also significantly explained some variation in community composition (2.6 and 3.0%, PERMANOVA P < 0.001) (Table 1; Fig. 1). However, the interactions among these parameters had no significant effect on the communities. Nematode α diversity was not significantly different between the control and neonicotinoid-treated samples (Shannon index mean ± SE 2.3 ± 0.05 versus 2.3 ± 0.06, P value = 0.94).

TABLE 1 Drivers of the soil nematode community composition variation in response to neonicotinoid seed treatment in a 3-year soybean-corn crop rotation in L'Acadie, Quebec, Canada^a

Variablesb	R2 (%)	F	Pr (>F)
Neonicotinoid seed treatment	2	2.5	<0.001***
Host species	2.6	3.1	<0.001***
Year	3	3.6	<0.001***
Month	NS	NS	NS
Neonicotinoid seed treatment × host species	NS	NS	NS
Neonicotinoid seed treatment × year	NS	NS	NS
Neonicotinoid seed treatment × month	NS	NS	NS
Neonicotinoid seed treatment × host species × month	NS	NS	NS
Neonicotinoid seed treatment × year × month	NS	NS	NS

^aResults are based on a PERMANOVA test on Bray-Curtis dissimilarities among samples to quantify variation in community composition (R²) explained by each factor. Significance levels for each variable are given: *** indicates P < 0.001; NS = P ≥ 0.05.

^bVariables and the interaction between variables.

View as image HTML

Download as PowerPoint

Impacts of neonicotinoid seed treatment on nematode taxonomic composition and trophic functions.

Neonicotinoid seed treatment had significant effects on the taxonomic composition of soil nematodes by modulating the relative abundance of certain families from different trophic groups. Bacterivores were the dominant trophic group in both control and neonicotinoid-treated samples (52.4 and 51.1%, respectively, of total relative abundance). Although not necessarily statistically significant, three bacterivorous families—Rhabditidae, Alaimidae, and Cephalobidae—were the most abundant nematode families in control samples at 16.4, 15.4, and 9.3%, respectively, of total relative abundance. Neonicotinoid-treated samples were also dominated by the same three families (Alaimidae 23.7%, Rhabditidae 12.1%, and Cephalobidae 9.9% of total relative abundance). Regarding nematode trophic functional groups, there was no significant differences between control and neonicotinoid-treated samples (Fig. 2A). However, out of 28 identified nematode families, the relative abundance of 2 families was significantly higher in control than neonicotinoid-treated samples: (i) the bacterivorous family Monhysteridae was only present in control samples (total relative abundance of 0.3%, Wilcoxon adjusted P < 0.05) (Fig. 2B) and (ii) the omnivorous family Dorylaimidae had a significantly higher total relative abundance in control than neonicotinoid-treated samples (total relative abundance of 2.5 versus 1.1%, Wilcoxon adjusted P < 0.05) (Fig. 2B). Overall, the nematode maturity index slightly increased in the neonicotinoid-treated samples compared with the control (mean ± standard deviation [SD] 2.84 ± 0.52 versus 2.54 ± 0.51, P < 0.05) but no temporal variation with years or months within each year was observed (P > 0.05). Conversely, the omnivore's metabolic footprint was reduced in the neonicotinoid-treated samples compared with the control (mean ± SD 0.07 ± 0.11 versus 0.13 ± 0.08, P < 0.05), whereas no difference was found among the other trophic groups. The nematode faunal analysis did not detect any overall or interannual impact of neonicotinoid seed treatment on the enrichment and structure indices of the soil food web. According to the standard faunal profile (Ferris et al. 2001), soil nematode families were highly enriched and structured in both control and neonicotinoid-treated samples (Supplementary Fig. S2).

Download as PowerPoint

Neonicotinoid seed treatment effects on nematode-bacteria cooccurrence networks.

We observed notable differences in microbial networks between control and neonicotinoid-treated soil samples following cooccurrence network analysis. We detected 391 significant nonzero correlations in control samples, including 259 positive and 132 negative correlations, and 1,110 significant nonzero correlations in neonicotinoid-treated samples (2.8 times more than control), including 579 positive and 531 negative correlations (|correlation coefficient| > 0.3, Benjamini and Hochberg [BH]-adjusted P < 0.01) (Benjamini and Hochberg 1995). The interkingdom cooccurrence relationships among nematodes and bacteria were affected by neonicotinoid seed treatment. In the control network, we observed three dominant bacterivorous nematode families (Cephalobidae, Neodiplogasteridae, and Rhabditidae) showing significant positive cooccurrence relationships and two dominant bacterivorous nematode families (Alaimidae and Neodiplogasteridae) showing significant negative cooccurrence relationships with several bacterial families (Fig. 3; Table 2). However, no significant cooccurrence between bacterivorous nematodes and bacteria was detected in the neonicotinoid-treated network. Overall, cooccurrence of bacterial families with different nematode trophic groups, including bacterivorous, fungivorous, and herbivorous nematode families (mostly from the cp1 to cp3 nematodes), was observed in the control network, while cooccurrence between omnivorous and predator nematode families (mostly cp4 nematodes) was detected in the neonicotinoid-treated network (Table 2).

Download as PowerPoint

TABLE 2 Interkingdom cooccurrences of nematode families from different trophic groups and bacterial families in control and neonicotinoid-treated soil samples in a 3-year soybean-corn crop rotation in L'Acadie, Quebec, Canada^a

			Correlations in network
			Control		Neonicotinoid-treated
Trophic group	Nematode families	cp	Positive	Negative	Positive	Negative
Bacterivore	Alaimidae	4	−	✓	−	−
	Cephalobidae	2	✓	−	−	−
	Neodiplogasteridae	1	✓	✓	−	−
	Rhabditidae	1	✓	−	−	−
Fungivore	Aphelenchoididae	2	−	✓	−	−
Herbivore	Merliniidae	3	−	✓	−	−
	Pratylenchidae	3	✓	−	−	−
	Tylenchidae	2	✓	✓	✓	−
Omnivore	Dorylaimidae	4	−	✓	−	−
	Nordiidae	4	−	−	✓	−
Predator	Mononchidae	4	−	−	−	✓

^aTrophic groups (Yeates et al. 1993) and colonizer-persister (cp) classes (Bongers 1990) of the nematode families that cooccurred with bacterial communities are presented. Correlations are calculated using the SparCC algorithm and only the significant positive and negative correlations (|correlation coefficient| > 0.3, Benjamini and Hochberg-adjusted P < 0.01) are presented, indicated by ✓ and –, respectively.

View as image HTML

Neonicotinoid seed treatment affected the structure and taxonomic composition of the cooccurrence networks. As mentioned earlier, our cooccurrence networks had almost the same number of nodes (microbial families; 109 in control versus 106 in the neonicotinoid-treated network). However, the number of cooccurrence relationships in the networks was more than doubled in response to neonicotinoid application (391 edges in control versus 1,110 edges in neonicotinoid-treated network) (Fig. 3; Table 3). These networks consisted of nine nematode families and 100 bacterial families from 14 phyla in control samples (Fig. 3A) versus five nematode families and 101 bacterial families from 13 phyla in neonicotinoid-treated samples (Fig. 3B). Null model analysis suggested that all measured network metrics in both control and neonicotinoid-treated networks were significantly different from random networks. In terms of topological network properties, whereas the number of network motifs was significantly higher than random null models for both networks and network modularity was significantly higher than random null models in control network and lower in neonicotinoid-treated network, as shown by the difference of their Z scores, the control network had fewer motifs and a higher modularity value than the neonicotinoid-treated network (observed motifs 3,110 versus 22,726, Z scores 10.21 and 29.32, P value < 0.05; observed modularity 0.46 versus 0.09, Z scores 10.42 and −4.34, P value < 0.05) (Supplementary Fig. S3; Table 3). Moreover, although the betweenness centrality and coreness of nodes in both networks were also significantly different than those of random networks based on their Z score differences, the control network nodes showed much higher average and maximum betweenness centrality values and lower average and maximum coreness values than neonicotinoid-treated network nodes (observed betweenness centrality average 240.4 versus 57.16, maximum 480.8 versus 272.9, Z scores 6.66 and 13.69, P values < 0.05; observed coreness average 4.2 versus 13.6, maximum 8.0 versus 24.0, various Z scores, P values < 0.05) (Supplementary Table S1). Neonicotinoid seed treatment also led to shifts in the taxonomic identity of the microbial families with the highest values for different network metrics. Rhodobacteraceae (Proteobacteria phylum) in the control network and Gaiellaceae (Actinobacteria phylum) in the neonicotinoid-treated network showed the highest betweenness centralities (Supplementary Table S1). The core microbial families, those with the highest coreness values associated with each network, were dominated by the families from the Actinobacteria phylum in both cooccurrence networks (Fig. 3; Supplementary Table S1). We identified 11 bacterial families, all from the Actinobacteria phylum, with the highest coreness values in both networks, while 24 bacterial families from different bacterial phyla and dominated by Proteobacteria and Actinobacteria phyla showed the same highest coreness value only in the neonicotinoid-treated network (Supplementary Table S1). None of the nematode families were among the microbial families with the highest observed network node metrics (Supplementary Table S1).

TABLE 3 Topological properties of cooccurrence networks and their respective null models obtained from microbial cooccurrence relationships in control and neonicotinoid-treated soil samples in a 3-year soybean-corn crop rotation in L'Acadie, Quebec, Canada

		Cooccurrence network	Null model
Metrics	Treatments	Observed	Z score	P value	Mean	SDa
Motifs	Control	3,110	10.21	<0.05	2,650.68	44.98
	Neonicotinoid-treated	22,726	29.32	<0.05	19,947.6	94.77
Modularity	Control	0.46	10.42	<0.05	0.29	0.02
	Neonicotinoid-treated	0.09	−4.34	<0.05	0.14	0.01
Nodes	Control	109	−	−	−	−
	Neonicotinoid-treated	106	−	−	−	−
Edges	Control	391	−	−	−	−
	Neonicotinoid-treated	1,110	−	−	−	−

^aStandard deviation.

View as image HTML

DISCUSSION

Neonicotinoid seed treatment had a small but significant impact on soil nematode community composition in a soybean-corn agroecosystem, which supports our first hypothesis. The effect of neonicotinoids was similar in magnitude to the variation explained by host species (corn versus soybean) or year. These results corroborate other studies showing nontarget impacts of neonicotinoids on different organisms, including nematodes (Seres et al. 2016), fungi (Hannula et al. 2019), and bacterial communities (Hannula et al. 2019; Parizadeh et al. 2021; Tarlera et al. 2008; Wieland et al. 2001). However, in contrast to our second hypothesis, there was no significant change in the nematode α diversity due to neonicotinoid seed treatment, indicating that neonicotinoid application changed the relative abundance and contribution of certain nematode families in soil without reducing the total number of species present.

Our results suggest that some bacterivorous (Monhysteridae) and omnivorous (Dorylaimidae) nematode families may be more responsive and sensitive to neonicotinoid seed treatment, supporting our third hypothesis about how neonicotinoid seed treatment affects the taxonomic composition of soil nematodes. Although these shifts in abundance are relatively small, they do correspond to a large shift in the abundance of certain rare taxa such as Monhysteridae. Intriguingly, many studies have used the Dorylaimidae family as an indicator of environmental disturbance because of its sensitivity to different stresses, including pesticide application (Gomes et al. 2003; Sohlenius and Sandor 1989; Thomas 1978). On the other hand, based on the faunal analysis of the soil food web, neither the nematode enrichment index nor the structure index was affected by neonicotinoid application. However, even with similar soil food web indices, the difference in nematode community composition could lead to functional consequences. For example, we observed a decrease in the omnivore's metabolic footprint in response to neonicotinoid seed treatment, indicating less contribution of this trophic group in ecosystem services and functions (such as regulating populations of opportunists) in the neonicotinoid-treated samples. A recent study on the impacts of insecticide application and pest management practices on the soil nematodes also showed no significant changes in nematode species richness but some effects on the metabolic footprint of specific trophic groups (Yang et al. 2020). Fiscus and Neher (2002) have studied in detail the direct and indirect effects of physical and chemical agricultural disturbances on soil nematodes and revealed that they differ greatly among genera. Surprisingly, the maturity index, a key index to monitor the ecological effect of disturbance on nematode communities and soil food web stability, was slightly higher in the neonicotinoid-treated samples. This is probably explained by the relative importance of the Alaimidae family in our dataset. This family is a relatively high-cp-value nematode (cp4) and has been reported to be positively sensitive to chemical disturbance (Neher 2001). Moreover, it is important to note that univariate indicators such as the maturity index are often poor indicators of subtle changes in nematode communities (Wall et al. 2002). Overall, we suggest that the soil ecosystem of our field was already in a relatively stable, enriched, and structured status, with low to moderate disturbance, as shown by the results of the faunal analysis. The fact that the field was a meadow for the 3 years before the experiment may explain the stability of the soil food web. Taken together, these results suggest that the effect of neonicotinoid application on soil nematodes could be specific to certain families and trophic groups and, thus, not evident when looking at community-wide measures such as taxon richness. It is also possible that this study did not capture the full nematode diversity and that the effects of neonicotinoids on nematode communities were slightly underestimated. Because we aimed to study the cooccurrence of nematode and bacterial communities in this work, it was important to extract the DNA of bacteria and nematodes from the same soil samples. To avoid stochasticity and increase the chance of capturing the whole diversity, we extracted DNA twice from composite samples. However, it is possible that this yielded slightly lower diversity than the classical Baermann pans extraction, which also has its bias. Further research on nematode functions is needed to better assess their variation in response to neonicotinoid seed treatment.

Ecological network analysis indicated that neonicotinoid seed treatment modulated nematode-bacteria cooccurrence patterns. This supported our last hypothesis. Despite the larger number of cooccurrence relationships of microbial families in the neonicotinoid-treated than the control network, these relationships were dominated by cooccurrences among bacterial families, and there were far fewer interkingdom cooccurrences among nematodes and bacteria in the neonicotinoid-treated network. A shift in the trophic group of the nematodes cooccurred with bacteria from bacterivorous nematodes in the control network to an omnivorous nematode in the neonicotinoid-treated network might be related to the variation in the taxonomic composition of the bacterial families in each of these networks. It has been hypothesized that measures of network structure such as network modularity are related to the sensitivity and resilience of networks faced with environmental stress and that networks with lower modularity are more sensitive to stress and more likely to spread disturbance across the network (Kharrazi et al. 2020). The lower network modularity observed in the neonicotinoid-treated network could suggest that this network has lower stability and resilience to stress than the control network. Conversely, a recent study reported a decrease in the number of nodes and edges and network complexity with increasing neonicotinoid concentrations (Zhang et al. 2021b). However, this study used much higher neonicotinoid concentrations and different soil types and physicochemical properties (Zhang et al. 2021b) compared with our study. Our results have generated hypotheses about nematode and bacterial interactions that will need to be followed up by future studies that directly characterize interactions among nematodes and bacteria and measure the stability and resilience of networks.

Although we observed shifts in nematode-bacteria cooccurrence patterns in soil exposed to neonicotinoid pesticides, there was no adequate evidence to confirm whether the variation in bacterial communities was related to the impacts of neonicotinoids on nematodes, as higher trophic groups that feed on bacteria or the direct effects of neonicotinoids on bacteria that can also affect nematodes. Thus, it is still not clear whether the nontarget effects of neonicotinoid seed treatment on bacterial community variation depend on the top-down effects of nematode predation or the other higher trophic groups such as microarthropods that feed on bacteria (Staley et al. 2015; Thakur and Geisen 2019). Another scenario, also suggested by a past study (Bradford et al. 2020), is that bottom-up effects of neonicotinoids on bacteria that, in turn, influence the nematodes by limiting the availability of bacteria as their food resources are what drive the patterns we observe. Furthermore, laboratory-based studies to examine how neonicotinoid seed treatment affects the behavioral interaction between mock bacterial communities and single species of nematodes are needed to investigate this matter. In addition to shifts in overall network structure, we found that several bacterial families, mostly from the Proteobacteria and Actinobacteria phyla, have the highest coreness values in the neonicotinoid-treated network but not in the control network. Our previous study showed that the bacteria favored by neonicotinoid seed treatment were dominated by Actinobacteria, including genera that are potentially involved in the biodegradation of neonicotinoids, and those that were negatively affected by neonicotinoid seed treatment were dominated by Proteobacteria, including the plant-growth-promoting rhizobacteria (Parizadeh et al. 2021). The effects of neonicotinoid seed treatment on these keystone bacteria could, in turn, influence the other members of the network, given their core position within the network of ecological interactions (Kitsak et al. 2010). Future research to evaluate the role of other beneficial eukaryotic microbial communities such as microarthropods in the microbial cooccurrence networks will also be required to improve our understanding of the impacts of pesticides on microbial trophic relationships.

Conclusions.

Despite the essential role of nematodes as bioindicators of soil quality and ecological functions, there have been only a few studies that focused on evaluating the effects of neonicotinoids on them, which were laboratory-based studies conducted on single species. As far as we know, the present study is the first experimental design that simulates real farming conditions while quantifying neonicotinoid effects on nematode communities. This work advances our understanding of the impacts of neonicotinoid seed treatment on soil nematode communities and the nematode-bacteria cooccurrence networks. We observed that neonicotinoids influence soil nematode community compositions and some families from different trophic groups, as well as nematode-bacteria cooccurrence networks. Future research to evaluate the role of other beneficial eukaryotic microbial communities such as microarthropods in the microbial cooccurrence networks and to investigate the effects of neonicotinoid seed treatment on top-down or bottom-up regulations of microbial communities will be required to improve our understanding of the impacts of pesticides on microbial trophic relationships.

Availability of data and materials.

We have deposited the raw sequences at the NCBI Sequence Read Archive (SRA accession number PRJNA781370). Our scripts to perform the current study analyses are available in the following GitHub repository: https://github.com/memoll/acadie_18s.

ACKNOWLEDGMENTS

We thank Michel Fortin, Éléonore Tremblay, Katherine Bisaillon, and other colleagues at Agriculture and Agri-Food Canada (AAFC) for assistance in crop cultivation and sample collection, processing, and sequencing; Gaston Mercier (AAFC) for performing soil physicochemical analysis; Etienne Lord and Pierre-Yves Véronneau (AAFC) for their help in data processing; Mathieu Landry (Université du Québec à Montréal) for his great ideas on network analyses; and Geneviève Labrie (Centre de recherche agroalimentaire de Mirabel) and Annie-Ève Gagnon (AAFC) for their valuable scientific advice and suggestions. This article was released as a preprint at SSRN (Parizadeh et al. 2022).