RESEARCHOpen Access icon OPENOpen Access license

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

    Authors and Affiliations
    • Mona Parizadeh1 2
    • Steven W. Kembel1
    • Benjamin Mimee2
    1. 1Département des Sciences Biologiques, Université du Québec à Montréal, 141 Avenue du Président-Kennedy, Montréal, Québec, H2X 1Y4, Canada
    2. 2Agriculture and Agri-Food Canada, 430 Gouin Boulevard, Saint-Jean-sur-Richelieu, Québec, J3B 3E6, Canada


    Neonicotinoid insecticides are widely used to control early-season and foliar-feeding pests. Some studies have revealed their nontarget impacts on pollinators and other invertebrates but few investigated their effects on soil microbiota. Given the crucial role of soil prokaryotic and eukaryotic microbial communities in agroecosystem regulation and their contribution to soil fertility, it is critical to understand their structure and changes in response to disturbances such as pesticide application. Among these communities, free-living nematodes have the potential to indicate the ecological changes in soil caused by environmental stress and have a key role in forming and modulating soil microbial composition and function by feeding on other soil microorganisms or interacting with them. Here, we used 18S ribosomal RNA gene amplicon sequencing to characterize the effects of neonicotinoids on soil nematode communities in a 3-year soybean-corn crop rotation in Quebec, Canada. We also quantified the changes in nematode-bacteria cooccurrence networks in soil exposed to neonicotinoids. We found that neonicotinoid seed treatment significantly explained variation in nematode community composition and affected the relative abundance of some nematode families, such as a decrease in the omnivorous family Dorylaimidae in neonicotinoid-treated samples. Moreover, neonicotinoids altered the patterns of nematode-bacteria cooccurrence, including the structure and taxonomic composition of the networks. However, it is unclear whether neonicotinoids affected bacterial cooccurrence networks directly, or indirectly by affecting nematodes that feed on bacteria. Further research is needed to understand how neonicotinoids affect nematodes and the role of nematodes in microbial network variation in soil exposed to neonicotinoids.

    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.


    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 ( 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.


    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, Canadaa

    Fig. 1.

    Fig. 1. Soil nematode community composition variation in response to neonicotinoid seed treatment. Principal coordinate analysis on Bray-Curtis dissimilarities demonstrates the composition variation of soil nematode communities between control (n = 58) and neonicotinoid-treated (n = 61) samples in a 3-year soybean (2016 and 2018) and corn (2017) rotation in L'Acadie, Quebec, Canada. Nematode communities varied among control (blue points) and neonicotinoid-treated (pink points) samples. Ellipses are shaded based on treatment (gray for control and red for neonicotinoid-treated samples) and represent a 95% confidence interval.

    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).

    Fig. 2.

    Fig. 2. Changes in soil nematode families and their relative abundance in response to neonicotinoid seed treatment. Changes in A, the total relative abundance of nematode families based on their trophic groups and B, the relative abundance of each nematode family in response to neonicotinoid seed treatment in a 3-year soybean-corn crop rotation in L'Acadie, Quebec, Canada. Box and whisker plots show the distribution of the nematodes in each trophic group (A) or family (B) based on their relative abundance (shown as square root transformed on the y-axis). Each color indicates one trophic group (blue: bacterivores, yellow: fungivores, red: herbivores, green: omnivores, and pink: predators). Boxes with lighter colors represent control, while darker boxes show the neonicotinoid-treated samples. In each box, the lower and upper hinges indicate the first and third quartiles, and the line in the middle shows the median. Lower and upper whiskers extend from the hinge to the minimum and maximum data points, respectively; at most, 1.5 times the interquartile of the hinge. Each point represents one soil sample. Asterisks indicate the groups whose relative abundance was significantly different between control and neonicotinoid-treated samples based on P values from a Wilcoxon rank-sum test. The y-axis is scaled using square roots to better visualize the groups with very low relative abundance. Taxa without taxonomic assignments at the assessed level have been considered in calculating relative abundance but have not been shown.

    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).

    Fig. 3.

    Fig. 3. Interkingdom nematode-bacteria cooccurrence network analysis. Cooccurrence networks among nematode (triangle) and bacterial (circle) communities in A, control (n = 54) and B, neonicotinoid-treated (n = 55) soil samples in a 3-year soybean-corn crop rotation in L'Acadie, Quebec, Canada. Each node represents one microbial family. Node size corresponds to microbial relative abundance and its color shows the microbial phylum. Only nodes with absolute correlations > 0.30 and BH-adjusted P value < 0.01 based on the SparCC algorithm are represented. Edges linking two nodes represent positive (solid lines) and negative (dashed lines) cooccurrences. Unknown or unidentified families are labeled with their corresponding order.

    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, Canadaa

    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


    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.


    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:


    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).

    The author(s) declare no conflict of interest.


    • Acinas, S. G., Sarma-Rupavtarm, R., Klepac-Ceraj, V., and Polz, M. F. 2005. PCR-induced sequence artifacts and bias: Insights from comparison of two 16s rRNA clone libraries constructed from the same sample. Appl. Environ. Microbiol. 71:8966-8969. Crossref, ISIGoogle Scholar
    • Agler, M. T., Ruhe, J., Kroll, S., Morhenn, C., Kim, S.-T., Weigel, D., and Kemen, E. M. 2016. Microbial hub taxa link host and abiotic factors to plant microbiome variation. PLoS Biol. 14:e1002352. Crossref, ISIGoogle Scholar
    • Albert, R., and Barabási, A.-L. 2002. Statistical mechanics of complex networks. Rev. Mod. Phys. 74:47-97. Crossref, ISIGoogle Scholar
    • Allesina, S., and Tang, S. 2012. Stability criteria for complex ecosystems. Nature 483:205-208. Crossref, ISIGoogle Scholar
    • Batagelj, V., and Zaversnik, M. 2003. An O(m) algorithm for cores decomposition of networks. arXiv 0310049. Google Scholar
    • Bender, S. F., Wagg, C., and van der Heijden, M. G. A. 2016. An underground revolution: Biodiversity and soil ecological engineering for agricultural sustainability. Trends Ecol. Evol. 31:440-452. Google Scholar
    • Benjamini, Y., and Hochberg, Y. 1995. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57:289-300. Google Scholar
    • Bokulich, N. A., Kaehler, B. D., Rideout, J. R., Dillon, M., Bolyen, E., Knight, R., Huttley, G. A., and Caporaso, J. G. 2018. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6:90. Crossref, ISIGoogle Scholar
    • Bolyen, E., Rideout, J. R., Dillon, M. R., Bokulich, N. A., Abnet, C. C., Al-Ghalith, G. A., Alexander, H., Alm, E. J., Arumugam, M., Asnicar, F., Bai, Y., Bisanz, J. E., Bittinger, K., Brejnrod, A., Brislawn, C. J., Brown, C. T., Callahan, B. J., Caraballo-Rodríguez, A. M., Chase, J., Cope, E. K., Da Silva, R., Diener, C., Dorrestein, P. C., Douglas, G. M., Durall, D. M., Duvallet, C., Edwardson, C. F., Ernst, M., Estaki, M., Fouquier, J., Gauglitz, J. M., Gibbons, S. M., Gibson, D. L., Gonzalez, A., Gorlick, K., Guo, J., Hillmann, B., Holmes, S., Holste, H., Huttenhower, C., Huttley, G. A., Janssen, S., Jarmusch, A. K., Jiang, L., Kaehler, B. D., Kang, K. B., Keefe, C. R., Keim, P., Kelley, S. T., Knights, D., Koester, I., Kosciolek, T., Kreps, J., Langille, M. G. I., Lee, J., Ley, R., Liu, Y.-X., Loftfield, E., Lozupone, C., Maher, M., Marotz, C., Martin, B. D., McDonald, D., McIver, L. J., Melnik, A. V., Metcalf, J. L., Morgan, S. C., Morton, J. T., Naimey, A. T., Navas-Molina, J. A., Nothias, L. F., Orchanian, S. B., Pearson, T., Peoples, S. L., Petras, D., Preuss, M. L., Pruesse, E., Rasmussen, L. B., Rivers, A., Robeson, M. S., Rosenthal, P., Segata, N., Shaffer, M., Shiffer, A., Sinha, R., Song, S. J., Spear, J. R., Swafford, A. D., Thompson, L. R., Torres, P. J., Trinh, P., Tripathi, A., Turnbaugh, P. J., Ul-Hasan, S., van der Hooft, J. J. J., Vargas, F., Vázquez-Baeza, Y., Vogtmann, E., von Hippel, M., Walters, W., Wan, Y., Wang, M., Warren, J., Weber, K. C., Williamson, C. H. D., Willis, A. D., Xu, Z. Z., Zaneveld, J. R., Zhang, Y., Zhu, Q., Knight, R., and Caporaso, J. G. 2019. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37:852-857. Crossref, ISIGoogle Scholar
    • Bongers, T. 1990. The maturity index: An ecological measure of environmental disturbance based on nematode species composition. Oecologia 83:14-19. Crossref, ISIGoogle Scholar
    • Bongers, T., and Bongers, M. 1998. Functional diversity of nematodes. Appl. Soil Ecol. 10:239-251. Crossref, ISIGoogle Scholar
    • Bongers, T., and Ferris, H. 1999. Nematode community structure as a bioindicator in environmental monitoring. Trends Ecol. Evol. 14:224-228. Crossref, ISIGoogle Scholar
    • Bonmatin, J.-M., Giorio, C., Girolami, V., Goulson, D., Kreutzweiser, D. P., Krupke, C., Liess, M., Long, E., Marzaro, M., Mitchell, E. A. D., Noome, D. A., Simon-Delso, N., and Tapparo, A. 2015. Environmental fate and exposure; neonicotinoids and fipronil. Environ. Sci. Pollut. Res. 22:35-67. Crossref, ISIGoogle Scholar
    • Bradford, B. R., Whidden, E., Gervasio, E. D., Checchi, P. M., and Raley-Susman, K. M. 2020. Neonicotinoid-containing insecticide disruption of growth, locomotion, and fertility in Caenorhabditis elegans. PLoS One 15:e0238637. Crossref, ISIGoogle Scholar
    • Bray, J. R., and Curtis, J. T. 1957. An ordination of the upland forest communities of southern Wisconsin. Ecol. Monogr. 27:325-349. Crossref, ISIGoogle Scholar
    • Callahan, B. J., Paul, J. McMurdie, M. J. Rosen, A. W. Han, A. Johnson, J. o. A., and Holmes, S. P. 2016. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13:581-583. Crossref, ISIGoogle Scholar
    • Cesarz, S., Reich, P. B., Scheu, S., Ruess, L., Schaefer, M., and Eisenhauer, N. 2015. Nematode functional guilds, not trophic groups, reflect shifts in soil food webs and processes in response to interacting global change factors. Pedobiologia 58:23-32. Crossref, ISIGoogle Scholar
    • Creer, S., Fonseca, V. G., Porazinska, D. L., Giblin-Davis, R. M., Sung, W., Power, D. M., Packer, M., Carvalho, G. R., Blaxter, M. L., Lambshead, P. J. D., and Thomas, W. K. 2010. Ultrasequencing of the meiofaunal biosphere: Practice, pitfalls and promises. Mol. Ecol. 19:4-20. Crossref, ISIGoogle Scholar
    • Csardi, G., and Nepusz, T. 2006. The igraph software package for complex network research. Int. J. Complex Syst. 1695:1-9. Google Scholar
    • Cycoń, M., Markowicz, A., Borymski, S., Wójcik, M., and Piotrowska-Seget, Z. 2013. Imidacloprid induces changes in the structure, genetic diversity and catabolic activity of soil microbial communities. J. Environ. Manage. 131:55-65. Crossref, ISIGoogle Scholar
    • Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., and Callahan, B. J. 2018. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome 6:226. Crossref, ISIGoogle Scholar
    • de Vries, F. T., and Wallenstein, M. D. 2017. Below-ground connections underlying above-ground food production: A framework for optimising ecological connections in the rhizosphere. J. Ecol. 105:913-920. Crossref, ISIGoogle Scholar
    • Doran, J. W., and Safley, M. 1997. Defining and assessing soil health and sustainable productivity. Pages 1-28 in: Biological Indicators of Soil Health. C. E. Pankhurst, B. M. Doube, and V. V. S. R. Gupta, eds. CAB International, Wallingford, U.K. Google Scholar
    • Ferris, H. 1999. NEMAPLEX: The nematode information system. A Virtual Encyclopedia on Soil and Plant Nematodes. Google Scholar
    • Ferris, H. 2010a. Contribution of nematodes to the structure and function of the soil food web. J. Nematol. 42:63-67. ISIGoogle Scholar
    • Ferris, H. 2010b. Form and function: Metabolic footprints of nematodes in the soil food web. Eur. J. Soil Biol. 46:97-104. Crossref, ISIGoogle Scholar
    • Ferris, H., Bongers, T., and de Goede, R. G. M. 2001. A framework for soil food web diagnostics: Extension of the nematode faunal analysis concept. Appl. Soil Ecol. 18:13-29. Crossref, ISIGoogle Scholar
    • Fiscus, D. A., and Neher, D. A. 2002. Distinguishing sensitivity of free-living soil nematode genera to physical and chemical disturbances. Ecol. Appl. 12:565-575. Crossref, ISIGoogle Scholar
    • Fournier, B., Santos, S. P. D., Gustavsen, J. A., Imfeld, G., Lamy, F., Mitchell, E. A. D., Mota, M., Noll, D., Planchamp, C., and Heger, T. J. 2020. Impact of a synthetic fungicide (fosetyl-al and propamocarb-hydrochloride) and a biopesticide (Clonostachys rosea) on soil bacterial, fungal, and protist communities. Sci. Total Environ. 738:139635. Crossref, ISIGoogle Scholar
    • Freeman, L. C. 1977. A set of measures of centrality based on betweenness. Sociometry 40:35-41. CrossrefGoogle Scholar
    • Friedman, J., and Alm, E. J. 2012. Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8:e1002687. Crossref, ISIGoogle Scholar
    • García-Callejas, D., Molowny-Horas, R., and Araújo, M. B. 2018. Multiple interactions networks: Towards more realistic descriptions of the web of life. Oikos 127:5-22. Crossref, ISIGoogle Scholar
    • Glöckner, F. O., Yilmaz, P., Quast, C., Gerken, J., Beccati, A., Ciuprina, A., Bruns, G., Yarza, P., Peplies, J., Westram, R., and Ludwig, W. 2017. 25 Years of serving the community with ribosomal RNA gene reference databases and tools. J. Biotechnol. 261:169-176. Crossref, ISIGoogle Scholar
    • Gomes, G. S., Huang, S. P., and Cares, J. E. 2003. Nematode community, trophic structure and population fluctuation in soybean fields. Fitopatol. Bras. 28:258-266. CrossrefGoogle Scholar
    • Goulson, D. 2013. Review: An overview of the environmental risks posed by neonicotinoid insecticides. J. Appl. Ecol. 50:977-987. Crossref, ISIGoogle Scholar
    • Hannula, S. E., Kielak, A. M., Steinauer, K., Huberty, M., Jongen, R., Long, J. R. D. e, Heinen, R., and Bezemer, T. M. 2019. Time after time: Temporal variation in the effects of grass and forb species on soil bacterial and fungal communities. mBio 10:e02635-19. Crossref, ISIGoogle Scholar
    • Hopewell, H., Floyd, K. G., Burnell, D., Hancock, J. T., Allainguillaume, J., Ladomery, M. R., and Wilson, I. D. 2017. Residual ground-water levels of the neonicotinoid thiacloprid perturb chemosensing of Caenorhabditis elegans. Ecotoxicology 26:981-990. Crossref, ISIGoogle Scholar
    • Iwasa, T., Motoyama, N., Ambrose, J. T., and Michael Roe, R. 2004. Mechanism for the differential toxicity of neonicotinoid insecticides in the honey bee, Apis mellifera. Crop Prot. 23:371-378. Google Scholar
    • Kagabu, S. 1996. Studies on the synthesis and insecticidal activity of neonicotinoid compounds. J. Pestic. Sci. 21:231-239. Crossref, ISIGoogle Scholar
    • Kennedy, A. C., and Stubbs, T. L. 2006. Soil microbial communities as indicators of soil health. Ann. Arid Zone 45:287-308. Google Scholar
    • Kharrazi, A., Yu, Y., Jacob, A., Vora, N., and Fath, B. D. 2020. Redundancy, diversity, and modularity in network resilience: Applications for international trade and implications for public policy. Curr. Res. Environ. Sustain 2:100006. CrossrefGoogle Scholar
    • Kitsak, M., Gallos, L. K., Havlin, S., Liljeros, F., Muchnik, L., Eugene Stanley, H., and Makse, H. A. 2010. Identification of influential spreaders in complex networks. Nat. Phys. 6:888-893. Crossref, ISIGoogle Scholar
    • Kurtz, Z. D., Müller, C. L., Miraldi, E. R., Littman, D. R., Blaser, M. J., and Bonneau, R. A. 2015. Sparse and compositionally robust inference of microbial ecological networks. PLoS Comput. Biol. 11:e1004226. Crossref, ISIGoogle Scholar
    • Labrie, G., Gagnon, A.-È., Vanasse, A., Latraverse, A., and Tremblay, G. 2020. Impacts of neonicotinoid seed treatments on soil-dwelling pest populations and agronomic parameters in corn and soybean in Quebec (Canada). PLoS One 15:e0229136. Crossref, ISIGoogle Scholar
    • Li, Y., An, J., Dang, Z., Lv, H., Pan, W., and Gao, Z. 2018. Treating wheat seeds with neonicotinoid insecticides does not harm the rhizosphere microbial community. PLoS One 13:e0205200. Crossref, ISIGoogle Scholar
    • McKenney, D. W., Pedlar, J. H., Lawrence, K., Papadopol, P., Campbell, K., and Hutchinson, M. F. 2014. Change and evolution in the plant hardiness zones of Canada. Bioscience 64:341-350. Crossref, ISIGoogle Scholar
    • Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., and Alon, U. 2002. Network motifs: Simple building blocks of complex networks. Science 298:824-827. Crossref, ISIGoogle Scholar
    • Nannipieri, P., Grego, S., and Ceccanti, B. 1990. Ecological significance of the biological activity in soil. In: Soil Biochemistry, vol. 6. J.-M. Bollag and G. Stotzky, eds. Marcel Dekker, Inc., New York. Google Scholar
    • Neher, D. A. 2001. Role of nematodes in soil health and their use as indicators. J. Nematol. 33:161-168. ISIGoogle Scholar
    • Neilson, J. A. D., Robertson, C. J., Snowdon, E. W., and Yevtushenko, D. P. 2020. Impact of fumigation on soil microbial communities under potato cultivation in southern Alberta. Am. J. Potato Res. 97:115-126. Crossref, ISIGoogle Scholar
    • Newman, M. E. J. 2006. Modularity and Community Structure in Networks. Proc. Natl. Acad. Sci. U.S.A. 103:8577-8582. Crossref, ISIGoogle Scholar
    • Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., O'Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., Szoecs, E., and Wagner, H. 2020. vegan: Community Ecology Package. R Package v2.5.7. Google Scholar
    • Parizadeh, M., Kembel, S. W., and Mimee, B. 2022. Neonicotinoid seed treatments influence soil nematode taxonomic composition and the soil microbial co-occurrence networks. SSRN. Google Scholar
    • Parizadeh, M., Mimee, B., and Kembel, S. W. 2021. Neonicotinoid seed treatments have significant non-target effects on phyllosphere and soil bacterial communities. Front Microbiol. 11:619827. Crossref, ISIGoogle Scholar
    • Paterson, E., Sim, A., Osborne, S. M., and Murray, P. J. 2011. Long-term exclusion of plant-inputs to soil reduces the functional capacity of microbial communities to mineralise recalcitrant root-derived carbon sources. Soil Biol. Biochem. 43:1873-1880. Crossref, ISIGoogle Scholar
    • Pisa, L. W., Amaral-Rogers, V., Belzunces, L. P., Bonmatin, J. M., Downs, C. A., Goulson, D., Kreutzweiser, D. P., Krupke, C., Liess, M., McField, M., Morrissey, C. A., Noome, D. A., Settele, J., Simon-Delso, N., Stark, J. D., Van der Sluijs, J. P., Van Dyck, H., and Wiemers, M. 2015. Effects of neonicotinoids and fipronil on non-target invertebrates. Environ. Sci. Pollut. Res. 22:68-102. Crossref, ISIGoogle Scholar
    • Porazinska, D. L., Giblin-Davis, R. M., Faller, L., Farmerie, W., Kanzaki, N., Morris, K., Powers, T. O., Tucker, A. E., Sung, W., and Kelley Thomas, W. 2009. Evaluating high-throughput sequencing as a method for metagenomic analysis of nematode diversity. Mol. Ecol. Resour. 9:1439-1450. Crossref, ISIGoogle Scholar
    • Porazinska, D. L., Sung, W. A. Y., Giblin-Davis, R. M., and Thomas, W. K. 2010. Reproducibility of read numbers in high-throughput sequencing analysis of nematode community composition and structure. Mol. Ecol. Resour. 10:666-676. Crossref, ISIGoogle Scholar
    • Qi, S., Wang, D., Zhu, L., Teng, M., Wang, C., Xue, X., and Wu, L. 2018. Effects of a novel neonicotinoid insecticide cycloxaprid on earthworm, Eisenia fetida. Environ. Sci. Pollut. Res. 25:14138-14147. Google Scholar
    • Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., Peplies, J., and Glöckner, F. O. 2013. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41:D590-D596. Crossref, ISIGoogle Scholar
    • R Core Team. 2020. R: A Language and Environment for Statistical Computing, ver. 4.0. 3. R Foundation for Statistical Computing, Vienna, Austria. Google Scholar
    • Samson-Robert, O., Labrie, G., Chagnon, M., and Fournier, V. 2014. Neonicotinoid-contaminated puddles of water represent a risk of intoxication for honey bees. PLoS One 9:e108443. Crossref, ISIGoogle Scholar
    • Samson-Robert, O., Labrie, G., Chagnon, M., and Fournier, V. 2017. Planting of neonicotinoid-coated corn raises honey bee mortality and sets back colony development. PeerJ 5:e3670. Crossref, ISIGoogle Scholar
    • Sapkota, R., and Nicolaisen, M. 2015. High-throughput sequencing of nematode communities from total soil DNA extractions. BMC Ecol. 15:3. Crossref, ISIGoogle Scholar
    • Seres, A., Hrács, K., Gyurcsó, G., Sárospataki, M., Szakálas, J., and Nagy, P. I. 2016. Laboratory studies on the effects of a neonicotinoid-containing seed treatment product on non-target soil animals. J. Agric. Environ. Sci. 3:7-14. Google Scholar
    • Sharma, S. K., Ramesh, A., Sharma, M. P., Joshi, O. m P., Govaerts, B., Steenwerth, K. L., and Karlen, D. L. 2010. Microbial community structure and diversity as indicators for evaluating soil quality. Pages 317-358 in: Biodiversity, Biofuels, Agroforestry and Conservation Agriculture. E. Lichtfouse, ed. Springer, Dordrecht, The Netherlands. CrossrefGoogle Scholar
    • Sieriebriennikov, B., Ferris, H., and de Goede, R. G. M. 2014. NINJA: An automated calculation system for nematode-based biological monitoring. Eur. J. Soil Biol. 61:90-93. Crossref, ISIGoogle Scholar
    • Sohlenius, B., and Sandor, A. 1989. Ploughing of a perennial grass ley—Effect on the nematode fauna. Pedobiologia 33:199-210. ISIGoogle Scholar
    • Staley, Z. R., Harwood, V. J., and Rohr, J. R. 2015. A synthesis of the effects of pesticides on microbial persistence in aquatic ecosystems. Crit. Rev. Toxicol. 45:813-836. Crossref, ISIGoogle Scholar
    • Tarlera, S., Jangid, K., Ivester, A. H., Whitman, W. B., and Williams, M. A. 2008. Microbial community succession and bacterial diversity in soils during 77,000 years of ecosystem development. FEMS Microbiol. Ecol. 64:129-140. Crossref, ISIGoogle Scholar
    • Thakur, M. P., and Geisen, S. 2019. Trophic regulations of the soil microbiome. Trends Microbiol. 27:771-780. Crossref, ISIGoogle Scholar
    • Thomas, S. H. 1978. Population densities of nematodes under seven tillage regimes. J. Nematol. 10:24-27. ISIGoogle Scholar
    • Tomizawa, M., and Casida, J. E. 2005. Neonicotinoid insecticide toxicology: Mechanisms of selective action. Annu. Rev. Pharmacol. Toxicol. 45:247-268. Crossref, ISIGoogle Scholar
    • Trivedi, P., Delgado-Baquerizo, M., Anderson, I. C., and Singh, B. K. 2016. Response of soil properties and microbial communities to agriculture: Implications for primary productivity and soil health indicators. Front. Plant Sci. 7:990. Crossref, ISIGoogle Scholar
    • van den Hoogen, J., Geisen, S., Routh, D., Ferris, H., Traunspurger, W., Wardle, D. A., de Goede, R. G. M., Adams, B. J., Ahmad, W., Andriuzzi, W. S., Bardgett, R. D., Bonkowski, M., Campos-Herrera, R., Cares, J. E., Caruso, T., de Brito Caixeta, L., Chen, X., Costa, S. R., Creamer, R., Mauro da Cunha Castro, J., Dam, M., Djigal, D., Escuer, M., Griffiths, B. S., Gutiérrez, C., Hohberg, K., Kalinkina, D., Kardol, P., Kergunteuil, A., Korthals, G., Krashevska, V., Kudrin, A. A., Li, Q., Liang, W., Magilton, M., Marais, M., Martín, J. A. R., Matveeva, E., Mayad, E. H., Mulder, C., Mullin, P., Neilson, R., Nguyen, T. A. D., Nielsen, U. N., Okada, H., Rius, J. E. P., Pan, K., Peneva, V., Pellissier, L., Carlos Pereira da Silva, J., Pitteloud, C., Powers, T. O., Powers, K., Quist, C. W., Rasmann, S., Moreno, S. S., Scheu, S., Setälä, H., Sushchuk, A., Tiunov, A. V., Trap, J., van der Putten, W., Vestergård, M., Villenave, C., Waeyenberge, L., Wall, D. H., Wilschut, R., Wright, D. G., Yang, J.-i., and Crowther, T. W. 2019. Soil nematode abundance and functional group composition at a global scale. Nature 572:194-198. Crossref, ISIGoogle Scholar
    • Wall, J. W., Skene, K. R., and Neilson, R. 2002. Nematode Community and trophic structure along a sand dune succession. Biol. Fertil. Soils 35:293-301. Crossref, ISIGoogle Scholar
    • Whitehorn, P. R., Norville, G., Gilburn, A., and Goulson, D. 2018. Larval exposure to the neonicotinoid imidacloprid impacts adult size in the farmland butterfly Pieris brassicae. PeerJ 6:e4772. Crossref, ISIGoogle Scholar
    • Wieland, G., Neumann, R., and Backhaus, H. 2001. Variation of microbial communities in soil, rhizosphere, and rhizoplane in response to crop species, soil type, and crop development. Appl. Environ. Microbiol. 67:5849-5854. Crossref, ISIGoogle Scholar
    • Wilcoxon, F. 1945. Individual comparisons by ranking methods. Biometrics 1:80. Crossref, ISIGoogle Scholar
    • Williams, R. J., Howe, A., and Hofmockel, K. S. 2014. Demonstrating microbial co-occurrence pattern analyses within and between ecosystems. Front Microbiol. 5:358. Crossref, ISIGoogle Scholar
    • Yang, B., Chen, Q., Liu, X., Chen, F., Liang, Y., Qiang, W., He, L., and Ge, F. 2020. Effects of pest management practices on soil nematode abundance, diversity, metabolic footprint and community composition under paddy rice fields. Front. Plant Sci. 11:88. Crossref, ISIGoogle Scholar
    • Yeates, G. W., Bongers, T., De Goede, R. G., Freckman, D. W., and Georgieva, S. S. 1993. Feeding habits in soil nematode families and genera-an outline for soil ecologists. J. Nematol. 25:315-331. ISIGoogle Scholar
    • Yilmaz, P., Wegener Parfrey, L., Yarza, P., Gerken, J., Pruesse, E., Quast, C., Schweer, T., Peplies, J., Ludwig, W., and Glöckner, F. O. 2014. The SILVA and “All-Species Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res. 42:D643-D648. Crossref, ISIGoogle Scholar
    • Yu, B., Chen, Z., Lu, X., Huang, Y., Zhou, Y., Zhang, Q. i, Wang, D., and Li, J. 2020. Effects on soil microbial community after exposure to neonicotinoid insecticides thiamethoxam and dinotefuran. Sci. Total Environ. 725:138328. Crossref, ISIGoogle Scholar
    • Zhang, H., Song, J., Zhang, Z., Zhang, Q., Chen, S., Mei, J., Yu, Y., and Fang, H. 2021a. Exposure to fungicide difenoconazole reduces the soil bacterial community diversity and the co-occurrence network complexity. J. Hazard. Mater. 405:124208. Crossref, ISIGoogle Scholar
    • Zhang, H., Zhang, Z., Song, J., Mei, J., Fang, H., and Gui, W. 2021b. Reduced bacterial network complexity in agricultural soils after application of the neonicotinoid insecticide thiamethoxam. Environ. Pollut. 274:116540. Crossref, ISIGoogle Scholar

    Current address for M. Parizadeh: Department of Physiology & Pharmacology, Cumming School of Medicine, University of Calgary, Calgary, AB, Canada.

    Author contributions: M.P., B.M., and S.W.K. conceived and designed the study; B.M. and S.W.K. obtained the funding; M.P. collected and analyzed data and wrote the manuscript; and all authors critically reviewed and edited the manuscript.

    Funding: This research was funded by Agriculture and Agri-Food Canada (grant AAFC J-001765 to B. Mimee), the Natural Sciences and Engineering Research Council of Canada Discovery Grants program (grant 213619) and the Canada Research Chair (to S. W. Kembel) and, in part, by support provided by Compute Canada and Calcul Quebec.

    The author(s) declare no conflict of interest.