O ﬀ -Target Stoichiometric Binding Identi ﬁ ed from Toxicogenomics Explains Why Some Species Are More Sensitive than Others to a Widely Used Neonicotinoid

: Neonicotinoids are currently licensed for use in 120 countries, making accurate nontarget species sensitivity predictions critical. Unfortunately, such predictions are fraught with uncertainty, as sensitivity is extrapolated from only a few test species and neonicotinoid sensitivities can di ﬀ er greatly between closely related taxa. Combining classical toxicology with de novo toxicogenomics could greatly improve sensitivity predictions and identify unexpectedly susceptible species. We show that there is a >30-fold di ﬀ erential species sensitivity (DSS) for the neonicotinoid imidacloprid between ﬁ ve earthworm species, a critical nontarget taxon. This variation could not be explained by di ﬀ erential toxicokinetics. Furthermore, comparing key motif expression in subunit genes of the classical nicotinic acetylcholine receptor (nAChR) target predicts only minor di ﬀ erences in the ligand binding domains (LBDs). In contrast, predicted dissimilarities in LBDs do occur in the highly expressed but nonclassical targets, acetylcholine binding proteins (AChBPs). Critically, the predicted AChBP divergence is capable of explaining DSS. We propose that high expression levels of putative nonsynaptic AChBPs with high imidacloprid a ﬃ nities reduce imidacloprid binding to critical nAChRs involved in vital synaptic neurotransmission. This study provides a clear example of how pragmatic interrogation of key motif expression in complex multisubunit receptors can predict observed DSS, thereby informing sensitivity predictions for essential nontarget species.


■ INTRODUCTION
Limiting the adverse effects of the >2 M tonnes/annum of pesticides used globally is of paramount importance for the preservation of ecosystem goods and services.Consequently, there is a critical need to better understand and predict the ecological impacts of pesticides on ecologically important nontarget species.Understanding what drives nontarget species sensitivity to agrochemicals can help to identify the species most at risk. 1 Neonicotinoids are a class of widely used insecticides that are licensed in 120 countries.Neonicotinoid compounds are known to impact on biomes throughout the globe, including through effects on pollinators, 2−4 freshwater communities, 5 and soil organisms. 6Differential species sensitivity (DSS) for individual neonicotinoids can range over 7 orders of magnitude. 7Differential species sensitivity to pesticides (and other chemicals) is dependent on both toxicokinetic (TK) 8 and toxicodynamic (TD) traits, 9 with the former determining the absorption, distribution, metabolism, and excretion (ADME) and the latter including chemical target receptor orthologue diversity, 10 ligand binding domain characteristics, and the subsequent downstream adverse outcome pathway responses. 11,12The potential effects of pesticides on critical nontarget organisms are often difficult to predict because TK and TD traits (e.g., receptor compliment) are often very poorly described.
Earthworms represent a critical nontarget group that have wide-ranging impacts on soil-based processes and provide vital soil ecosystem services. 13By combining classical toxicology with de novo toxicogenomics, we here explain the basis for DSS for the neonicotinoid imidacloprid, the most widely used insecticide worldwide, both between earthworms and a sequenced insect of known sensitivity (Drosophila melanogaster) and also between five earthworm species tested here.As well as including the standard test species (Eisenia fetida), the chosen earthworm species represent a wide phylogenetic range (i.e., Lumbricidae and Megascolecidae) and multiple ecotypes (i.e., epigeic and endogeic species).Previous to our work, the toxicity of neonicotinoids to earthworms has mainly been assessed in Eisenia fetida.However, studies conducted with insecticides have shown that earthworms can vary in sensitivity, both compared to other soil taxa 14 and also between species. 15Hence, the species in this taxon provide an ideal case study in which to assess the completeness and breadth of our combined approach for attributing a mechanistic cause to the observation of interspecies difference in sensitivity.
Species-specific sensitivity relating to TD (e.g., receptor) traits can be predicted by receptor analyses that compare primary amino acid sequences, functional domains, and individual residues of target receptors to assess the critical interactions with chemicals and ligand binding domains (LBD). 16Imidacloprid causes dysfunction by acting as a partial agonist of the nicotinic acetylcholine receptor (nAChR). 17Given that considerable electrophysiological and structural analyses have already identified specificl o o p sequences and residues critical to nAChR−imidacloprid binding, 18−20 nAChRs should be particularly amenable to a toxicogenomic predictive approach.Indeed, residue changes identified using electrophysiology that reduce nAChR responses to imidacloprid have been associated with lower imidacloprid sensitivity in ticks and aphids. 21,22Therefore, applying a predictive toxicogenomic approach to nAChRs that accounts for amino acid identity and expression in loop residues critical to LBD−imidacloprid interactions has great potential to improve neonicotinoid sensitivity predictions.

■ MATERIALS AND METHODS
Study Species and Genotyping.Five earthworm species (four Lumbricidae, Eisenia fetida, Lumbricus rubellus, Aporrectodea caliginosa, and Dendrobaena octaedra; one Megascolecidae, Amynthas gracilis) were tested for imidacloprid sensitivity.E. fetida were taken from an in-house culture, while the remaining species were collected from the field populations in the UK (L.rubellus, A. caliginosa, and D. octaedra)o rS a õ Miguel Island, The Azores (A. gracilis) (see the Supporting Information (SI) methods section for site details).Earthworm morphospecies are known to comprise distinct clades for E. fetida, 23 L. rubellus, 24,25 A. caliginosa, 26 and A. gracilis 27 but not D. octaedra. 28For this reason, the mitochondrial cytochrome oxidase 1 (CO1) locus was sequenced for 18−20 individuals to assess the population genetic homogeneity.All animals used for transcriptome expression analysis (see below) were also CO1 genotyped (see the SI Methods).
Earthworm Toxicity Testing.The full test method is detailed in the SI Methods.Briefly, tests with the five earthworm species were conducted in a natural soil (Kettering loam) amended with 3% w/w organic matter (composted bark from LBS Horticultural, Colne, UK).The soil was described in more detail by Robinson et al. 29 Earthworm numbers per replicate were modified (minimum five individuals per rep) on the basis of their availability from collection and also to avoid crowding effects. 30Imidacloprid was spiked as an aqueous solution at seven concentrations (including a control) up to 10 or 30 mg/kg dry weight soil with four replicates per treatment.The added pesticide was left for 24 h before the earthworms were added.This preincubation time allowed for an initial chemical binding to soil but avoided the potential for substantial degradation.Spiked horse manure was supplied as food.As the earthworm species differ in their temperature preference, tests were conducted at 20 ± 2 °C for E. fetida and A. gracilis and 13 ± 2 °C for A. caliginosa, L. rubellus, and D. octaedra.These temperatures were selected on the basis of the conditions used in culturing or annual median temperatures at their field collection sites and were chosen to ensure optimal cocoon production while limiting any background control mortality (as can occur at higher temperatures in temperate soil species 31 ).Earthworm survival was recorded, the laid cocoons were counted after 28 days exposure, and LC 50 and EC 50 -reproduction values were calculated in the DRC package in R 1.2.1335. 32midacloprid Uptake and Distribution.At the end of the test, earthworm tissue and a soil sample were collected for imidacloprid determination.The resource constraints meant that we could not analyze all samples, but instead we had to adopt a targeted approach.For imidacloprid in soil, we measured concentrations in 15% of samples from across the 28 day sampled soils.An analysis of these samples would provide both confirmation of no gross errors in dosing stock solution preparations (i.e., measured concentration within an order of magnitude of nominal) and also an indication of imidacloprid degradation.A general linear model (GLM) was conducted in Minitab v18 with log value to assess the concentration and temperature effects on the measured concentrations using concentration and temperature as fixed factors.
Tissue imidacloprid concentrations were measured in the tail samples of two earthworms from each replicate in all treatments for L. rubellus to assess uptake over the concentration range.To compare uptake between the species, imidacloprid was measured in two individuals per replicate in the 0 and 0.37 mg/kg treatments.As the toxicokinetics of imidacloprid was not known, the earthworms were not depurated as this takes 48 h, over which time a proportion of the chemical could have been eliminated. 33Hence, the measured concentrations contain both truly internalized chemical and also gut content associated residues.Imidacloprid measurements were conducted on the basis of the study by Woodcock et al. 2 by liquid chromatography−mass spectrometry (LC−MS) (see the SI Methods).The imidacloprid tissue concentrations between the species were compared.The log transformed data were tested for normality using a Shapiro−Wilk test before performing an one-way analysis of variance (ANOVA) with Tukey's post hoc test in Minitab v18.Toxicokinetic rates were assessed in time-series radiolabeled imidacloprid uptake and elimination studies with the most sensitive species A. gracilis, least sensitive L. rubellus, and the standard test species E. fetida (second least sensitive).These species were selected to provide comparisons of earthworms with different and similar sensitivities.Exposures to 14 C imidacloprid were conducted in the same soil as the toxicity bioassays with the labeled chemical at 0.0625 mg/kg (approximately half of the lowest EC 50 -reproduction), giving an activity of 481 kBq/kg (see the SI Methods).The activity was measured after tissue homogenization in acetonitrile by following Carter et al. 34 Individual measurements for each species were fitted without statistical transformation using a Environmental Science & Technology one compartment TK model to derived the uptake (k 1 ) and elimination (k 2 )r a t ec o n s t a n t s( ±SE) in GENSTAT.Imidacloprid half-life in tissues (T 1/2 ) and bioaccumulation factors were derived from these k 1 and k 2 rates.
To assess the imidacloprid tissue distribution, separate exposures were conducted for 35 individual L. rubellus, E. fetida,a n dA.gracilis exposed to 0.0625 mg/kg 14 Cimidacloprid for 168 h.All earthworms were depurated for 24 h on wet filter paper before being dissected to obtain nerve cord, seminal vesicle, and post clitellum body wall tissues for scintillation counting. 34Samples were pooled to give five replicates each with tissues from seven individuals (see the SI Methods).
Animal Dissection, RNA Isolation, and Transcriptome Production.Tissues were dissected from the same three (CO1 genotyped) earthworms used for the TK study as pairs showing different and similar sensitivities (L.rubellus, five animals; A. gracilis, five animals; E. fetida, six animals) and isolated RNA were prepared and indexed for paired-end multiplexed sequencing (TruSeq Stranded mRNA Library Prep Kit, Illumina, San Diego, California) before being sequenced (150 bp paired-end) using various Illumina sequencing platforms (see the SI Methods).Raw reads were trimmed using Trimmomatic24 (v.0.36) 35 and assembled using Trinity (v.2.2.0). 36An assembly was created for each species using all reads for all tissues in addition to single tissue assemblies.The resulting species lists of contigs from both the all tissue and single tissue assemblies were combined and processed using EvidentialGene: tr2aacds, mRNA Transcript Assembly Software (v.2017.12.21) 37 into the most biologically useful "best" set of mRNA and classified into primary and alternate transcripts.
Earthworm nAChR/AChBP Orthologue Identification.Paired reads (L.rubellus 437 million, E. fetida 819 million, A. gracilis 182 million, including ∼41 million L. rubellus nerve cord, ∼ 66 million E. fetida nerve cord and ∼66 million A. gracilis nerve cord) were mapped to the transcriptome assembly (RSEM v. 1.3.0),and counts were normalized with weighted trimmed mean of M-values (TMM) to calculate the relative expression levels.Earthworm nAChRs and AChBPs were determined using a Reciprocal Best Hit BLAST method 38 and sequence alignments, trimming, and trees built using Geneious (v.9.1.8).Alignments against D. melanogaster nAChR subunit sequences revealed earthworm sequences containing all loops relevant for acetylcholine/imidacloprid binding (such sequences were termed "full-length").
nAChR/AChBP Interaction Loop Residue Expression.The lists of full-length nAChR α, non-α, and AChBP subunits were used to determine the expression of specific amino acid residues in interaction loops critical to imidacloprid binding (i.e., residues experimentally associated with differential binding, EC 50 or maximal response values).As some earthworm populations are genetically heterologous and some species present higher allelic variation, the transcriptome assembly produces a different number of contig sequences for any given gene or gene family (including nAChR α, non-α, and AChBP subunits).This is an issue when comparing the expression between species from genetically homogeneous (e.g., Amynthas gracilis) with heterologous populations (e.g., Lumbricus rubellus).This is because sequencing reads map to a single contig (uniquely mapping reads) and the expression of a given gene will be split across a greater number of distinct contigs in a species with a higher genetic variation.To allow for species comparisons, we clustered all full-length contigs into either nAChR α, non-α, and AChBP subunits and summed the total normalized expression for each subunit type in each species across each tissue, giving the total expression by tissue for the three subunit types.This approach also allows for comparisons in the expression proportions of specific residue identities (in specific tissues) between species across the subunit types (calculated as percentage of total tissue expression).
Comparative Binding "Affinity" Scores for nAChR and AChBP Subunit Residues.The expression proportions for each residue (associated with binding) in each species was used to create comparative "affinity" and "sensitivity" scores (see the SI Results and Discussion for rules followed at each binding site, full sequences in SI Data S1−S3, and detailed residue quantifications in SI Data S4).Briefly, within a species, the percentage expression (in nerve cord) of residues experimentally associated with higher/lower imidacloprid binding were divided to give a higher/lower binding fraction at each site.For example, if a critical position is represented by a range of nonpolar, polar, and acidic residues across all fulllength nAChR α subunits, existing information is used to categorize the relationships of the various residues with imidacloprid binding.That is, they are associated with a higher binding or lower binding or, in cases where no confident assignment is possible, were considered to have neutral roles.In this example, if mutational experiments reveal that the replacement of nonpolar or polar residues with acidic residues reduces imidaclopr i db i n d i n ga n dv i c ev e r s a (indicated by altered imidacloprid concentrations activating the receptor to half its maximal level), then nonpolar/polar and acidic residues are grouped into higher and lower binding categories, respectively, and their relative expression percentages were used to produce the higher/lower binding fraction for that position.Combining residue identities possessing similar relationships to binding allows for the simultaneous consideration of multiple residue identities at the same position.The exact details of various decisions (and the justifying literature) for each residue across all subunit types are included in SI Data S4.Then, for the earthworm− earthworm comparisons, the high/low fractions for each residue in A. gracilis were then divided by the mean high/ low fractions for the two lumbricid species.Likewise, the fractions for D. melanogaster were divided by the A. gracilis values (in the D. melanogaster−A.gracilis comparison) and by the mean high/low fractions for the two lumbricid species (in the D. melanogaster−lumbricid comparison).These calculations gave a comparative affinity score for each site of interest, where a score <1 = relatively greater proportion of high-affinity residues in denominator species and >1 = relatively greater proportion of high-affinity residues in numerator species.At sites where there is no obvious higher/lower affinity relationship between residues but the expression of a particular residue is thought to enhance the imidacloprid binding (e.g., structural analysis suggests tyrosine at a site promotes imidacloprid ligand binding), the affinity score is produced directly by dividing the nerve cord residue expression percentages (of tyrosine in this case).
Earthworm nAChRs are assumed to be composed of a 3:2 ratio of α/non-α subunits (see the Supporting Information for discussion and justification of this assumption).Such pentamers are thought to contain three LBDs, two at the α/ non-α interfaces and one at the single α/α interface. 39This Environmental Science & Technology means residues within specificl o o p sm a yn o te q u a l l y contribute to ligand binding.Specifically, while binding-critical residues in loops A, B, and C in α subunits are predicted to contribute to all three LBDs, residues in loops G, D, E, and F in α subunits contribute to the single α/α interface (one-third of all LBDs), whereas loops D, E, and F in non-α subunits contribute to both α/non-α interfaces (two-thirds of all LBDs).To account for these differences, the affinity scores were adjusted.Initially, the reciprocal was calculated (1/affinity score) for scores predicting greater proportions of high-affinity residues in any denominator group.Then, the scores for loops G, D, E, and F (if a score was present) in α subunits were adjusted, ([affinity score − 1] × 0.333), and loops D, E, and F in non-α subunits, ([affinity score − 1] × 0.666).No adjustment was made to scores for loops A, B, and C in nAChR α subunits and all AChBP loops on the basis that they contributed to all LBDs.The subsequent residue scores for interacting loops were combined as follows: Loops associated with multiple residue scores but with all scores indicating higher affinity for either the denominator or numerator species, [(sum residue affinity scores) + 1].When residue scores for a single loop gave a combination of both higher denominator and numerator species scores, the following steps were taken: first, sum scores for numerator species and denominator species separately.The larger summed value reveals the overall affinity direction for the whole loop (i.e., toward higher affinity in the numerator or denominator species).Then, [(larger value − smaller value) + 1].For loops associated with only a single score: ([residue affinity score] + 1).The log 2 values of the resulting loop scores were calculated and represented the final loop affinity score.
Comparative Imidacloprid Sensitivity Scores for nAChR and AChBP Subunit Residues.For nAChR subunits, the affinity scores (prior to the log 2 conversion) for all loops in nAChR α and non-α subunits were combined.The loop affinity scores suggesting a higher affinity in the numerator species were summed, as were the loop affinity scores predicting a higher affinity in the denominator species.The category represented by the larger value was taken as the overall sensitivity direction (i.e., higher sensitivity in the numerator or denominator species).The smaller value was then subtracted from the larger value and log 2 converted to give the overall comparative nAChR sensitivity score.The process was completed for AChBP loops but with the exception that the category (either numerator or denominator species) represented by the smaller value was taken as the overall sensitivity direction (on the basis that high binding to an off-target protein reduces sensitivity).The expression of AChBPs critical to the sensitivity score were validated using reverse transcription polymerase chain reaction (RT-PCR) (see the SI Methods).

■ RESULTS AND DISCUSSION
Comparative Sensitivity of Earthworm Species.Imidacloprid measurements confirmed the exposure concentrations within an expected range, although indicating up to 35% biodegradation over 28 days (see Figure S1).GLM analysis confirmed that the measured concentrations were significantly dependent on nominal values (F = 352, P < 0.001).There was no significant effect of temperature (F = 0.01, p > 0.05) or any nominal value temperature interaction (F = 0.23, p > 0.05) on the measured concentrations.This indicates no temperature effects on imidacloprid fate (see Figure S1) that could explain any observed DSS.
LC 50 s for imidacloprid in the five species varied by a factor of >14.4 from the most sensitive A. gracilis to the least sensitive L. rubellus and A. caliginosa (Figure 1A).The DSS determined using EC 50 s for reproduction was broadly consistent with the LC 50 values, indicating DSS in the order A. gracilis > D. octaedra > A. caliginosa > E. fetida > L. rubellus (Figure 1A).Comparatively, A. gracilis was 31.7-foldmore sensitive than L. rubellus.Temperature can have a significant effect on toxicity, 40,41 including on imidacloprid in short-term bioassays (although decreasing with exposure time). 42In this study, however, the attribution of species difference to exposure temperature is not supported.For example, the two species tested at 20 °C(A.gracilis and E. fetida) were found to show a high DSS, being the most sensitive and second least sensitive species.Further, the order and 3-fold difference found here for sensitivity for E. fetida and D. octaedra tested at 20 and 13 °C, respectively, is consistent (indeed in a magnitude slightly less) with the 5-fold difference found at the same test temperature by Kreutzweiser et al. 43 Hence, exposure temperature clearly cannot explain DSS.We, therefore, looked for a physiological cause.
Toxicokinetics.Imidacloprid in tissues of exposed L. rubellus showed tissue/soil bioaccumulation factors of 1.6−2.7 based on nominal concentrations and 2.76−7.51based on the final measured soil concentrations.Bioaccumulation factors were lower at the higher exposure concentrations (Figure S2).When exposed at 0.37 mg/kg, the tissue concentrations varied significantly (ANOVA F = 15,34, p < 0.001; Tukey's test, P < 0.05) between species (4-fold between the lowest, A. gracilis, and the highest, D. octaedra)( Figure 1B).The ranking of imidacloprid accumulation and magnitude of species difference did not correlate with the sensitivity.Indeed, the most sensitive species, A. gracilis, showed the lowest tissue imidacloprid followed by the least sensitive, L. rubellus (Figure 2A,B).There was a pattern for higher tissue concentrations in smaller (e.g., D. octaedra) rather than larger (e.g., L. rubellus, A. gracilis) Environmental Science & Technology earthworms (Figure 1B).This may be explained by the sizerelated surface area to body mass ratio uptake limitation in the larger species.From this data, we conclude that the DSS cannot be explained by imidacloprid bioaccumulation factors alone.
To investigate the attribution of DSS to TK rates, we conducted 14 C-imidacloprid uptake and elimination experiments with the most sensitive species A. gracilis, the least sensitive L. rubellus, and the standard test species E. fetida (second least sensitive).The accumulation rates (k 1 ) determined from one compartment TK model fits varied 2.8fold (Figure 2A).The predicted elimination rates were low in all species (k 2 =6× 10 −5 − 4.5 × 10 −4 h −1 , estimated half-lives >1000 days all cases), resulting in near linear uptake and slow loss (Figure 2A,B).Slow earthworm imidacloprid elimination rates contrast with those for the teleost fish Oncorhynchus mykissa half-life of 60−70 h 44 and amphipod Gammarus pulexa half-life of 82 h. 45Earthworms (E.fetida) have previously shown to have moderate to slow elimination of a range of lipophilic organic compoundshalf-lives of 2−28 days 46,47 suggesting that earthworms possess a low biotransformation capacity for multiple organic chemicals, including imidacloprid.
The pattern of comparatively high tissue accumulation of imidacloprid by E. fetida and low accumulation by L. rubellus as detected by LC−MS (Figure 1B) was confirmed by radiochemical analysis; however, the kinetic study places A. gracilis closer to E. fetida than L. rubellus (Figure 2A,B).The difference in TK traits indicated by the radiolabel and unlabeled LC−MS measurements may provide indications of possible imidacloprid metabolism in species.In L. rubellus, uptake is constrained by the low assimilation rate, as confirmed by tissue measurements.In A. gracilis, the assimilation rate is higher but measurements indicate low parent compound concentrations.As the 14 C-imidacloprid radiolabel is placed in the stable (pyridyl-2,6) ring, this may suggest some degree of metabolic biotransformation by A. gracilis with the radiolabel retained as a radiolabeled metabolite.Such metabolism might be expected to reduce sensitivity, which is demonstrably not the case for A. gracilis.For E. fetida, a high assimilation rate leads to high tissue concentrations, but this is not associated with a high sensitivity.The tissue localization of 14 Cimidacloprid indicates higher concentrations in the putative target nerve cord than in the body wall or seminal vesicle.Nerve cord concentrations are higher in A. gracilis, but concentrations in the lumbricids also differ (Figure 2C), which is not consistent with their comparable sensitivities.From these results, we conclude that neither TK dynamics or distributions to critical nerve tissues can wholly explain imidacloprid DSS.This conclusion is in agreement with those of Van Den Berg et al. 48and Dalhof et al., 49 respectively, for the organophosphate chlorpyrifos and pyrethroid cypermethrin.Both of these studies further recommended a greater consideration of TD traits, which was done here through receptor analysis.
Toxicodynamics through Receptor Analysis.We investigated the potential contribution of divergent nAChRs and AChBPs to DSS.This requires consideration of how to convert the nAChR/AChBPs subunit gene sequence and expression differences into sensitivity predictions.The nAChR consists of α and non-α subunits combined into heteropentamer and α subunit homopentamer conformations 17 (Figure 3A,B), while AChBPs are homopentamers formed from α-like subunits (Figure 3C). 50In nAChR hetropentamers, imidacloprid is thought to bind at α/non-α subunit interfaces, as well as the single α/α interface present in 3:2 α/

Environmental Science & Technology
non-α ratio heteropentamers (Figure 3A). 39In heteropentamers, each ligand binding domain (LBD) is composed of several interacting loops: A, B, C, and G provided by the αsubunit and D, E, and F present on the non-α subunit (Figure 3D).In nAChR or AChBP homopentamers, all interacting loops are provided by the adjacent α-subunits (Figure 3D).As non-α subunits lack the necessary loop C sequences, 18 they cannot form homopentamers and only create binding sites when adjacent to α subunits with necessary loop sequences.−20 Characterizing the expression of residues across all nAChR subunits is needed to predict functionality and, by extension, DSS.
Reciprocal best BLAST identified 10, 11, and 10 AChBP subunit sequences for L. rubellus, E. fetida, and A. gracilis, (after manual collation to remove allelic variants).AChBPs have been previously described only in molluscs, 50,53 polychaetes, 54 and spiders 55−57 and have not to date been reported in insects.AChBPs are homopentamers formed from α-like subunits (Figure 3C) that share high amino acid homology with the extracellular domains of nAChRs. 53Each AChBP possesses five sites (Figure 3C) capable of binding acetylcholine 58 and imidacloprid 59 but no functional transmembrane domain. 53,57he functional role of AChBPs is still unknown, 57 but they are hypothesized to function away from the synapse, perhaps in nonsynaptic cholinergic communication. 50Our analyses indicate that earthworms have an expanded (∼10 AChBP subunit gene) repertoire compared to those of molluscs (1− 2 53,60 ) and spiders (<5 57 ).Like nAChRs, AChBP subunit expression is nerve cord biased (SI and Figure S3).AChBPs are highly expressed (∼21× greater than nAChR subunits, SI and Figure S3), with the expression dominated by a few very highly expressed orthologues.For these highly expressed AChBPs, we validated the expression using PCR and found high agreement with our RNASeq data (SI and Figure S4).Further, the striking correlation between the total AChBP subunit nerve cord expression and 14 C-imidacloprid uptake in each species nerve cord supports AChBP−imidacloprid interaction (Figure 4).Thus, highly expressed AChBP orthologues, in addition to nAChR targets, may play important roles in mediating imidacloprid toxicity in earthworms.
Electrophysiological investigations of mutant nAChRs exposed to imidacloprid reveal the influence of critical loop residues on both the receptor peak current (I max ) and the imidacloprid concentration necessary to activate the receptor to the half maximal level (EC 50 value).In theory, changes to either property could influence neonicotinoid action. 17owever, for imidacloprid, a high positive correlation has been found between neurotoxic activity and both its nAChR binding affinity and EC 50 value, 61,62 but no such relationship has been found for I max . 62Therefore, while residues altering I max were characterized (Figures S6−S8 and SI Data S4), our predictions focused on residues shown either to change EC 50 or be critical for imidacloprid binding. 18(Figure S5 and SI Data S4).Residues at critical sites with uncertain properties were assumed neutral for imidacloprid binding.

Environmental Science & Technology
Using full-length contigs, the expression of critical loop residues was determined (see the SI Methods and Results and Discussion approach and detailed quantitative loop analyses).Our analysis covered the α (SI and Figures S6 and S7) and non-α subunits (SI and Figure S8) and all relevant critical residues for each of the less sensitive earthworms L. rubellus and E. fetida (collectively termed the lumbricids), the sensitive earthworm A. gracilis, and the known highly sensitive insect D. melanogaster 63 using the published genome and gene expression data for this species.Given their structural homology to nAChRs, 18,57 we also consider AChBPs as potential off-target binding sites (SI and Figures S9 and S10).For AChBPs, our focus was limited to residues affecting EC 50 , as residues affecting I max have no obvious relevance to AChBPs as they are not functioning ion channels.The expression characterization of both nAChRs and AChBPs reveals variations in the expression of critical loop residues between differentially sensitive species.To validate the extent of expression stability at critical residues across multiple samples, we compared expression patterns between the L. rubellus nerve cord and closely related neural ganglia, finding an effectively identical expression across nAChR and AChBP subunits (Figure S11).
The relative expression of high and low imidacloprid binding residues was used to quantitatively predict a logarithmic affinity score for each interacting loop.A comparison of Drosophila to lumbricids (Figure 5A) and A. gracilis (Figure 5B) suggests, with a single exception, that Drosophila expresses a greater proportion of high-binding residues in nAChR subunit loops, most strikingly for loops C and E in α and nonα subunits (Figure 5A ,B).In contrast, the lumbricid to A. gracilis comparison shows only relatively minor expression differences in high-and low-binding residues across all loops (Figure 5C).Further, even the small expression differences seen are inconsistent and an affinity score in one loop is largely nullified by differences in another (e.g., loops D and E in non-α subunits).Unlike for the nAChRs, there are striking differences in the expression of key loop residues in AChBPs between the lumbricids and A. gracilis (Figure 5C).The lumbricids express a greater proportion of AChBPs with high-binding residues, with loops F and C presenting the most dramatic differences (Figure 5C).This means that, while A. gracilis may express the highest level of AChBP in the nerve cord, it is predicted to have a weaker binding affinity for imidacloprid.
For nAChRs, it is reasonable to assume that a higher imidacloprid affinity equates to a higher sensitivity, 61,62 meaning that the calculated affinity score (created by summing the comparative expression proportions of nAChR loop residues with known effects on imidacloprid binding) equates to a sensitivity score.Calculated scores suggest that differences in nAChR residues between the three earthworm species and Drosophila (Figure 6A,B) explain DSS between the two taxa.Between the lumbricids and A. gracilis,n A C h Rr e s i d u e differences affecting affinity and, thereby, sensitivity appear incapable of accounting for the >30× difference in imidacloprid sensitivity (SI and Figure 6).Unlike nAChR, AChBP do not play a critical role in synaptic transmission, instead likely being expressed away from this critical target site with only an ancillary role in synaptic function.It is, therefore, reasonable to assume that a higher AChBP−imidacloprid affinity may equate to a lower sensitivity, as a high affinity to AChBPs will mitigate toxicity via off-target binding, which will in turn reduce nAChR−imidacloprid to reduce nerve  Environmental Science & Technology excitation and resulting neurotoxicity.Therefore, for AChBPs, the inverse of the summed affinity score represents a sensitivity score and predicts a high sensitivity of A. gracilis relative to that of the lumbricids (Figure 6C), consistent with our observations of actual species sensitivity to imidacloprid (Figure 1A).
Taken overall, the data suggests that nonsynaptic AChBPs play a critical role in determining earthworm sensitivity to imidacloprid.The high expression of AChBPs and the correlation between AChBP expression and the relative 14Cimidacloprid uptake suggest that these proteins dominate the pool of potential imidacloprid targets in earthworms (Figure S3 and Figure 4).We hypothesize that possessing a pool of imidacloprid binding receptors associated with nonsynaptic functionality protects the critical synaptic nAChR target.If correct, this AChBP pool may, in part, contribute to the lower sensitivity of earthworms relative to insects that lack AChBPs.However, it is important to note that we predict the classic nAChR target also contributes to the differential sensitivity between the taxa (Figure 6), as the comparative affinity scores across the nAChR subunit loops suggest Drosophila nAChRs bind imidacloprid with a greater affinity than their earthworm equivalents (Figure 5).In contrast, our predictions suggest that differences in the classic nAChR target do not underlie the differential sensitivity between earthworms (Figure 6), as the earthworm species are predicted to possess nAChR repertoires with effectively identical imidacloprid binding affinities (Figure 5).However, unlike the nAChRs, the differential sensitivity between the earthworm species can be explained by the predicted differential binding capacities of their highly expressed AChBP pools.Specifically, the AChBPs in A. gracilis are predicted to have a notably weaker binding affinity for imidacloprid than that of the relatively insensitive lumbricids (Figure 5), meaning that sensitivity is increased (Figure 6)a s imidacloprid in nerve tissues of A. gracilis can be more readily given up to interact with the primary nAChR target to cause neurotoxic effects.
The potential for stoichiometric binding of imidacloprid to AChBPs to mitigate primary receptor binding, and as a result toxicity, identified for earthworms could also be important mechanisms for determining DSS in other nontarget taxa.Arachnids also possess AChBPs. 55,57Comparative studies have shown that the oribatid mites Oppia nitens has a low sensitivity to imidacloprid compared to insect species, 14 and neonicotinoids are also known to have a low efficacy against spider mites. 64Further, an analogous off-target stoichiometric binding mechanism is proposed to explain resistance development in malarial mosquitos via increased expression of another nontarget protein, the sensory appendage protein SAP2. 65ence, off-target binding may be a common but, until now, unrecognized driver of species DSS of relevance for de novo sensitivity predictions.
Past approaches for predicting sensitivity have focused on the separate roles of TK traits (e.g., internal biotransformation rate) 66,67 and TD-receptor traits (e.g., presence of known target orthologues). 10,68Further, LaLone et al. 9 have shown that not only just orthologue presence but also the amino acid composition at key motifs can influence DSS.These authors proposed the SeqaPASS approach, which uses the orthologue presence and motif and amino acid sequence information to predict sensitivity.While it is reasonable to suppose that species possessing known target orthologues that interact strongly with a chemical ligand may be sensitive, our study reveals the need for extended analyses, specifically the presence and expression of off-target binding as a potential mechanism for mitigating primary target site interactions.For many cases, truly accurate DSS predictions will require the development of tools capable of comparing and modeling complex, multisubunit, and target (and nontarget) receptor repertoires that are expressed from large and divergent gene families like nAChR subunits.In this paper, we see clear agreement between the observed sensitivity and toxicogenomically generated predictions.This reveals that our innovative but pragmatic approach, exploiting the expression patterns of critical residues in complex receptors, can predict the sensitivity of critical nontarget taxon to environmental toxicants.
Discussions of earthworm toxicity testing, imidacloprid uptake and distribution, earthworm genotyping, animal dissection and RNA isolation, RT-PCR validation of AChBP expression, assumptions about extrapolating residue function from insects/vertebrates to earthworms, nAChR α subunit and nAChR non-α subunit analyses, summary of nAChR and AChBP interaction loop expression differences, AChBP subunit analysis, and off-target binding could mitigate toxicity, table of primer sequences used, and figures of measured versus intended nominal concentrations of imidacloprid, measured imidacloprid in the tissues of L. rubellus exposed for 28 days at different soil concentrations, expression patterns of nAChR and AChBP subunits, PCR validation of the nerve cord expression of AChBPs, alignments of N-terminus LBD sequences from Drosophila melanogaster nAChR subunits with earthworm nAChR and AChBP subunits, relative expression proportions of residues in loops A, B, and C of nAChR α subunits, relative expression proportions of residues in loops G, D, E, and F of nAChR α subunits, relative expression proportions of residues in nAChR non-α subunits, relative expression proportions of residues in AChBP subunits, and comparison of expression proportions for residues critical to differential imidacloprid binding in the neural ganglion and nerve cord of L. rubellus (PDF) Dataset of all L. rubellus α nAChR subunits found using reciprocal best BLAST (PDF) Dataset of all E. fetida α nAChR subunits found using reciprocal best BLAST (PDF) Dataset of all A. gracilis α nAChR subunits found using reciprocal best BLAST (PDF) Dataset of summary of all nAChR subunits and nAChBP loop residue effects on imidacloprid I max and EC 50 values (XLSX)

Figure 1 .
Figure 1.Sensitivity and tissue accumulation of imidacloprid in five earthworm species.(A) Toxicity expressed as the 28 day LC 50 and EC 50 -reproduction (both with 95% CIs) for imidacloprid (milligrams per kilogram of soil) for five earthworm species derived from Probit (LC 50 ) and logistic (EC 50 ) model fits.(B) Tissue concentrations of imidacloprid in four replicates each with two individuals (when available) of five earthworm species exposed at 0.37 mg/kg soil; letter indicated difference between treatments (Tukey's Test p < 0.05).

Figure 2 .
Figure 2. 14 C-Imidacloprid toxicokinetic and tissue localization in two insensitive lumbricid and one sensitive megascolecid earthworm.(A) Uptake and elimination rate constants (k 1 , k 2 ± SE), imidacloprid half-life in tissues (T 1/2 ), and bioaccumulation factors (BAF) from the one compartment TK model fits for three earthworm species.(B) Uptake and elimination of 14 C-labeled imidacloprid for three earthworm species exposed in soil at a concentration of 0.0615 mg/kg soil for 336 h, followed by a 336 h elimination period in clean soil; the broken gray line indicates the start of the elimination period, the circles are measured data points, and the solid line is a one-compartment TK model fit.(C) Tissue localization of 14 C imidacloprid in the body wall (BW), nerve cord (NC), and seminal vesicle (SV) of four individuals from three earthworm species exposed at 0.0615 mg/kg soil for 168 h.

Figure 3 .
Figure 3. nAChRs and AChBPs in various pentamer conformations and interaction loops contributing to ligand binding domains (LBDs).(A) Schematic of nAChR heterompentamer with a 3:2 ratio of α/ non-α subunits (box shows the region expanded in D), with three LBDs shown as orange circles.(B) Schematic of nAChR α subunit homopentamer with five LBDs.(C) Schematic of AChBP homopentamer formed of AChBP α-like subunits.(D) Enhanced representation of boxed region in A, showing LBDs that form between loops/residues located at α/non-α and α/α subunit interfaces.At α/ non-α LBDs, the α-subunit provides loops A, B, and C and the non-α subunit provides D, E, and F. In α−α LBDs, loops D, E, and F are provided by the equivalent regions of α-subunits (with a further loop, G, potentially playing a role at these LBDs).The LBDs in homopentamers, whether nAChR or AChBP, occur in a manner similar to the α−α interface shown.

Figure 4 .
Figure 4. Relative levels of 14 C-imidacloprid uptake and AChBP expression in earthworm nerve cords.

Figure 5 .
Figure 5. Comparative imidacloprid logarithmic affinity-score predictions for nAChR and AChBP subunits based on the expression of variant residues that modulate imidacloprid binding.(A) D. melanogaster versus lumbricid earthworms (L.rubellus and E. fetida) predictive comparison across loops in nAChR subunit ligand binding domains (LBDs).(B) D. melanogaster versus A. gracilis predictive comparison across nAChR subunit LBD loops.(C) Lumbricid earthworms versus A. gracilis predictive comparison across nAChR and AChBP LBD loops.

Figure 6 .
Figure 6.Comparative imidacloprid sensitivity-score predictions based on the affinity scores for nAChR and AChBP subunits.(A) D. melanogaster versus lumbricid earthworms (L.rubellus and E. fetida) sensitivity score on the basis of nAChR subunit comparisons.(B) D. melanogaster versus A. gracilis sensitivity score on the basis of nAChR subunit comparisons.(C) Lumbricid earthworms versus A. gracilis on the basis of both nAChR and AChBP subunit comparisons.Sensitivity scores assuming highly expressed AChBPs represent a less critical target than nAChRs; the low-affinity AChBPs in A. gracilis are predicted to offer less protection for the nAChR, resulting in a predicted higher sensitivity for A. gracilis.