- Open Access
Up-regulation of apoptotic- and cell survival-related gene pathways following exposures of western corn rootworm to B. thuringiensis crystalline pesticidal proteins in transgenic maize roots
BMC Genomics volume 22, Article number: 639 (2021)
Resistance of pest insect species to insecticides, including B. thuringiensis (Bt) pesticidal proteins expressed by transgenic plants, is a threat to global food security. Despite the western corn rootworm, Diabrotica virgifera virgifera, being a major pest of maize and having populations showing increasing levels of resistance to hybrids expressing Bt pesticidal proteins, the cell mechanisms leading to mortality are not fully understood.
Twenty unique RNA-seq libraries from the Bt susceptible D. v. virgifera inbred line Ped12, representing all growth stages and a range of different adult and larval exposures, were assembled into a reference transcriptome. Ten-day exposures of Ped12 larvae to transgenic Bt Cry3Bb1 and Gpp34/Tpp35Ab1 maize roots showed significant differential expression of 1055 and 1374 transcripts, respectively, compared to cohorts on non-Bt maize. Among these, 696 were differentially expressed in both Cry3Bb1 and Gpp34/Tpp35Ab1 maize exposures. Differentially-expressed transcripts encoded protein domains putatively involved in detoxification, metabolism, binding, and transport, were, in part, shared among transcripts that changed significantly following exposures to the entomopathogens Heterorhabditis bacteriophora and Metarhizium anisopliae. Differentially expressed transcripts in common between Bt and entomopathogen treatments encode proteins in general stress response pathways, including putative Bt binding receptors from the ATP binding cassette transporter superfamily. Putative caspases, pro- and anti-apoptotic factors, as well as endoplasmic reticulum (ER) stress-response factors were identified among transcripts uniquely up-regulated following exposure to either Bt protein.
Our study suggests that the up-regulation of genes involved in ER stress management and apoptotic progression may be important in determining cell fate following exposure of susceptible D. v. virgifera larvae to Bt maize roots. This study provides novel insights into insect response to Bt intoxication, and a possible framework for future investigations of resistance mechanisms.
The efficacy of pesticidal agents that control feeding damage on agriculturally-important crop plants become reduced following repeated exposures and selection for resistance within target arthropod pest populations [1,2,3]. A diversity of bacterial pore-forming pesticidal proteins have been described including B. thuringiensis (Bt) crystalline (Cry), toxin-10 like (Tpp) and AeGerolysin like pesticidal proteins (Gpp) . Transgenic maize hybrids that express pesticidal proteins are widely used by growers in the United States and several other countries worldwide . The number of insect species with documented resistance continues to increase [3, 6].
The western corn rootworm, Diabrotica virgifera virgifera (Insecta: Coleoptera: Chrysomelidae), causes extensive damage to cultivated maize, Zea mays, throughout the major production regions of North America and Europe [7, 8]. This univoltine species diapauses over the winter months as eggs in the soil. At high population densities, maize root feeding by larvae which hatch during early summer can reduce yields through both physiological and mechanical damage [9, 10]. Adult beetles emerge and feed mainly on maize silk and pollen. Larval damage was historically controlled by soil insecticides applied at planting . In some regions, adults are sprayed aerially with insecticides to reduce egg-laying and hence larval populations the following year . After its first detection in 2009 , a high proportion of D. v. virgifera larvae in field populations in the United States now show resistance to Cry3Bb1  along with cross-resistance to the structurally similar mCry3A  and eCry3.1Ab proteins expressed by maize hybrids [16, 17]. Resistance to transgenic Cry34/35Ab1 (Gpp34/Tpp35Ab1 according to new nomenclature by Crickmore et al. (2021))  maize is also documented in D. v. virgifera field populations, but these phenotypes show no cross-resistance to Cry3 proteins [18,19,20]. Analogous lack of Gpp34/Tpp35Ab1cross-resistance with mCry3Aa was shown in laboratory selected D. v. virgifera . This resistance occurred to transgenic hybrids that express a “low-dose” of Bt proteins  and adversely impacts crop production [23, 24].
Ingested pesticidal proteins cause disruption of gut cell integrity in susceptible insects, leading to lethargic behavior, cessation of feeding, and death . Midgut cells of susceptible D. v. virgifera larvae fed diet containing Cry3Aa1 and Gpp34/Tpp35Ab1 swell and shed microvilli and other cell debris into the gut lumen . A proposed mode of action involves sequential binding of Bt pesticidal proteins to membrane-bound protein receptors on the apical side of midgut epithelial cell (enterocyte) membranes [27,28,29]. Within this model, receptor binding precedes the formation of an ion pore channel that generates an osmotic imbalance due to an influx of extracellular Ca2+ , leakage of gut contents, and eventual death of the insect. This model further proposes that ingested monomeric Cry toxins interact with the midgut apical membrane-bound receptor protein, cadherin, causing a conformational change in the pesticidal protein that leaves it liable to cleavage by gut proteases. Cleavage in turn facilitates subsequent oligomerization into a pre-pore structure  which inserts into the enterocyte membrane following interactions with cell surface aminopeptidase N  or alkaline phosphatase . ATP binding cassette (ABC) transporters are also implicated as receptors for Cry proteins [34,35,36], and hypothesized to facilitate Cry protein oligomerization and membrane pore insertion [37, 38]. Additional gut proteins have been identified as putative Bt protein receptors in insects [39,40,41,42], including the glycosyl hydrolase, α–amylase, from Tenebrio molitor , but their roles in Bt intoxication remains unknown.
An alternative model for Cry protein mode of action proposes a Mg2+-dependent G protein-mediated cell signalling pathway. In this model, binding of Cry proteins to cadherin directly triggers an intracellular adenylate cyclase signalling cascade that activate protein kinase A (PKA) and cell death (apoptosis) [44, 45]. This model hypothesizes a mechanism independent of other receptors or requirments for pore formation .
Mechanisms of resistance often implicate structural or functional changes in a midgut receptor protein that putatively disrupts the Bt mode of action, but are mainly from studies on Bt resistant species of Lepidoptera. This includes alteration of cadherin receptor protein structure by transposon-mediated insertional knockout or point mutations, that result in reduced Cry1A binding [47,48,49,50]. An amino acid change in a transmembrane domain of tetraspanin 1 is associated with Helicoverpa armigera resistance to Cry1Ac . Cry1A resistance is also associated with reduced expression of one or more aminopeptidase N paralog in Spodoptera exigua , Trichoplusia ni  and Ostrinia nubilalis , and alkaline phosphatase in Heliothis virescens . The ABC transporter subfamily C member, abcc2, is linked or associated with lepidopteran resistance to Cry1Ac in H. virescens , Plutella xylostella , Bombyx mori , and T. ni , Cry2Ab2 in Pectinophora possypiella , and Cry1Fa in O. nubilalis  and Spodoptera frugiperda . An ABC subfamily B member in in the coleopteran Chrysomela tremula is linked to Cry3Aa resistance, and is capable of mediating cell disruption via ectopic expression in S. frugiperda Sf9 cells . An abcb gene is also located in proximity to a single QTL for Cry3Bb1 resistance in D. v. virgifera .
Alternatively, enhanced repair of damaged midgut cells in response to Cry protein-mediated damage contributes to Cry1Ac resistance in H. virescens [64, 65], suggesting that stress-induced cell regeneration or degradation mechanisms are involved in physiological responses . Conversely, transcripts encoding proteins involved in apoptosis (programmed cell death) are significantly up-regulated in Manduca sexta cells by Cry1Ac , and in the midgut of S. exigua by exposure to the Bt vegetative insecticidal protein (Vip) 3A protein . Apoptosis was induced in Choristoneura fumiferana cells by a mechanism involving the mitogen activated protein kinase (MAPK) protein p38 following Cry1Ac exposure , and RNAi-mediated knockdown of MAPK p38 in Chilo suppressalis led to increased larval susceptibility to Cry1Ca .
Despite these insights, Bt mode(s) of action, mechanisms of cellular intoxication, and intracellular responses are not fully understood. In D. v. virgifera, mCry3A binding to midgut receptors is reduced among larvae selected for mCry3A resistance . Kinetic data demonstrate that Cry3Bb binds strongly to specific domains of the cadherin protein and enhances toxicity in D. v. virgifera  and other Coleoptera . However, cadherin is not considered a receptor in vivo since RNAi-based transcript knockdown does not alter D. v. virgifera susceptibility to Cry3Aa or Gpp34/Tpp35Ab1 . Estimates of differential gene expression shows no significant induction of cadherin in susceptible compared to resistant larvae fed Cry3Bb1 , eCry3.1Ab , or between susceptible larvae exposed and not exposed to Cry3Bb1  or Gpp34/Tpp35Ab1 . Although a suite of ABC transporters and aminopeptidase N transcripts are differentially-regulated in constitutive or induced fashions between Bt resistant and susceptible D. v. virgifera larvae to Cry3Bb1  and eCry3.1Ab , paralogs are also differentially-regulated in susceptible larvae exposed to Cry3Bb1 and Gpp34/Tpp35Ab1 [74, 76]. Furthermore, these transcriptome-wide comparisons have implicated a relatively large number of differentially expressed transcripts including those encoding proteins in general stress response pathways (e.g. cytochrome P450 monooxygenases, esterases, oxidases, and peroxidases) and those with transporter function .
A better understanding of how Bt intoxication affects gene expression among susceptible arthropods may reveal points at which mechanistic disruption could lead to resistance. To this end, we developed an inbred strain of Bt susceptible D. v. virgifera, Ped12, and used it for assembly of a comprehensive reference transcriptome, which was then applied to estimate transcript quantity differences following exposures to Cry3Bb1, Gpp34/Tpp35Ab1, and non-Bt maize within a common genetic background. Furthermore, transcripts differently expressed by Ped12 larvae following exposures to one or both Cry3Bb and Gpp34/Tpp35Ab1 maize, and exposures to an entomopathogenic fungus and nematode, were associated with generalized stress and immune response pathways. A filtered set of differentially expressed transcripts not shared with entomopathogens encoding candidate Bt receptor proteins, metabolic and detoxification enzymes, and proteins putatively determining cell fate (pro-survival or -death) were focused upon. This study contributes to an understanding of mechanisms potentially involved in determining cell fate (death or survival) which may inform future research into processes involved in Bt intoxication or resistance mechanisms.
Samples, treatments, and collections
Samples of Ped12 D. v. virgifera were collected from all developmental stages and different exposure conditions to create a reference transcriptome assembly (C1 to C20; Table 1). Among replicate treatments used in the downstream analysis of differential gene expression (T1 to T8; Table 2), Ped12 2nd instars were recovered from transgenic Cry3Bb1 (VT3; T8) and Gpp34/Tpp35Ab1 (Hx; T5) maize treatments after 48 h exposure. Approximately 1/3 of larvae were dissected from inside Cry3Bb1-expressing roots (T8), whereas none were found burrowing inside the roots of Gpp34/Tpp35Ab1 hybrid maize (T5). All D. v. virgifera larvae in the control non-Bt maize Corteva Pioneer hybrid 38B85 treatment (T7) were feeding within roots. Larvae exposed for 48 h to M. anisopliae (T3) and H. bacteriophora (T6) were lethargic, but none were moribund.
Complementary DNA libraries, sequencing and data processing
A total of 769.6 million reads were obtained after trimming, and nearly 600 million (78%) remained PE reads (Table 3a). Among all reads, 249.7 million (97.4 paired and 55.0 singletons) were produced from the normalized Pooled library (Supplementary Table S1). 21.7 ± 8.5 (mean + SE; range: 8.3 to 34.2) million reads were generated from among replicates of the non-normalized RNA-seq libraries. All resulting trimmed reads were used in the construction of the D. v. virgifera reference transcriptome, and the non-normalized libraries were used to estimate read counts as a proxy for predicting differentially expressed transcripts (Fig. 1). All raw Illumina read data were submitted to GenBank Sequence Read Archive (SRA) database under accessions ERR2791371 to ERR2791395 (Table 2; BioProject PRJEB28633; BioSample SAMEA4896309).
De novo reference transcriptome assembly and annotation
The D. v. virgifera reference transcriptome was assembled from all trimmed reads from the normalized cDNA library conditions (C1 to C20: Table 1) and non-normalized RNA-seq treatments (T1 to T8: Table 2), that used 82.5 million in silico normalized sequences (Table 3a). The multiple k-mer assembly approach was used (k = 61, 71, 81, 91) to generate a total of 228,885 redundant transcripts (Fig. 1a). After reduction in sequence redundancy and implementing a size cutoff ≥200 bp, a total of 116,070 transcripts comprised the final D. v. virgifera reference transcriptome (Fig. 1a; Table 3a), and were submitted to the NCBI Transcript Sequence Assembly (TSA) database (accession ERZ1775117.1).
Completeness was estimate by the presence of Core Eukaryotic Genes Mapping Approach (CEGMA) genes and BUSCO scores. Among the 248 CEGs, 98% were present in D. v. virgifera reference transcriptome (n = 243; mean 2.16 transcripts per CEG; 5 CEGs present as partial transcripts). A BUSCO score of 91.7% was obtained: 978 complete (821 complete single copy and 157 complete duplicated), 35 fragmented, and 53 missing BUSCOs were predicted.
Protein coding sequences were predicted within 56,656 of the 116,070 non-redundant D. v. virgifera transcripts, of which 23,728 received PFAM annotations (Table 3a). Among the 116,070 transcripts, 410 (0.35%), 107 (0.09%), and 65 (0.06%) were predicted to be putative Wolbachia, Heterorhabditis and Metarhizium transcripts, respectively (remaining data not shown). BLASTx searches resulted in matches to protein models from D. melanogaster (n = 17,322), T. castaneum (n = 24,543), D. ponderosae (n = 26,965), and A. glabripennis (n = 25,340), as well as in the SwissProt database (n = 20,961) (Table 3b; Supplementary Fig. S1). These results predicted that 7568 distinct D. v. virgifera transcripts to be “putative full-length” (proportional coverage lengths = 1.0) and 14,574 were “near complete (proportional coverage lengths ≥ 0.8 and < 1.0; Table 3b).
Estimates of quantitative differences in transcript expression
A total of 520.2 million reads (86.3%) were aligned to transcripts within the D. v. virgifera reference transcriptome (65.0 ± 21.5 million across treatments; 21.7 ± 8.5 million across replicates within treatments). Reads with multiple alignments, and those for which the mate-pair was aligned to a different target were discarded. Approximately half (52.1%) mapped properly (range 0.4036 and 0.5906; Supplementary Table S2). From these alignments, estimates of significant differences in read count (proxy for gene expression) for each transcript were generated from among replicates between control maize (Cn) relative to exposure treatments Cry3Bb1 (VT3; T8; Supplementary Table S3) and Gpp34/Tpp35Ab1 (Hx; T5; Supplementary Table S4), as well as entomopathogens H. bacteriophora (Hb; T6; Supplementary Table S5) and M. anisopliae (Ma; T3; Supplementary Table S6). Data from the comparison between T7 (control maize; Cn) and T8 (Cry3Bb1 maize; VT3) treatment groups, produced a DESeq2 adjusted read count for each transcript fitted to a dispersion around an empirically estimated mean (Supplementary Fig. 2A) and this was used to determine the significance of differences between treatments. Outliers within this dispersion were not fitted and not used in further DESeq2 analyses. MA-plots of estimated mean read counts normalized by size factor and transformed on a Log2(fold-change) scale showed that 2710 transcripts had differences in quantity that surpassed a Benjamini and Hochberg adjusted significance threshold (FDR) of ≥ 0.05. Among these transcripts, a greater number were up-regulated (n = 2503) than down-regulated (n = 207; Supplementary Fig. S2B; Supplementary Table S3). EdgeR estimated 1228 differentially-expressed transcripts for the same comparison between T7 and T8 with a greater number up-regulated (n = 1014) than down-regulated (n = 214). There was a strong correlation between Log2(fold-change) estimates between DESeq2 and EdgeR methods (R2 = 0.9806; Fig. 2a).
The comparison of read count data between replicate libraries from control maize (Cn; T7) versus transgenic Gpp35/Tpp35Ab1 maize exposed larvae (Hx; T5) similarly resulted in adjusted dispersions fitted to the mean (Supplementary Fig. S2C), and a final distribution of Log2(fold-change) in an MA plot (Supplementary Fig. S2D). DESeq2 output predicted 2389 transcripts that surpassed an adjusted significance threshold (FDR ≤ 0.05; Supplementary Table S4). EdgeR analysis of the same data predicted 1690 differentially expressed transcripts. Overall estimates of Log2(fold-change) were highly correlated between EdgeR and DESeq2 (R2 = 0.9758; Fig. 2a). Significant levels of differential expression (FDR ≤ 0.05) were predicted between Cn and Hb treatments for 1818 and 2376 transcripts using DESeq2 and EdgeR, respectively (Supplementary Table S5). Similarly, significant levels of differential expression (FDR ≤ 0.05) were predicted between Cn and Ma treatments for 1583 and 1684 transcripts by DESeq2 and EdgeR, respectively (Supplementary Table S6). No other comparisons were conducted or evaluated.
Differential expression following Cry3Bb1 and Gpp34/Tpp35Ab1 exposure
Our pipeline defined differential expression occurring only among transcripts having an adjusted P-value (FDR) ≤ 0.05 in both DESeq2 and EdgeR analyses (Fig. 1b). Furthermore, any differentially expressed transcripts after entomopathogen exposures and pesticidal protein treatments were subtracted to account for putative general stress response genes. Specifically, 1562 transcripts showed significant differential expression among D. v. virgifera exposed to H. bacteriophora compared to control larvae in both DESeq2 and EdgeR analyses, of which 1269 and 293 were up- and down-regulated, respectively (Table 4; Supplementary Table S5). Analogously, 1199 differentially expressed transcripts were predicted between M. anisopliae exposure and control treatments, with 639 and 563 transcripts up- and down-regulated, respectively (Table 4; Supplementary Table S6). Among the most prevalent PFAM annotations assigned to predicted differentially-regulated transcripts in both Hb and Ma treatments were those with cytochrome P450, transporter, and protease and protease inhibitor domains (Supplementary Table S7; Supplementary Table S8).
Comparison between T7 (control maize; Cn) and T8 (Cry3Bb1 maize; VT3) treatments predicted a total of 1064 differentially expressed transcripts by both DESeq2 and EdgeR (Table 4; Supplementary Table S3). Subsequent BLASTx results showed 4 and 5 of these transcripts putatively originated from Wolbachia sp., and the entomopathogens Hb and Ma, respectively. Among the remaining 1055 transcripts 942 and 113 were up- and down-regulated, respectively (Table 4; Fig. 2a). Comparison showed that 257 endogenous transcripts differentially expressed between Cry3Bb1 and controls were also differentially expressed in Hb and/or Ma treatments compared to controls (Fig. 2b). GO enrichment analyses of these shared transcripts predicted secretory vesicle (category CC), C-N bonding forming ligase activity (MF), and purine compound biosynthesis process (BP) among the most significantly overrepresented (Supplementary Fig. S3A). After removal of these 257 transcripts shared with Hb and Ma treatments 798 endogenous transcripts were unique to the Cry3Bb1 response. Of the 775 unique PFAM domain annotations assigned to 609 of these 798 differentially expressed transcripts (76.3%), cathepsin inhibitor (Inhibitor_I29), papain (Peptidase_C1) and trypsin proteases functional domains were most prevalent (Table 5). Transcripts encoding alkaline phosphatase, aminopeptidase, or cadherin domains were not among those that were differentially expressed. Annotations did indicate that three transcripts up-regulated in Cry3Bb1 treatments encoded putative ABC transporter-like proteins from subfamilies ABCG (n = 1) and ABCC (n = 2) (Table 6), while a single transcript was annotated with a tetraspanin domain (DIAVI021979). Transcripts putatively encoding six apoptosis-related proteins, including two caspases, one IAP, and the BAX-domain containing protein BI-1 were up-regulated by Cry3Bb1 exposed larvae (Table 7). The most significantly enriched GO terms assigned at level 2 to differentially expressed transcripts in the Cry3Bb1 treatment were in extracellular space and microbody (category CC), coenzyme binding and channel regulator activities (MF), and small molecule catabolism and drug metabolism processes (BP) (Fig. 3a).
Comparisons between control (Cn; T7) and Gpp35/Tpp35Ab1 maize (Hx; T5) treatments predicted significant differences for 1387 transcripts by both DESeq2 and EdgeR (Table 4; Supplementary Table S4). Among these transcripts, 13 showed top BLASTx hits to Wolbachia sp. (n = 7), H. bacteriophora (n = 1), or other microbial protein database sources (n = 5), and were removed from the dataset. This filtered set of 1374 D. v. virgifera transcripts (Table 4; Fig. 2a) contained 295 that were also differentially expressed in one or both of the Hb and/or Ma treatments (Fig. 2b). GO enrichment analyses predicted the most significant over-representation was for secretory vesicle (category CC), coenzyme binding activity (MF), and organophosphate biosynthesis process (BP) (Supplementary Fig. S3B). Following removal of these 295 transcripts shared with Hb and Ma responses, a total of 1079 were retained in the dataset of endogenous transcripts uniquely responding to Gpp34/Tpp35Ab1 (Fig. 2b). A total of 1066 PFAM domain annotations were assigned to 837 of these 1079 differentially-expressed transcripts (77.6%), with sugar transporter (Sugar_tr), cytochrome P450 (p450), and lectin C-type domain (Lectin_C) numerically greatest (Table 5). No alkaline phosphatase, aminopeptidase, or cadherin domains were annotated among differentially expressed transcripts. PFAM domains for ABC transporter were assigned to four up-regulated transcripts, one assigned to each of the ABCB, C, G, and E subfamilies (Table 6). By comparison, a total of 6 and 5 transcripts encoding ABCC subfamily members were differentially regulated in H. bacteriophora and M. anisopliae treatments, respectively (Table 6), but none were predicted in common with those in the Cry3Bb1 or Gpp34/Tpp35Ab1 treatments. A set of nine apoptosis-related protein-encoding transcripts were up-regulated in Gpp34/Tpp35Ab1 exposed larvae, of which five were in common with those also up-regulated in Cry3Bb1 maize (Table 7). Transcripts uniquely up-regulated following exposure to Gpp34/Tpp35Ab1 maize included one putatively encoding a BI-1 ortholog, and a different IAP1 isoform (X1) than the IAP1 isoform X5 up-regulated in Cry3Bb1 treatments. Enrichment at GO level 2 showed greatest significance within in microbody and secretory vesicle (category CC), coenzyme binding activity (MF), and small molecule catabolism processes (BP) in the Gpp34/Tpp35Ab1 treatment (Fig. 3b).
A set of 696 differentially-expressed endogenous transcripts were shared between both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments (Supplementary Table S9; Fig. 4a). The Log2(fold-change) estimates for these transcripts were highly correlated between DESeq2 and EdgeR analyses (R2 ≥ 0.8524; Fig. 4a). Following removal of 133 transcripts that were also differentially expressed in Hb and/or Ma treatments, a set of 535 endogenous transcripts were shared and responding to both Cry3Bb1 and Gpp34/Tpp35Ab1 (Supplementary Table S9; Fig. 4b). Predicted PFAM domains showed that of putative cytochrome P450 (p450), cathepsin inhibitor (Inhibitor_I29), sugar transporter (Sugar_tr), papain cysteine protease (Peptidase_C1), and carbohydrate binding module (CBM_14) were the five most prevalent (Table 5). GO enrichment at level 2 for differentially expressed transcripts shared in pesticidal protein exposures was greatest within CC categories microbody, secretory vesicle, and extracellular space (Fig. 5). Additionally, enrichment was greatest for coenzyme binding activity (MF), and small molecule catabolism, monocarboxylic acid metabolism, and drug metabolism processes (BP) (Fig. 5). An ABCC subfamily member putatively orthologous to the T. castaneum ABCC-ST gene was significantly up-regulated in the Cry3Bb1 and Gpp34/Tpp35Ab1 treatment (Table 6). Five apoptosis-related proteins were up-regulated in both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments (Table 7). The transcript DIAVI004770, encoding a putative tetraspanin-like protein, was up-regulated in both Cry3Bb1 and Gpp34/Tpp35Ab1 exposure treatments.
Overall, these results show that after filtering, differentially-expressed transcripts in separate Cry3Bb1 and Gpp34/Tpp35Ab1 treatments encode proteins putatively most enriched for those localized in the extracellular space, secretory vesicles, and microbodies, and having coenzyme binding function and are involved in small molecule catabolism (Fig. 3). These same categories are also most significantly enriched among differentially expressed transcripts share in both treatments (Fig. 5). Of the putative Bt receptor proteins, only two were differentially expressed across Bt maize exposure larvae (Table 6), whereas five transcripts putatively encoding apoptosis and cell stress-related proteins were upregulated in both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments (Table 7).
Phylogenetics and structural annotations
BLASTp results showed that 7 D. v. virgifera transcripts (this study) and 9 D. v. virgifera gene models surpassed E-value thresholds against D. melanogaster caspases: DRONC, death regulator Nedd2-like caspase (FlyBase ID: FBgn0026404); DRED, Death related ced-3/Nedd2-like caspase (FBgn0020381), DAMM (FBgn0033659), STRICA, Ser/Thr-rich caspase (FBgn0033051), DECAY, Death executioner caspase related to Apopain/Yama (FBgn0028381), DCP-1, death caspase 1 (FBgn0010501), and DRICE, death related ICE-like caspase (FBgn001997; data not shown). A 197 amino acid consensus alignment was generated for a partial enzymatic domain region among seven D. melanogaster and putative D. v. virgifera caspases (Supplementary Fig. S4). Percent identity among aligned sequences ranged from 17.12 to 100.00, with catalytic histidine (H) and cysteine (C) residues 100% conserved. The LG + G + I model maximized the BIC score at 4496.352, and was implemented as the “Best Model” for subsequent phylogenetic reconstruction. The subsequent unrooted ML-based phylogeny had a G of 2.7553 and I of 6.77% that minimized at the log likelihood score of − 2224.28, resulting in a consensus tree of total branch length of 10.157 (Fig. 6). Two major clades comprised of effector and initiator caspases were supported by 65% of 1000 bootstrap pseudo-replicates of the data. Clade was defined based on phylogenetic position of D. melanogaster effector (DRICE, DCP-1 and DECAY) and initiator caspases (DRONC, DRED, STICA and DAMM), and predicted the effector DECAY and initiators STRICA or DAMM as nearest orthologs to up-regulated D. v. virgifera caspases DIAVI027204 and DIAVI022989, respectively (Fig. 6). Overall, each transcript within the assembled transcriptome showed a one-to-one phylogenetic correspondence with a Dvir_v2 protein model, with the exception of a lack of homologous transcripts for XP_028140498.1 and XP_028140499.1.
Local BLASTn and BLASTp searches of D. v. virgifera gene models using corresponding sequences from putative apoptosis- or cell stress-related proteins encoded by the differentially expressed transcripts, IAP, BI-1, and LFG, resulted in identification of nearest homologs. Queries of the CDD with proteins encoded by DIAVI007715 and DIAVI011972 differentially expressed following D. v. virgifera exposure to Cry3Bb1 or Gpp34/Tpp35Ab1, respectively, allowed classification of both as IAP family 1 (IAP1) members based on presence of two baculoviral inhibition of apoptosis protein repeat (BIR) domains (Supplementary Fig. S5A). BLASTn results indicated that DIAVI011972 and DIAVI007715 were derived from Dvir_v2.0 LOC114335598, and represent splice variants XM_028285849.1 (isoform X1; protein translation XP_028141650.1; E-value = 0.0, %ID = 100.0) and XM_028285854.1 (isoform X5; protein XP_028141655.1; E-value = 0.0, %ID = 100.0), respectively. By comparison, FPAM and CDD results predicted that the non-differentially expressed transcript DIAVI011430 (homolog of XP_028283750.1; Dvir_v2.0 LOC114333756 E-value = 0.0, %ID = 100.0), encoded three BIR domains which is structurally similar to Drosophila DIAP2 orthologs (Supplementary Fig. S5B). Alignments of D. v. virgifera IAP1 and IAP2 BIR domains with D. melanogaster orthologs DIAP1 and DIAP2 showed ≥ 22.39% and ≥ 26.47% identities, respectively (Supplementary Fig. S5C). Analogous results were obtained for BAX inhibitor-like, BI-1 (Supplementary Fig. S6A) and Lifeguard 4 (LFG4) (Supplementary Fig. S6B), and SERP2-like proteins (Supplementary Fig. S7).
These results showed that in Cry3Bb1 and Gpp34/Tpp34Ab1 treatments the two upregulated caspases belonged in two separate clades corresponding to initiator and effector caspases (Fig. 6), that function at different stages to propagate apoptosis-related proteolytic cascades. Transcripts upregulated in both treatments also included a IAP2 class protein that correspondingly repressed progression of caspase cascade, and the remaining upregulated transcripts BI-1, LGF4 and SERP2 show homology to related peptides with putative function in repression of cell death and promotion of cell survival.
Transcriptome responses to Bt pesticidal protein exposures
The modes by which Bt intoxication evokes changes that elicit cell death or recovery, tissue damage or repair, and organismal mortality or survival are not yet fully understood [46, 66, 78]. Several studies demonstrate a role for gut membrane-bound protein receptors in mediating Bt intoxication, but uncertainties remain regarding molecular mechanisms underlying subsequent physiological and cellular changes. Specifically, pore-formation/osmotic imbalance, signal transduction, or hybrid models suggest different causes . Studies comparing transcriptome-wide expression differenes between susceptible insects exposed or unexposed to Bt proteins can provide valuable mechanistic insights, but such studies tend to identify hundreds or thousand of differentially-expressed genes [74,75,76, 79,80,81,82]. This phenomenon is also observed as part of insect responses to chemical insecticide exposures  and pathogen infections . Such outcomes not only impose challenges for interpreting results, but also in formulating hypotheses to guide future reseach.
In the current study we aimed to reduce the genetic component of variation in responses of D. v. virgifera by using an inbred strain Ped12. Additionally, we removed any transcripts from further consideration that were not identified by both DESeq2 and EdgeR packages, resulting in substantial reductions of 63.1 and 49.1% in candidate transcripts for Cry3Bb1 and Gpp34/Tpp35Ab1 exposures, respectively. Of the remaining transcripts, we performed in silico subtraction to remove transcripts that also were differentially expressed following exposures to entomopathogens (Hb or Ma), which was hypothesized to eliminate transcripts in shared general stress response pathways. These subtraction measures reduced the number of significantly differentially expressed transcripts in susceptible D. v. virgifera larvae exposed to transgenic Cry3Bb1 and Gpp34/Tpp35Ab1 maize to 798 and 1079, respectively (≥ 24.4% reduction) (Table 4; Fig. 2). Similarly, among the 696 differentially expressed transcripts shared between the two Bt maize treatments, an additional 133 (19.1%) were also differentially regulated in Hb and Ma treatments, leaving 563 after removal (Fig.4b). Among the filtered transcripts, are those encoding candidate Bt receptor proteins, metabolic and detoxification enzymes, and proteins putatively involved in cell fate (pro-survival or -death) were over represented, and these were focused upon in our further investigations.
Up-regulation of putative B. thuringiensis receptors
Transcripts predicted to encode previously identified receptor proteins were differentially-expressed among susceptible D. v. virgifera exposed to Cry3Bb1 and/or Gpp34/Tpp35Ab1 compared to controls. For example, tetraspanin encoding transcripts were up-regulated in response to Cry3Bb1 (transcript DIAVI021979) and Gpp34/Tpp35Ab1 (transcripts DIAVI004770 and DIAVI021979; Table 6). A non-synonymous change in the transmembrane domain in the H. armigera tetraspanin gene, HaTSPAN1, was linked to Cry1Ac resistance in strain AY2 . The HaTSPAN1 transcript levels are increased 2.7-fold in AY2, but did not alter Cry1Ac binding. The means by which structural changes and up-regulation of tetraspanin interrupts the Bt mode of action in H. armigera AY2 , or role of tetraspanin-like transcripts in Cry3Bb1 and Gpp34/Tpp35Ab1 responses by susceptible D. v. virgifera remains unknown. Alternate midgut receptors in resistant insects have been proposed to sequester Bt proteins, thereby reducing binding to membrane-bound proteins functionally involved in pore formation [85, 86], such as ABC transporters, cadherin, aminopeptidase N or alkaline phosphatases. Regardless, a potential protein sequestering role of tetraspanin remains to be investigated.
The involvement of ABC transporters in Bt intoxication has been demonstrated through linkage or association with resistance among species of Lepidoptera to subfamily members abcc2 [34, 56, 57, 60], abca2 , and abcg [88, 89]. In the current study of Bt susceptible D. v. virgifera larvae, five transcripts encoding putative ABC transporters were significantly up-regulated in response to Bt maize proteins (Table 6), and agree with results from a prior study for a eCry3.1Ab resistant strain of this species . These contrast with other studies of D. v. virgifera where ABC transporter expression was not significantly different following exposures of susceptible larvae to Cry3Bb1 or Gpp34/Tpp35Ab1 pesticidal proteins , or where transcripts were not detectable in gut tissues . Differences in ABC transporter transcription may be dependent upon environmental conditions or genetic background of strains being compared, although using an inbred strain may have minimized effects of the latter in our study.
Evidence from other systems indicate that ABC transporters are involved in pro-survival stress response mechanisms among mammals [90, 91]. Specifically, ABCC subfamily members function with glutathione S-transferases (GSTs) and UDP-galactosyl transferases (UGTs) to enhance drug and carcinogen efflux in cellular maintenance of homeostasis in human and mouse [92,93,94]. Our predicted up-regulation of transcripts encoding GST, UGT, and ABCC transporter domains in D. v. virgifera following Cry3Bb1 and Gpp34/Tpp35Ab1 exposure (Table 5), and significant enrichment for GO MF category drug metabolism and BP category drug metabolism (Fig. 3; Fig. 5), could suggest increased cellular transport may be involved in responses to Bt intoxication. ABC transporters also have other cellular functions within insects including immune responses , consistent with up-regulation of unique ABCC members in responses to entomopathogens (Table 6). Although mutations in a specific ABC transporter may inhibit pesticidal protein pore formation in resistant insects, proteins in the same superfamily may mediate other cellular responses following pore formation in susceptible insects. It may be possible that modulation of ABC transporter expression impacts cellular efflux capacities in response to increased solute influx through pores, but this hypothesis requires additional investigation.
Modulation of metabolic and detoxification pathways
Our analyses showed that the most significantly enriched GO terms among differentially expressed transcripts following exposure of D. v. virgifera to Cry3Bb1 maize encompassed extracellular space and microbody (category CC), binding and transport functions (MF), and small molecule catabolism and drug metabolism (BP) (Fig. 3a). Transcripts significantly differentially-expressed in the Gpp34/Tpp35Ab1 maize treatment were most enriched for terms microbody, ER membrane, and extracellular space (category CC), coenzyme binding (MF), and small molecule catabolism (BP) (Fig. 3b). In conjunction with prior studies [75, 76], our results suggest an increase in metabolic processes may be a common response among D. v. virgifera to Bt intoxication. In other organisms, metabolic rates increase during instances of cellular repair and survival [96, 97], which might be conserved and potentially explanatory of the observed increase in metabolic pathway gene expression in susceptible D. v. virgifera responses to Cry protein intoxication. Also, enrichment of GO terms for extracellular space shared between Cry3Bb1 and Gpp34/Tpp35Ab1 treatments here and in eCry3.1Ab resistant D. v. virgifera , might suggest a role of secreted factors in intoxication responses.
Even though no GO enrichment analyses were conducted in a prior investigation of Cry3Bb1 maize exposure among susceptible D. v. virgifera larvae , these authors described significant differential expression of transcripts putatively encoding proteins involved in xenobiotic stress responses and detoxification (cytochrome P450 monooxygenase, esterase, oxidase, and peroxidase) functions, and transporter activities. Our study similarly predicted PFAM domains for sugar transporter and detoxification enzymes (cytochrome P450, carboxyesterase, glutathione S transferase, and UDP glycosyltransferases) are prevalent among transcripts differentially regulated in Cry3Bb1 exposed larvae (Table 5). Our data also show that sugar transporter and cytochrome P450 domains encoded by transcripts differentially-regulated in the Gpp34/Tpp35Ab1 treatment are shared between the Cry3Bb1 treatment. Coincidence of these functions between this prior study  and our current study suggest a role for these proteins in cellular intoxication response. Cytochrome P450s are involved in a large breadth of insect cellular functions , including regulation of insect ecdysone and juvenile hormone pathways , cuticle formation , and xenobiotics detoxification . Uncoupling of P450 oxygenation reactions results in production of reactive oxygen species (ROS), hydrogen peroxide and superoxide . Cellular homeostasis can become disrupted during times of cell stress when excess ROS is produced due to high P450 activity. ROS can trigger apoptosis , or act as second messengers that modulate other cellular processes . Stress responses triggered by ROS during high metabolic states are intimately tied to increased ABC transporter activities , suggesting a possible basis for our predicted up-regulation of ABC transporters in Cry protein and entomopathogen responses (Table 6). This also highlights the value of identifying transcripts putatively involved in general (i.e. not Bt-specific) cellular stress. Moreover, we hypothesize that increases in metabolic, transport, and detoxification pathways following Bt exposure of susceptible insects may be connected with the increased energy demands of cellular repair or death/survival processes. Despite the tantalizing connections, further research is required to demonstrate the roles of these pathway components for cellular or organismal survival.
Up-regulation of cell survival pathways
The current study predicts that two D. v. virgifera transcripts putatively encoding orthologs of D. melanogaster effector caspase DECAY and initiator caspase STRICA or DAMM (Fig. 6) are up-regulated following feeding on Cry3Bb1 and Gpp34/Tpp35Ab1 maize roots (Table 7). Although DRONC and effectors DRICE and DCP1 are the main caspases involved in apoptosis, DECAY and DAMM may represent redundancies or have yet unknown functions . Regardless, the up-regulation of caspase-encoding transcripts in D. v. virgifera following Bt intoxication could suggest an apoptotic response. This is in partial agreement with prior studies showing significant up-regulation of transcripts encoding caspases in Manduca sexta cells following Cry1Ac exposure  and in midgut tissues of S. exigua following a sublethal exposure to the Bt Vip3A protein .
Because caspase activation by the apoptosome is a critical step in determining cell fate, this process is tightly regulated. In mammals, caspase translation occurs as inactive pro-peptides that undergo autocatalysis in response to pro-apoptotic stimuli . However, evidence suggests D. melanogaster caspases are translated in active forms but are suppressed following reversible binding by inhibitor of apoptosis protein (IAP) family members . Although several mechanisms function to regulate pro-apoptotic signals, only the evolutionarily conserved IAPs inhibit caspase function through direct binding [108, 109]. Paralogs IAP1 or IAP2 are each sufficient to inhibit cell death in lepidopteran cells , suggesting their functional conservation. IAPs are up-regulated in response to cell stress conditions , as was observed here following exposures of susceptible D. v. virgifera to Cry3Bb1 and Gpp34/Tpp35Ab1 (Table 7). Interestingly, the up-regulation of caspase- and IAP-encoding transcripts was concurrent, perhaps suggesting an intimate balance between pro- and anti-apoptotic signals within cells following Bt intoxication. For example, caspase-mediated apoptosis could be partially suppressed through the up-regulation of IAPs. Although Drosophila IAP1 transcription is regulated by factors including those in the Hippo pathway , the conservation of this regulatory framework across arthropods remains unknown. Future investigation of the basis and consequences of IAP up-regulation in D. v. virgifera may clarify the role of apoptosis in determining cell fate, and organismal survival following low-dose Cry intoxication.
Transcripts encoding anti-apoptotic proteins including structurally and functionally conserved transmembrane B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor 1 motif (TMBIM)-containing protein family members, the BAX inhibitor 1 (BI-1) and Lifeguard 4 (LFG4) , are significantly up-regulated in Cry3Bb1 and Gpp34/Tpp35Ab1 exposed susceptible D. v. virgifera larvae (Table 7). Similarly, the serine protease inhibitor, stress-associated ER protein 2 (SERP2), is also up-regulated. These factors function in cellular adaptation to stress in the mitochondrion, Golgi or endoplasmic reticulum (ER), thereby suppressing apoptosis. Specifically, Bcl-2 protein family members are critical for regulation of the intrinsic pathway of apoptosis in mammals . In this system, pro-apoptotic Bcl-2 proteins, including the BAX protein, can oligomerize under cell stress conditions to form mitochondrial outer membrane pores that release cytochrome c and other factors, which in turn cause caspase activation. In Drosophila, there is a single pro-apoptotic Bcl-2 protein, decbl, that is a functional BAX ortholog , but evidence suggests it has a limited role in triggering apoptosis . TMBIM proteins modulate stress through several intracellular mechanisms . For example, BI-1 does not inhibit BAX via direct protein-protein interaction [118, 119], but instead inhibits ROS production in the mitochondrion. In the ER, BI-1 remediates the unfolded protein response (UPR) and H+ antiporter activity, counteracting cytosolic acidification characteristic of apoptosis . This may be important because prolonged presence or high accumulation of misfolded proteins can lead to pro-apoptotic signaling, and the UPR can restore ER homeostasis . Therefore, BI-1 promotes cell survival pathways by suppressing factors that would otherwise promote apoptotic signaling . SERP also functions within the UPR by enhancing protein stability [123, 124]. Somewhat analogously, LFG4 remediates stress in the Golgi and ER , where the single ortholog in D. melanogaster interacts with the anti-apoptotic Bcl-2 protein, buffy, resulting in organismal survival via repression of pro-apoptotic decbl . BI-1 and LFG4 structure and function are remarkably conserved, where ectopic expression of viral orthologs can rescue knockdown phenotypes in mammalian cells , suggesting retention of ortholog function in D. v. virgifera. Although the up-regulation of transcripts encoding TMBIM members BI-1 and LFG4, as well as SERP4, in susceptible D. v. virgifera following Cry3Bb1 and Gpp34/Tpp35Ab1 maize exposure (Table 7) suggests activation of pathways to counteract cell death by apoptosis or related mechanisms, additional investigations are required to determine individual contributions and impacts on cellular and organismal outcomes.
A greater understanding of cellular responses among susceptible insects to pesticidal protein exposure can inform future research into mechanisms that lead to resistance. To date mutations in cell surface receptors have mainly been implicated in pesticidal protein resistance among insects, although some evidence suggests alteration of intracellular signaling or cell recovery mechanisms may have a role (see Introduction). Only a few prior studies demonstrated apoptotic pathways in Bt protein response among insects, specifically involving up-regulation of caspases [67, 68] or the MAPK p38 pathway [69, 70]. We provide evidence for an apoptotic response in susceptible D. v. virgifera larvae following exposure to transgenic maize roots expressing either Bt Cry3Bb1 or Gpp34/Tpp35Ab1. Interestingly, the patterns we observed in differential expression suggest counteracting pathways may simultaneously remediate cell stress and suppress apoptosis, possibly through modulating a balance between opposing pathways to determine cell fate. Because we used whole larvae, we cannot predict any possible tissue- or cell type-specific responses. In general, exposure level has a role in cell response to pore forming proteins suggesting a “high-dose” leads to death by oncosis (cell swelling and blebbing due to osmotic imbalance) as opposed to a “low-dose” that tends to trigger apoptosis . Therefore, the responses of D. v. virgifera to the “low-dose” pesticidal protein exposures characteristic of current transgenic Bt maize hybrids commercialized for their control may not be comparable to responses among species of Lepidoptera that feed on crop tissues that provide a “high-dose”. Regardless, this work provides a framework for understanding cellular responses to Bt pesticidal protein exposure in the most devastating maize pest in the United States, and suggests that mechanisms promoting and counteracting apoptosis may characterize these responses.
Samples, treatments, and collections
A Bt susceptible non-diapausing strain of D. v. virgifera  maintained at the United States Department of Agriculture, Agricultural Research Service, North Central Agricultural Research Laboratory (USDA-ARS, NCARL), Brookings, SD, USA was used to generate an inbred line, Ped12. Ped12 was initiated from a single mated pair, followed by inbreeding for 9 generations: single pair full-sib mating in the F1 to F5; (generations G1 to G5), followed by en masse mating among siblings within G6 and G7, then single pair full-sib mating in the subsequent generations (G7 to G9). Ped12 was maintained thereafter as a closed colony of approximately 1000 individuals for 8 generations prior to use in this study.
Ped12 individuals were sampled at different developmental stages or following different exposure conditions. Different developmental stages (conditions C1 through C12; Table 1) were sampled among Ped12 individuals reared using standard laboratory methods  at USDA-NCARL. For conditions C13 to C16 (Table 1) maize seedling mats were germinated from seeds of hybrids expressing modified mCry3A (Event MIR604; Syngenta, Basel, Switzerland;), Cry3Bb1 (Event Mon88017 by Bayer Crop Sciences/Monsanto Company, St. Louis, MO, USA in VT TriplePro, VT3 hybrid DCK61–69 from DeKalb Seed, DeKalb, IL, USA), Gpp34/Tpp35Ab1 (Event DAS-59122-7 by Corteva/Pioneer, Indianapolis, IN, USA in HerculexXtra hybrid 2 T789 from Mycogen Seeds, San Diego, CA, USA), or no Bt (hybrid 38B85, Corteva/Pioneer). In brief, recently laid Ped12 D. v. virgifera eggs were suspended in a 0.15% agar solution and dispensed into 15 by 10-cm plastic containers (708 ml; The Glad Products Company, Oakland, CA) at a rate of 500 eggs per container. Containers with eggs were filled with 20 ml of water and 150 ml of the soil mixture described previously. After 1 wk., 50 maize seeds were added to containers and covered with an additional 300 ml of soil mixture and 80 ml of water. Containers were held in a controlled environmental chamber (Powers Scientific Inc., Pipersville, PA) at constant 25 °C and a photoperiod of 14:10 (L:D) h, as described previously . Four weeks after infestation of eggs, larvae of each treatment were recovered from seedling mats, placed in tin foil packages, flash frozen in liquid nitrogen, and stored at − 80 °C. Seedling mats in treatments C13 to C15 received regular watering. Condition C16 simulating drought stress of hybrid 38B85 (Corteva/Pioneer) received no water (Table 1).
Additionally, Ped12 larvae were separately exposed to the entomopathogenic nematode Heterorhabditis bacteriophora (Hb) strain BU (Becker Underwood, Ames, IA, USA; C17 in Table 1), and the entomopathogenic fungus, Metarhizium anisopliae (Ma) strain F52 (provided by Stefan Jaronski, USDA-ARS; C18 in Table 1). In brief, 300 μl of a 510 Hb ml− 1 suspension was aliquoted onto filter paper in 10 cm Petri plates, and six 2nd and 3rd instars were added and exposed for 48 h. Analogously, 300 μl of 1.5 × 107 Ma conidia ml− 1 in 0.1% Tween 20 was aliquoted onto filter paper in 10 cm Petri plates, and six 2nd and 3rd instars exposed for 48 h (C18). Adults were exposed in a modified glass scintillation vial exposure assay . For this, 20 ml vials were coated with sublethal levels of the neonicotinoid insecticide, thiamethoxam (Poncho, Bayer Crop Sciences, Leverkusen, Germany; C19 in Table 1) or the adult attractant cucurbitacin (Sigma-Aldrich, St. Louis, MO, USA; C20 in Table 1). Ten Ped12 adults were exposed per vial for 24 h. In all cases, individuals were pooled as a single replicate per condition, flash frozen in liquid nitrogen, and stored at − 80 °C.
Ped12 individuals were also subjected to different treatments in triplicate and sampled for quantitative analyses of transcript expression (Table 2) and inclusion within the reference transcriptome assembly (Fig. 1). Pooled samples were collected from Ped12 eggs in diapause (treatment 1; T1), and 1st (T2) and 3rd instars (T7; Table 2) reared under standard laboratory conditions  at USDA-NCARL. Second instar Ped12 were exposed to maize seedling roots expressing Gpp34/Tpp35Ab1 (Herculex® XTRA Corteva/Pioneer Event DAS-59122-7 in hybrid 2 T789 from Mycogen Seed, San Diego, CA, USA; T5), Cry3Bb1 (VT TriplePro; Bayer Crop Sciences/Monsanto Company; T8), or the non-Bt hybrid 38B85 (Corteva/Pioneer; T7) as described above. Treatment 4 (T4; Table 2) consisted of larval midgut tissue dissected and pooled from approximately 50 3rd instars feeding on the non-Bt hybrid 38B85. Treatments 3 (T3) and 6 (T6) were from Ped12 larvae exposed to M. anisopliae and H. bacteriophora respectively, as described above. All pooled samples were flash frozen in liquid nitrogen separately within replicates for each treatment, and stored at − 80 °C.
Complementary DNA libraries, sequencing, and data processing
Total RNA was purified by replicate for each condition (C1 to C20 in Table 1), and triplicates within each treatment (T1 to T8 in Table 2). Each sample was ground in liquid nitrogen and then 10.0 mg of tissue added to 250 μl TriZol Reagent (Life Technologies, Grand Island, NY, USA). Purification was conducted using the Direct-zol RNA MiniPrep kit (Zymo Research, Irvine, CA, USA), which included a 15 min DNase I digestion performed according to manufacturer instructions. Total RNA extracts were quantified using Qubit dsDNA HS Assay Kits (Life Technologies) on a Qubit 2.0 Fluorimeter (Thermo-Fisher, Wilmington, DE), and quality determined by electrophoresis on a 10 cm 2% denaturing agarose gel in 1X MOPS buffer run at 70 V for 45 min.
Normalized cDNA was prepared from an equimolar pool of total RNA isolated across 20 different conditions (C1 to C20; Table 1) using the Trimmer-2 cDNA Normalization Kit according to manufacturer instructions (Evrogen, Moscow, Russia). Resulting cDNA was quantified on a Nanodrop2000 (Thermo-Fisher). Non-normalized cDNA was prepared separately from 1.0 μg of purified RNA for each of the three replicates per treatment (T1 to T8 in Table 2) using the SMARTer cDNA Synthesis Kit (Life Technologies) according to manufacturer instructions. Long-range PCR of SMARTer cDNA products was carried out using Advantage Taq Polymerase according to manufacturer instructions (Life Technologies, Carlsbad, CA), which included 18 amplification cycles according to manufacturer instructions on a Tetrad2 thermocycler (BioRad, Hercules, CA, USA) with 10 min extension at 65 °C. All amplified cDNAs were quantified using dsDNA HS Assay Kits (Life Technologies) on a Qubit 2.0 Fluorometer (Thermo-Fisher).
The normalized and non-normalized cDNA samples were sent to the Laboratoire de Sequencage, Institute de Genomique, France, where uniquely indexed RNA-seq libraries were prepared using TruSeq RNA library Prep Kit V.1 (Illumina Inc., San Diego, CA, USA). Sequencing was performed on an Illumina HiSeq2000 platform (Illumina, Inc.) in 100 bp paired-end (PE) mode. Non-normalized libraries (8 conditions × 3 biological replicates per condition; n = 24; Table 2) were run on three different lanes of the same flow cell, with one replicate from each condition per lane. The normalized pooled library was run on a single lane of a separate Illumina HiSeq2000 flow cell. All raw reads were trimmed to remove Illumina adaptor and low-quality sequences with Phred33 bases quality scores (q) < 20 using Trimmomatic v 0.32 . Long poly(A/T) stretches were removed using SEQCLEAN (Dana-Farber Cancer Institute, Boston, MA, USA; https://sourceforge.net/projects/seqclean/files/) with the -l option to retain trimmed reads ≥30 bp. Ambiguous nucleotide bases (‘N’) were removed using a custom in-house PERL script.
De novo reference transcriptome assembly and annotation
The D. v. virgifera reference transcriptome assembly pipeline used an iterative approach with in silico read normalization (Fig. 1a). Specifically, trimmed reads from the pooled normalized library (n = 1) and each replicate of all non-normalized RNA-seq libraries (n = 24) were combined into PE and single end (SE) read groups. In silico normalization was conducted separately for each using the script insilico_read_normalization.pl (available in TRINITY v2014-07-17 software package  with parameter --max_cov 30. This normalization step was performed to reduce the coverage of highly expressed transcripts, thereby improving the ability to assemble transcripts expressed at low levels. In silico normalized trimmed PE and SE reads were then merged across all libraries into a single fastq file, and used in a multiple k-mer assembly approach with VELVET 1.2.03  and OASES 0.2.06 . K-mer hashes were prepared for k = 61, 71, 81 and 91 with in silico normalized SE and PE reads applied as Short and Short Paired fastq input classes, respectively. After each single k-mer assembly, CD-HIT-EST  was used to remove shorter redundant transcripts when entirely covered by other transcripts with > 99% identity. Transcripts ≥ 100 bp from the four k-mer runs were then merged into a combined assembly using kmerge. CD-HIT-EST was once again used to remove the shorter redundant transcripts (same parameters as previously). To remove additional residual redundancy, we used a custom tool which merges two sequences using successive iterations of CAP3  and NRCL tools with a progressive reduction in the stringency parameters . The final threshold included shared ≥ 94% identity with ≥ 100 bp overlap and overhangs < 40 bp. The final set of transcripts were filtered to remove those < 200 bp to satisfy NCBI Transcriptome Shotgun Assembly (TSA) minimum requirements.
Completeness of the de novo D. v. virgifera reference transcriptome was evaluated by ortholog identification using the Core Eukaryotic Genes Mapping Approach (CEGMA) . Transcripts were mapped to 248 CEGs, and the set of 1066 universal single copy orthologs from Arthropoda obtained from OrthoDB v 9  as identified using BUSCO v 3  (E-value cutoff ≤ 0.05).
Assembled transcripts were annotated by comparing to the SWISSPROT database, and to predicted protein databases from gene models for insects (Tribolium castaneum v3.0, Dendroctonus ponderosae, A. glabripennis, Drosophila melanogaster r5.46), nematode (H. bacteriophora), fungus (M. anisopliae), and NCBI Wolbachia protein sequences. In all cases the BLASTx algorithm  was used for querying, and the results filtered for those with an E-value ≤ 10− 7 and high-scoring segment pairs (HSP) length ≥ 25 amino acids. Putative transcript origin of H. bacteriophora, M. anisopliae, or Wolbachia were predicted according to the lowest E-value obtained from the set of database query results. Transcripts were considered “full-length” if query match lengths covered the entire sequence of the best-hit protein. Twenty missing amino acids were allowed in HSPs at both ends of the subject if the query length was compatible with the complete the sequence. Among the remaining D. v. virgifera transcripts, a “near complete” bin was defined as those with a best HSP that covered ≥ 80% of the subject length. Putative peptide-coding region (CDS) and derived protein sequences were predicted among all assembled D. v. virgifera transcripts using FRAMEDP . Putative protein family domains were assigned by searches of the PFAM A database v. 27.0 [143, 144] using the program HMMSEARCH  with derived D. v. virgifera proteins as the queries.
Estimates of quantitative differences in transcript expression
Transcript abundances (gene expression levels) were estimated by comparing non-normalized read counts among treatments (T1 to T8; Table 2), which relied on an experimental design and analysis pipeline that accounted for variance among replicates within and between treatment (Fig. 1b). To accomplish this, trimmed non-normalized read data from 3 replicates within treatment (T1 to T8; Table 2; n = 24) were mapped separately to the de novo assembled D. v. virgifera reference transcriptome (above) using the BOWTIE2 aligner v2.1.0  with parameters --all (report all alignments), --end-to-end (entire read must align), and --sensible (0 mismatches allowed in a seed of 22 bp). Mapped reads were filtered using “high stringency” parameters (PE and SE reads were retained if mapped properly on only one transcript. In the case of PE reads, both left and right reads were required to map in opposite directions on the same transcript at a distance compatible with the expected mean size of the fragment). The resulting output files were parsed by an in-house script to count the number of reads that mapped to each transcript.
Significant differences in aligned read counts (gene expression) were predicted for the comparisons of treatment T7 (Cn maize) with treatments T3 (Ma), T6 (Hb), T5 (Hx) or T8 (VT3) using the R packages DESeq2 v1.6.3  and EdgeR v3.8.6 . Other possible comparisons were not conducted. DESeq2 fit the dispersion of read counts for each transcript to an empirical mean, and performed independent filtering (default alpha = 0.1) to remove transcripts with low read counts (baseMean), which are subject to greater uncertainty and influence the results of multiple testing. P-values were adjusted in both packages for a false discovery rate (FDR) among multiple comparisons using the Benjamini-Hochberg (BH) method . For each comparison with an FDR ≤ 0.05 were considered significant for DESeq2 and EdgeR, and the final set of differentially expressed genes was defined as those with an FDR ≤ 0.05 using both statistical methods.
Differential expression following Cry3Bb1 and Gpp34/Tpp35Ab1exposure
Linux awk and uniq commands were used to compile a dataset of transcripts differentially expressed in both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments compared to controls (Cn). The same methods were used to identify and filter transcripts that were differentially expressed in entomopathogen Hb and/or Ma treatments compared to controls, and that were also shared with those differentially expressed in one or more of the treatments 1) Cry3Bb1 vs Cn 2) Gpp34/Tpp35Ab1 vs Cn, and 3) both Cry3Bb1 and Gpp34/Tpp35Ab1 vs Cn. This filtering excluded transcripts not specific to Bt toxin Cry responses.
Putative cellular component (CC), molecular function (MF), and biological process (BP) categories were assigned to differentially expressed transcripts using the PFAM domains to retrieve corresponding gene ontologies (GOs). This was accomplished using the dcGOR 1.0.6 package . Subsequent enrichment analyses by dcGOR 1.0.6 at GO level 2 applied significance thresholds of E-values ≤ 1.0− 5 for CC and MF terms, and ≤ 1.0− 6 for BP. Corresponding E-values and total number of transcripts represented within each GO category at level 2 were plotted for the comparisons 1) Cry3Bb1 vs Cn, 2) Gpp34/Tpp35Ab1 vs Cn, and 3) both Cry3Bb1 and Gpp34/Tpp35Ab1 vs Cn. Differentially expressed transcripts with ABC transporter PFAM domains were assigned putative orthologs as described previously . Putative tetraspanin-like proteins translated from differentially expressed transcripts DIAVI004770 and DIAVI021979 were used as BLASTp queries to the NCBI nr protein database and flybase . This was repeated for H. armigera TSPAN1.
Phylogenetics and structural annotations
Protein sequences of Drosophila melanogaster caspases were downloaded from FlyBase  by gene name search: DRONC, DRED, DAMM, STRICA, DECAY, DCP-1, and DRICE, . The 28,061 RefSeq proteins derived from gene models in the annotated draft D. v. virgifera genome assembly GCA_003013835.2 Dvir_v2 (NCBI GenBank Accession PXJM00000000.2) were downloaded in fasta format, and loaded into a local BLAST database. This database was queried with protein sequences for D. melanogaster caspases using the BLASTp algorithm  (E-value cutoffs ≤ 10− 20). An analogous search was conducted against the 56,656 proteins predicted from all assembled transcripts. A multiple protein sequence alignment was generated among D. melanogaster and D. v. virgifera caspase catalytic domains using the Clustal W algorithm  within the MEGA8.0 alignment utility  (default parameters). The LG + G + I model of protein sequence evolution  maximized the BIC score when alignment gaps were ignored, and was subsequently used to construct a ML-based phylogeny with a consensus built from among 1000 bootstrap pseudo-replicates using MEGA8.0 .
Sequences from up-regulated transcripts encoding putative inhibitor of apoptosis proteins (IAPs; DIAVI011972 and DIAVI007715), B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor (BI) proteins BI-1 (DIAVI026079) and Lifeguard 4-like (LFG4; DIAVI029891), and stress-induced endoplasmic reticulum protein 2 (SERP2) were used to query the local database of 28,061 Dvir_v2 RefSeq proteins using the BLASTp algorithm  (E-value cutoffs ≤ 10− 60). In addition to PFAM and GO annotations (above), further structural annotation and position of functional domains and residues were predicted by query of the conserved domain database (CDD)  using default parameters. Multiple protein sequence alignments were also generated for D. v. virgifera IAP, and BI-1 and LFG4 proteins with orthologs from D. melanogaster, and other coleopteran species (Tribolium castaneum, Dendroctonus ponderosae, A. glabripennis and Leptinotarsa decemlineata) using the Clustal W algorithm. Classification of D. v. virgifera IAP protein families were based on that defined for D. melanogaster orthologs . Prediction of transmembrane region (TMR) for D. v. virgifera BI-1 and LFG4 applied a Hidden Markov Model using the default parameters of the application, TMHMM 2.0 .
Availability of data and materials
The datasets used and/or analysed during the current study are available in the National Center for Biotechnology Information (NCBI) under BioProject PRJEB28633 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB28633) containing Sequence Read Archive (SRA) runs ERX2800462 to ERX2800486 (https://www.ncbi.nlm.nih.gov/sra?linkname=bioproject_sra_all&from_uid=491890). Assembled transcripts are available from the NCBI Transcriptome Shotgun Assembly (TSA; accession ERZ1775117.1). Additional data generated or analyzed in this study are included in this published article and its supplementary information files.
B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor-1
AeGerolysin-like pesticidal protein
stress-associated ER protein 2
toxin-10-like pesticidal protein
Roush RT, McKenzie JA. Ecological genetics of insecticide and acaricide resistance. Annu Rev Entomol. 1987;32(1):361–80. https://doi.org/10.1146/annurev.en.32.010187.002045.
Mallet J. The evolution of insecticide resistance, have the insects won? Trends Ecol Evol. 1989;4(11):336–40. https://doi.org/10.1016/0169-5347(89)90088-8.
Tabashnik BE, Carrière Y. Surge in insect resistance to transgenic crops and prospects for sustainability. Nat Biotechnol. 2017;35(10):926–35. https://doi.org/10.1038/nbt.3974.
Crickmore N, Berry C, Panneerselvam S, Mishra R, Connor TR, Bonning BC. A structure-based nomenclature for bacillus thuringiensis and other bacteria-derived pesticidal proteins. J Invertebr Pathol. 2021. https://doi.org/10.1016/j.jip.2020.107438.
James, C. Global status of commercialized biotech/GM crops. International Service for the Acquisition of Agribiotech applications brief 53. NY, USA: Ithaca; 2017. https://www.isaaa.org/resources/publications/briefs/53/download/isaaa-brief-53-2017.pdf
Peralta C, Palma L. Is the insect world overcoming the efficacy of bacillus thuringiensis? Toxins. 2017;1:39.
Sappington TW, Siegfried BD, Guillemaud T. Coordinated Diabrotica genetics research, accelerating progress on an urgent insect pest problem. Am Entomol. 2006;52(2):90–7. https://doi.org/10.1093/ae/52.2.90.
Dillen K, Mitchell PD, Looy TV, Tollens E. The western corn rootworm, a new threat to European agriculture, opportunities for biotechnology? Pest Mant Sci. 2010;66(9):956–66.
Gray ME, Steffey KL. Corn rootworm (Coleoptera: Chrysomelidae) larval injury and root compensation of 12 hybrids: an assessment of the economic injury index. J Econ Entomol. 1998;91(3):723–40.
Urías-López MA, Meinke LJ. Influence of western corn rootworm (Coleoptera: Chrysomelidae) larval injury on yield of different types of maize. J Econ Entomol. 2001;94(1):106–11. https://doi.org/10.1603/0022-0493-94.1.106.
Gray ME, Sappington TW, Miller NJ, Moeser J, Bohn MO. Adaptation and invasiveness of western corn rootworm: intensifying research on a worsening pest. Annu Rev Entomol. 2009;54:303–21.
Meinke LJ, Siegfried BD, Wright RJ, Chandler LD. Adult susceptibility of Nebraska western corn rootworm (Coleoptera, Chrysomelidae) populations to selected insecticides. J Econ Entomol. 1998;91(3):594–600.
Gassmann AJ, Petzold-Maxwell JL, Keweshan RS, Dunbar MW. Field-evolved resistance to Bt maize by western corn rootworm. PLoS One. 2011;6(7):e22629.
Shrestha RB, Dunbar MW, French BW, Gassmann AJ. Effects of field history on resistance to Bt maize by western corn rootworm, Diabrotica virgifera virgifera LeConte (Coleoptera: Chrysomelidae). PLoS One. 2018;13(7):e0200156. https://doi.org/10.1371/journal.pone.0200156.
Gassmann AJ, Petzold-Maxwell JL, Clifton EH, Dunbar MW, Hoffmann AM, Ingber DA, et al. Field-evolved resistance by western corn rootworm to multiple bacillus thuringiensis toxins in transgenic maize. Proc Natl Acad Sci U S A. 2014;111(14):5141–6. https://doi.org/10.1073/pnas.1317179111.
Jakka SR, Shrestha RB, Gassmann AJ. Broad-spectrum resistance to bacillus thuringiensis toxins by western corn rootworm (Diabrotica virgifera virgifera). Sci Reports. 2016;6:27860.
Zukoff SN, Ostlie KR, Potter B, Meihls LN, Zukoff AL, French L, et al. Multiple assays indicate varying levels of cross resistance in Cry3Bb1-selected field populations of the western corn rootworm to mCry3A, eCry3.1Ab, and Cry34/35Ab1. J Econ Entomol. 2016;109(3):1387–98. https://doi.org/10.1093/jee/tow073.
Wangila DS, Gassmann AJ, Petzold-Maxwell JL, French BW, Meinke LJ. Susceptibility of Nebraska western corn rootworm populations (Coleoptera, Chrysomelidae) populations to Bt corn events. J Econ Entomol. 2015;108(2):742–51.
Gassmann AJ, Shrestha RB, Jakka SR, Dunbar MW, Clifton EH, Paolino AR, et al. Dound, J., St. Clair CR. Evidence of resistance to Cry34/35Ab1 corn by western corn rootworm (Coleoptera, Chrysomelidae), root injury in the field and larval survival in plant-based bioassays. J Econ Entomol. 2016;109(4):1872–80.
Ludwick DC, Meihl LN, Ostlie KR, Potter BD, French L, Hibbard BE. Minnesota field population of western corn rootworm (Coleoptera, Chrysomelidae) shows incomplete resistance to Cry34Ab1/Cry35Ab1 and Cry3Bb1. J Appl Entomol. 2017;141(1–2):28–40.
Zhao JZ, O'neal MA, Richtman NM, Thompson SD, Cowart MC, Nelson ME, et al. mCry3A-selected western corn rootworm (Coleoptera, Chrysomelidae) colony exhibits high resistance and has reduced binding of mCry3A to midgut tissue. J Econ Entomol. 2016;109(3):1369–77. https://doi.org/10.1093/jee/tow049.
Narva KE, Siegfried BD, Storer NP. Transgenic approaches to western corn rootworm control. In: Vilcinskas A, editor. Yellow biotechnology II. Springer, Berlin: Heidelberg; 2013. p. 135–62. https://doi.org/10.1007/10_2013_195.
Cullen EM, Gray ME, Gassmann AJ, Hibbard BE. Resistance to Bt corn by western corn rootworm (Coleoptera, Chrysomelidae) in the US Corn Belt. J Integr Pest Manage. 2013;4(3):D1–6.
Gassmann AJ, Shrestha RB, Kropf AL, St Clair CR, Brenizer BD. Field-evolved resistance by western corn rootworm to Cry34/35Ab1 and other bacillus thuringiensis traits in transgenic maize. Pest Manag Sci. 2020;76(1):268–76.
de Maagd RA, Bravo A, Crickmore N. How bacillus thuringiensis has evolved specific toxins to colonize the insect world. Trends Genet. 2001;17(4):193–9. https://doi.org/10.1016/S0168-9525(01)02237-5.
Bowling A, Pence H, Li H, Tan S, Evans S, Narva K. Histopathological effects of Bt and TcdA insecticidal proteins on the midgut epithelium of western corn rootworm larvae (Diabrotica virgifera virgifera). Toxins. 2017;9(5):156.
Pigott CR, Ellar DJ. Role of receptors in bacillus thuringiensis crystal toxin activity. Microbiol Mol Biol Rev. 2007;71(2):255–81.
Bravo A, Likitvivatanavong S, Gill SS, Soberón M. Bacillus thuringiensis, a story of a successful bioinsecticide. Insect Biochem Mol Biol. 2011;41(7):423–31. https://doi.org/10.1016/j.ibmb.2011.02.006.
Vachon V, Laprade R, Schwartz JL. Current models of the mode of action of bacillus thuringiensis insecticidal crystal proteins, a critical review. J Invertebr Pathol 2012;111(1):1–12, 1, DOI: https://doi.org/10.1016/j.jip.2012.05.001.
Potvin L, Laprade R, Schwartz JL. Cry1Ac, a bacillus thuringiensis toxin, triggers extracellular Ca2+ influx and Ca2+ release from intracellular stores in Cf1 cells. J Exp Biol. 1998;201(12):1851–8.
Gómez I, Sánchez J, Miranda R, Bravo A, Soberón M. Cadherin-like receptor binding facilitates proteolytic cleavage of helix α-1 in domain I and oligomer pre-pore formation of bacillus thuringiensis Cry1Ab toxin. FEBS Lett. 2002;513(2–3):242–6.
Bravo A, Gómez I, Conde J, Munoz-Garay C, Sánchez J, Miranda R, Zhuang M, Gill SS, Soberón, M. Oligomerization triggers binding of a Bacillus thuringiensis Cry1Ab pore-forming toxin to aminopeptidase N receptor leading to insertion into membrane microdomains. Biochimica et Biophysica Acta (BBA)-Biomembranes 2004;1667(1):38–46.
Arenas I, Bravo A, Soberón M, Gómez I. Role of alkaline phosphatase from Manduca sexta in the mechanism of action of bacillus thuringiensis Cry1Ab toxin. J Biol Chem. 2010;285(17):12497–503.
Gahan LJ, Pauchet Y, Vogel H, Heckel DG. An ABC transporter mutation is correlated with insect resistance to bacillus thuringiensis Cry1Ac toxin. PLoS Genet. 2010;6(12):e1001248. https://doi.org/10.1371/journal.pgen.1001248.
Heckel DG. Learning the ABCs of Bt, ABC transporters and insect resistance to bacillus thuringiensis provide clues to a crucial step in toxin mode of action. Pesticide Biochem Physiol. 2012;104(2):103–10. https://doi.org/10.1016/j.pestbp.2012.05.007.
Tanaka S, Miyamoto K, Noda H, Jurat-Fuentes JL, Yoshizawa Y, Endo H, et al. The ATP-binding cassette transporter subfamily C member 2 in Bombyx mori larvae is a functional receptor for cry toxins from bacillus thuringiensis. FEBS J. 2013;280(8):1782–94.
Ocelotl J, Sánchez J, Gómez I, Tabashnik BE, Bravo A, Soberón M. ABCC2 is associated with bacillus thuringiensis Cry1Ac toxin oligomerization and membrane insertion in diamondback moth. Sci Reports. 2017;7(1):2386.
Stevens T, Song S, Bruning JB, Choo A, Baxter SW. Expressing a moth abcc2 gene in transgenic Drosophila causes susceptibility to Bt Cry1Ac without requiring a cadherin-like protein receptor. Insect Biochem Mol Biol. 2017;80:61–70. https://doi.org/10.1016/j.ibmb.2016.11.008.
Fernandez-Luna MT, Lanz-Mendoza H, Gill SS, Bravo A, Soberon M, Miranda-Rios J. An alpha-amylase is a novel receptor for bacillus thuringiensis ssp. israelensis Cry4Ba and Cry11Aa toxins in the malaria vector mosquito Anopheles albimanus (Diptera, Culicidae). Environ Microbiol. 2010;12(3):746–57. https://doi.org/10.1111/j.1462-2920.2009.02117.x.
Zhang Q, Hua G, Bayyareddy K, Adang MJ. Analyses of α-amylase and α-glucosidase in the malaria vector mosquito, Anopheles gambiae, as receptors of Cry11Ba toxin of bacillus thuringiensis subsp. jegathesan. Insect Biochem Mol Biol. 2013;43(10):907–15.
Ochoa-Campuzano C, Real MD, Martinez-Ramirez AC, Bravo A, Rausell C. An Adam metalloprotease is a Cry3Aa bacillus thuringiensis toxin receptor. Biochem Biophys Res Commun. 2007;362(2):437–42. https://doi.org/10.1016/j.bbrc.2007.07.197.
Yamaguchi T, Bando H, Asano S. Identification of a bacillus thuringiensis Cry8Da toxin-binding glucosidase from the adult japanese beetle. Popillia japonica J Invertebr Pathol. 2013;113(2):123–8.
Bulushova NV, Zhuzhikov DP, Lyutikova LI, Kirillova NE, Zalunin IA, Chestukhina GG. Toxin-binding proteins isolated from yellow mealworm Tenebrio molitor and wax moth Galleria mellonella. Biochem. (Moscow) 2011;76(2):202–208.
Zhang X, Candas M, Griko NB, Rose-Young L, Bulla LA Jr. Cytotoxicity of bacillus thuringiensis Cry1Ab toxin depends on specific binding of the toxin to the cadherin receptor BT-R1 expressed in insect cells. Cell Death Differ. 2005;12(11):1407–16. https://doi.org/10.1038/sj.cdd.4401675.
Zhang X, Candas M, Griko NB, Taussig R, Bulla LA. A mechanism of cell death involving an adenylyl cyclase/PKA signaling pathway is induced by the Cry1Ab toxin of bacillus thuringiensis. Proc Natl Acad Sci U S A. 2006;103(26):9897–902. https://doi.org/10.1073/pnas.0604017103.
Soberón M, Gill SS, Bravo A. Signaling versus punching hole, how do bacillus thuringiensis toxins kill insect midgut cells? Cell Mol Life Sci. 2009;66(8):1337–49.
Gahan L, Gould F, Heckel DG. Identification of a gene associated with Bt resistance in Heliothis virescens. Science. 2001;293(5531):857–60. https://doi.org/10.1126/science.1060949.
Zhang H, Wu S, Yang Y, Tabashnik BE, Wu Y. Non-recessive Bt toxin resistance conferred by an intracellular cadherin mutation in field-selected populations of cotton bollworm. PLoS One. 2012;7(12):e53418. https://doi.org/10.1371/journal.pone.0053418.
Wang L, Ma, Y, Wan, P, Liu, K, Xiao, Y, Wang, J, Cong, S, Xu, D, Wu, K, Fabrick, J. A, Li, X. Resistance to Bacillus thuringiensis linked with a cadherin transmembrane mutation affecting cellular trafficking in pink bollworm from China. Insect Biochem Mol Biol 2018;94:28–35, DOI: https://doi.org/10.1016/j.ibmb.2018.01.004.
Wang L, Wang J, Ma Y, Wan P, Liu K, Cong S, et al. Transposon insertion causes cadherin mis-splicing and confers resistance to Bt cotton in pink bollworm from China. Sci Reports. 2019;9(1):1–10.
Jin L, Wang J, Guan F, Zhang J, Yu S, Liu S, et al. Dominant point mutation in a tetraspanin gene associated with field-evolved resistance of cotton bollworm to transgenic Bt cotton. Proc Natl Acad Sci U S A. 2018;115(46):11760–5. https://doi.org/10.1073/pnas.1812138115.
Herrero S, Gechev T, Bakker PL, Moar WJ, de Maagd RA. Bacillus thuringiensis Cry1Ca-resistant Spodoptera exigua lacks expression of one of four aminopeptidase N genes. BMC Genomics. 2005;6(1):96. https://doi.org/10.1186/1471-2164-6-96.
Tiewsiri K, Wang P. Differential alteration of two aminopeptidases N associated with resistance to bacillus thuringiensis toxin Cry1Ac in cabbage looper. Proc Natl Acad Sci. 2011;108(34):14037–42. https://doi.org/10.1073/pnas.1102555108.
Coates BS, Sumerford DV, Siegfried BD, Hellmich RL, Abel, C.A. Unlinked genetic loci control the reduced transcription of aminopeptidase N 1 and 3 in the European corn borer and determine tolerance to bacillus thuringiensis Cry1Ab toxin. Insect Biochem Mol Biol 2013;43(12):1152–1160.
Jurat-Fuentes JL, Karumbaiah L, Jakka SRK, Ning C, Liu C, Wu K, et al. Reduced levels of membrane-bound alkaline phosphatase are common to lepidopteran strains resistant to cry toxins from bacillus thuringiensis. PLoS One. 2011;6(3):e17606.
Baxter SW, Badenes-Pérez FR, Morrison A, Vogel H, Crickmore N, Kain W, et al. Parallel evolution of bacillus thuringiensis toxin resistance in Lepidoptera. Genetics. 2011;189(2):675–9. https://doi.org/10.1534/genetics.111.130971.
Atsumi S, Miyamoto K, Yamamoto K, Narukawa J, Kawai S, Sezutsu H, et al. Single amino acid mutation in an ATP-binding cassette transporter gene causes resistance to Bt toxin Cry1Ab in the silkworm, Bombyx mori. Proc Natl Acad Sci U S A. 2012;109(25):E1591–8. https://doi.org/10.1073/pnas.1120698109.
Yang X, Chen W, Song X, Ma X, Cotto-Rivera RO, Kain W, et al. Mutation of ABC transporter ABCA2 confers resistance to Bt toxin Cry2Ab in Trichoplusia ni. Insect Biochem Mol Biol. 2019;112:103209.
Mathew LG, Ponnuraj J, Mallappa B, Chowdary LR, Zhang J, Tay WT, et al. ABC transporter mis-splicing associated with resistance to Bt toxin Cry2Ab in laboratory-and field-selected pink bollworm. Sci Reports. 2018;8(1):1–15.
Coates BS, Siegfried BD. Linkage of an ABC transporter to a single QTL that controls Ostrinia nubilalis larval resistance to the bacillus thuringiensis Cry1Fa toxin. Insect Biochem Mol Biol. 2015;63:86–96.
Flagel LE, Lee YW. Wanjugi H, Swarup S, Brown a, Wang J, Kraft E, Greenplate J, Simmons J, Adams N., Wang Y. mutational disruption of the ABCC2 gene in fall armyworm, Spodoptera frugiperda, confers resistance to the Cry1Fa and Cry1A.105 insecticidal proteins. Sci Reports. 2018;8(1):1–11.
Pauchet Y, Bretschneider A, Augustin S, Heckel DG. A P-glycoprotein is linked to resistance to the bacillus thuringiensis Cry3Aa toxin in a leaf beetle. Toxins. 2016;8(12):362. https://doi.org/10.3390/toxins8120362.
Flagel LE, Swarup. S, Chen M, Bauer C, Wanjugi H, Carroll M, hill P, Tuscan M, Bansal R, Flannagan R, Clark TL. Genetic markers for western corn rootworm resistance to Bt toxin. G3: genes. Genomes, Genetics. 2015;5(3):399–405.
Forcada C, Alcacer E, Garcera MD, Tato A, Martinez R. Resistance to bacillus thuringiensis Cry1Ac toxin in three strains of Heliothis virescens, proteolytic and SEM study of the larval midgut. Arch Insect Biochem Physiol. 1999;42(1):51–63. https://doi.org/10.1002/(SICI)1520-6327(199909)42:1<51::AID-ARCH6>3.0.CO;2-6.
Martinez-Ramirez AC, Gould F, Ferre J. Histopathological effects and growth reduction in a susceptible and a resistant strain of Heliothis virescens (Lepidoptera, Noctuidae) caused by sublethal doses of pure Cry1A crystal proteins from bacillus thuringiensis. Biocontrol Sci Tech. 1999;9(2):239–46. https://doi.org/10.1080/09583159929811.
Castagnola A, Jurat-Fuentes JL. Intestinal regeneration as an insect resistance mechanism to entomopathogenic bacteria. Curr Opin Insect Sci. 2016;15:104–10. https://doi.org/10.1016/j.cois.2016.04.008.
Porta H, Muñoz-Minutti C, Soberón M, Bravo A. Induction of Manduca sexta larvae caspases expression in midgut cells by bacillus thuringiensis Cry1Ab toxin. Psyche. 2011;938249:1–7.
Hernández-Martínez P, Gomis-Cebolla J, Ferré J, Escriche B. Changes in gene expression and apoptotic response in Spodoptera exigua larvae exposed to sublethal concentrations of Vip3 insecticidal proteins. Sci Reports. 2017;7(1):1–12.
Portugal L, Muñóz-Garay C, de Castro DLM, Soberón M, Bravo A. Toxicity of Cry1A toxins from bacillus thuringiensis to CF1 cells does not involve activation of adenylate cyclase/PKA signaling pathway. Insect Biochem Mol Biol. 2017;80:21–31. https://doi.org/10.1016/j.ibmb.2016.11.004.
Qiu L, Fan J, Liu L, Zhang B, Wang X, Lei C, et al. Knockdown of the MAPK p38 pathway increases the susceptibility of Chilo suppressalis larvae to bacillus thuringiensis Cry1Ca toxin. Sci Reports. 2017;7:43964.
Park Y, Abdullah MA, Taylor MD, Rahman K, Adang MJ. Enhancement of bacillus thuringiensis Cry3Aa and Cry3Bb toxicities to coleopteran larvae by a toxin-binding fragment of an insect cadherin. Appl Environ Microbiol. 2009;75(10):3086–92. https://doi.org/10.1128/AEM.00268-09.
Park Y, Hua G, Ambati S, Taylor M, Adang MJ. Binding and synergizing motif within coleopteran cadherin enhances Cry3Bb toxicity on the Colorado potato beetle and the lesser mealworm. Toxins. 2019;11(7):386. https://doi.org/10.3390/toxins11070386.
Tan SY, Rangasamy M, Wang H, Vélez AM, Hasler J, McCaskill D, et al. RNAi induced knockdown of a cadherin-like protein (EF531715) does not affect toxicity of Cry34/35Ab1 or Cry3Aa to Diabrotica virgifera virgifera larvae (Coleoptera, Chrysomelidae). Insect Biochem Mol Biol. 2016;75:117–24. https://doi.org/10.1016/j.ibmb.2016.06.006.
Rault LC, Siegfried BD, Gassmann AJ, Wang H, Brewer GJ, Miller NJ. Investigation of Cry3Bb1 resistance and intoxication in western corn rootworm by RNA sequencing. J Appl Entomol. 2018;142(10):921–36. https://doi.org/10.1111/jen.12502.
Zhao Z, Meihls LN, Hibbard BE, Ji T, Elsik CG, Shelby KS. Differential gene expression in response to eCry3.1Ab ingestion in an unselected and eCry3.1Ab-selected western corn rootworm (Diabrotica virgifera virgifera LeConte) population. Sci Reports. 2019;9(1):4896.
Wang H, Eyun SI, Arora K, Tan SY, Gandra P, Moriyama E, et al. Patterns of gene expression in western corn rootworm (Diabrotica virgifera virgifera) neonates, challenged with Cry34Ab1, Cry35Ab1 and Cry34/35Ab1, based on next-generation sequencing. Toxins. 2017;9(4):124.
Le SQ, Gascuel O. An improved general amino acid replacement matrix. Mol Biol Evol. 2008;25(7):1307–20. https://doi.org/10.1093/molbev/msn067.
Adang MJ, Crickmore N, Jurat-Fuentes JL. Diversity of bacillus thuringiensis crystal toxins and mechanism of action. Adv Insect Physiol. 2014;47:39–87.
Oppert B, Dowd SE, Bouffard P, Li L, Conesa A, Lorenzen MD, et al. Transcriptome profiling of the intoxication response of Tenebrio molitor larvae to bacillus thuringiensis Cry3Aa protoxin. PLoS One. 2012;7(4):e34624.
Bel Y, Jakubowska AK, Costa J, Herrero S, Escriche B. Comprehensive analysis of gene expression profiles of the beet armyworm Spodoptera exigua larvae challenged with bacillus thuringiensis Vip3Aa toxin. PLoS One. 2013;8:e81927.
Sparks ME, Blackburn MB, Kuhar D, Gundersen-Rindal DE. Transcriptome of the Lymantria dispar (gypsy moth) larval midgut in response to infection by bacillus thuringiensis. PLoS One. 2013;8(5):e61190.
Vellichirammal NN, Wang H, Eyun S, Moriyama E, Coates BS, Miller NJ, Siegfried, B.D. Transcriptional analysis of susceptible and resistant European corn borer strains and their response to Cry1F protoxin. BMC Genomics 2015;16(1):558, DOI: https://doi.org/10.1186/s12864-015-1751-6.
Seong KM, Coates BS, Sun W, Clark JM, Pittendrigh BR. Changes in neuronal signaling and cell stress response pathways are associated with a multigenic response of Drosophila melanogaster to DDT selection. Genome Biol Evol. 2017;9(12):3356–72. https://doi.org/10.1093/gbe/evx252.
Zhang W, Chen J, Keyhani NO, Zhang Z, Li S, Xia Y. Comparative transcriptomic analysis of immune responses of the migratory locust, Locusta migratoria, to challenge by the fungal insect pathogen. Metarhizium acridum BMC Genomics. 2015;16(1):867. https://doi.org/10.1186/s12864-015-2089-9.
Gunning RV, Dang HT, Kemp FC, Nicholson IC, Moores GD. New resistance mechanism in Helicoverpa armigera threatens transgenic crops expressing bacillus thuringiensis Cry1Ac toxin. Appl Environ Microbiol. 2005;71(5):2558–63. https://doi.org/10.1128/AEM.71.5.2558-2563.2005.
Ma G, Rahman MM, Grant W, Schmidt O, Asgari S. Insect tolerance to the crystal toxins Cry1Ac and Cry2Ab is mediated by the binding of monomeric toxin to lipophorin glycolipids causing oligomerization and sequestration reactions. Dev Compar Immunol. 2012;37(1):184–92. https://doi.org/10.1016/j.dci.2011.08.017.
Tay WT, Mahon RJ, Heckel DG, Walsh TK, Downes S, James WJ, et al. Insect resistance to bacillus thuringiensis toxin Cry2Ab is conferred by mutations in an ABC transporter subfamily a protein. PLoS Genet. 2015;11(11):e1005534.
Guo Z, Kang S, Zhu X, Xia J, Wu Q, Wang S, et al. Down-regulation of a novel ABC transporter gene (Pxwhite) is associated with Cry1Ac resistance in the diamondback moth, Plutella xylostella (L.). Insect Biochem Mol Biol. 2015;59:30–40.
Zhang T, Coates BS, Wang Y, Wang Y, Wang Z, He K. Down-regulation of aminopeptidase N and ABC transporter subfamily G transcripts in Cry1Ab and Cry1Ac resistant Asian corn borer, Ostrinia furnacalis (Lepidoptera, Crambidae). Int J Biol Sci. 2017;13(7):835–51. https://doi.org/10.7150/ijbs.18868.
Lynch J, Fukuda Y, Krishnamurthy P, Du G, Schuetz JD. Cell survival under stress is enhanced by a mitochondrial ATP-binding cassette transporter that regulates hemoproteins. Cancer Res. 2009;69(13):5560–7. https://doi.org/10.1158/0008-5472.CAN-09-0078.
Krishnamurthy P, Schuetz JD. The ABC transporter Abcg2/Bcrp, role in hypoxia mediated survival. Biometals. 2005;18(4):349–58.
Leslie EM. Arseniceglutathione conjugate transport by the human multidrug resistance proteins (MRPs/ABCCs). J Inorg Biochem. 2012;108:141–9.
Liu W, Feng Q, Li Y, Ye L, Hu M, Liu Z. Coupling of UDP-glucuronosyltransferases and multidrug resistance-associated proteins is responsible for the intestinal disposition and poor bioavailability of emodin. Toxicol Appl Pharmacol. 2012;265(3):316e324.
Sau A, Pellizzari Tregno F, Valentino F, Federici G, Caccuri AM. Glutathione transferases and development of new principles to overcome drug resistance. Arch Biochem Biophys. 2010;500(2):116e122.
Dermauw W, Van Leeuwen T. The ABC gene family in arthropods, comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Biol. 2014;45:89–110. https://doi.org/10.1016/j.ibmb.2013.11.001.
Mason EF, Rathmell JC. Cell metabolism, an essential link between cell growth and apoptosis. Biochimica et Biophysica Acta (BBA)-Mol. Cell Res. 2011;1813(4):645–54.
Eming SA, Wynn TA, Martin P. Inflammation and metabolism in tissue repair and regeneation. Science. 2017;356(6342):1026–30. https://doi.org/10.1126/science.aam7928.
Scott JG, Wen Z. Cytochromes P450 of insects, the tip of the iceberg. Pest Manag Sci. 2001;57(10):958–67. https://doi.org/10.1002/ps.354.
Iga M, Kataoka H. Recent studies on insect hormone metabolic pathways mediated by cytochrome P450 enzymes. Biol Pharm Bull. 2012;35(6):838–43. https://doi.org/10.1248/bpb.35.838.
Sztal T, Chung H, Berger S, Currie P. D, Batterham P, Daborn PJ. A cytochrome P450 conserved in insects is involved in cuticle formation. PLoS One 2012;7(5):e36544.
Feyereisen R. Insect P450 enzymes. Annu Rev Entomol. 1999;44(1):507–33. https://doi.org/10.1146/annurev.ento.44.1.507.
Denisov IG, Makris TM, Sligar SG, Schlichting I. Structure and chemistry of cytochrome p450. Chem Rev. 2005;105(6):2253–77. https://doi.org/10.1021/cr0307143.
Redza-Dutordoir M, Averill-Bates DA. Activation of apoptosis signaling pathways by reactive oxygen species. Biochimica et Biophysica Acta (BBA)-Mol. Cell Res. 2016;1863(12):2977–92.
Khan AU, Wilson T. Reactive oxygen species as cellular messengers. Chem Biol. 1995;2(7):437–45.
Tian J, Hu J, Liu G, Yin H, Chen M, Miao P, et al. Altered gene expression of ABC transporters, nuclear receptors and oxidative stress signaling in zebrafish embryos exposed to CdTe quantum dots. Environ Pollution. 2019;244:588–99. https://doi.org/10.1016/j.envpol.2018.10.092.
Cooper DM, Granville DJ, Lowenberger C. The insect caspases. Apoptosis. 2009;14(3):247–56. https://doi.org/10.1007/s10495-009-0322-1.
Shi Y. Caspase activation, revisiting the induced proximity model. Cell. 2004;117(7):855–8. https://doi.org/10.1016/j.cell.2004.06.007.
Kornbluth S, White K. Apoptosis in Drosophila, neither fish nor fowl (nor man, nor worm). J Cell Sci. 2005;118(9):1779–87. https://doi.org/10.1242/jcs.02377.
Hay BA. Understanding IAP function and regulation, a view from Drosophila. Cell Death Differ 2000;7(11):1045–1056.
Harvey AJ, Soliman H, Kaiser WJ, Miller LK. Anti-and pro-apoptotic activities of baculovirus and Drosophila IAPs in an insect cell line. Cell Death Differ. 1997;4(8):733–44.
Dong,Z, Venkatachalam MA, Wang J, Patel Y, Saikumar P, Semenza GL, Force T, Nishiyama J. Up-regulation of apoptosis inhibitory protein IAP-2 by hypoxia HIF-1-independent mechanisms. J Biol Chem 2001;276(22):18702–18709.
Orme M, Meier P. Inhibitor of apoptosis proteins in Drosophila, gatekeepers of death. Apoptosis. 2009;14(8):950–60. https://doi.org/10.1007/s10495-009-0358-2.
Gubser C, Bergamaschi D, Hollinshead M, Lu X, van Kuppeveld FJM, Smith GL. A new inhibitor of apoptosis from vaccinia virus and eukaryotes. PLoS Pathol. 2007;3(2):17.
Moldoveanu T, Follis AV, Kriwacki RW, Green DR. Many players in BCL-2 family affairs. Trends Biochem Sci. 2014;39(3):101–11. https://doi.org/10.1016/j.tibs.2013.12.006.
Colussi PA, Quinn LM, Huang DC, Coombe M, Read SH, Richardson H, et al. Debcl, a proapoptotic Bcl-2 homologue, is a component of the Drosophila melanogaster cell death machinery. J Cell Biol. 2000;148(4):703–14. https://doi.org/10.1083/jcb.148.4.703.
Galindo KA, Lu WJ, Park JH, Abrams JM. The Bax/Bak ortholog in Drosophila, Debcl, exerts limited control over programmed cell death. Development. 2009;136(2):275–83. https://doi.org/10.1242/dev.019042.
Rojas-Rivera D, Hetz C. TMBIM protein family, ancestral regulators of cell death. Oncogene. 2015;34(3):269–80.
Xu Q, Reed JC. Bax inhibitor-1, a mammalian apoptosis suppressor identified by functional screening in yeast. Mol Cell. 1998;1(3):337–46.
Chae HJ, Kim,HR, Xu C, Bailly-Maitre B, Krajewska M, Krajewski S, Banares S, Cui J, Digicaylioglu M, Ke N, Kitada S. BI-1 regulates an apoptosis pathway linked to endoplasmic reticulum stress. Mol Cell 2004;15(3):355–366, DOI: https://doi.org/10.1016/j.molcel.2004.06.038.
Robinson KS, Clements A, Williams AC, Berger CN, Frankel G. Bax inhibitor 1 in apoptosis and disease. Oncogene. 2011;30(12):2391–400. https://doi.org/10.1038/onc.2010.636.
Walter P, Ron D. The unfolded protein response, from stress pathway to homeostatic regulation. Science. 2011;334(6059):1081–6. https://doi.org/10.1126/science.1209038.
Lee GH, Kim HK, Chae SW, Kim DS, Ha KC, Cuddy M, et al. Bax inhibitor-1 regulates endoplasmic reticulum stress-associated reactive oxygen species and heme oxygenase-1 expression. J Biol Chem. 2007;282(30):21618–28.
Hori O, Miyazaki M, Tamatani T, Ozawa K, Takano K, Okabe M, et al. Deletion of SERP1/RAMP4, a component of the endoplasmic reticulum (ER) translocation sites, leads to ER stress. Mol Cell Biol. 2006;26(11):4257–67. https://doi.org/10.1128/MCB.02055-05.
Li J, Yu Q, Zhang B, Xiao C, Ma T, Yi X, et al. Stress-associated endoplasmic reticulum protein 1 (SERP1) and Atg8 synergistically regulate unfolded protein response (UPR) that is independent on autophagy in Candida albicans. Int J Med Microbiol. 2018;308(3):378–86.
M’Angale PG, Staveley BE. Knockdown of the putative lifeguard homologue CG3814 in neurons of Drosophila melanogaster. Genet Mol Res GMR. 2016;15(4):1–11.
Cancino-Rodezno A, Porta H, Soberón M, Bravo A. Defense and death responses to pore forming toxins. Biotechnol Genet Engineering Rev. 2009;26(1):65–82. https://doi.org/10.5661/bger-26-65.
Branson TF. The selection of a non-diapause strain of Diabrotica virgifera (Coleoptera, Chrysomelidae). Entomol Exper et Appl. 1976;19(2):148–54.
Branson TF, Jackson JJ, Sutter GR. Improved method for rearing Diabrotica virgifera virgifera (Coleoptera, Chrysomelidae). J Econ Entomol. 1988;81:410–4.
Meihls LN, Higdon ML, Ellersieck MR, Hibbard BE. Selection for resistance to mCry3A-expressing transgenic corn in western corn rootworm. J Econ Entomol. 2011;104(3):1045–54. https://doi.org/10.1603/EC10320.
Scharf ME, Meinke LJ, Siegfried BD, Wright RJ, Chandler LD. Carbaryl susceptibility, diagnostic concentration determination, and synergism for US populations of western corn rootworm (Coleoptera, Chrysomelidae). J Econ Entomol. 1999;92(1):33–9.
Bolger AM, Lohse M, Usadel B. Trimmomatic, a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20. https://doi.org/10.1093/bioinformatics/btu170.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protocols. 2013;8(8):1494.
Zerbino DR, Birney E. Velvet, algorithms for de novo short read assembly using de Bruijn graphs. Genome Res. 2008;18(5):821–9. https://doi.org/10.1101/gr.074492.107.
Schulz MH, Zerbino DR, Vingron M, Birney E. Oases, robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28(8):1086–92.
Li W, Godzik A. Cd-hit, a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.
Huang X, Madan A. CAP3, a DNA sequence assembly program. Genome Res. 1999;9(9):868–77. https://doi.org/10.1101/gr.9.9.868.
Duhoux A, Carrère S, Gouzy J, Bonin L, Délye C. RNA-Seq analysis of rye-grass transcriptomic response to an herbicide inhibiting acetolactate-synthase identifies transcripts linked to non-target-site-based resistance. Plant Mol Biol. 2015;87(4–5):473–87. https://doi.org/10.1007/s11103-015-0292-3.
Parra G, Bradnam K, Korf I. CEGMA, a pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics. 2007;23(9):1061–7. https://doi.org/10.1093/bioinformatics/btm071.
Zdobnov EM, Tegenfeldt F, Kuznetsov D, Waterhouse RM, Simao FA, Ioannidis P, et al. OrthoDB v9. 1, cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs. Nucl Acids Res. 2016;45(D1):D744–9. https://doi.org/10.1093/nar/gkw1119.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO, assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2. https://doi.org/10.1093/bioinformatics/btv351.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DL. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10. https://doi.org/10.1016/S0022-2836(05)80360-2.
Gouzy J, Carrere S, Schiex T, Frame DP. Sensitive peptide detection on noisy matured sequences. Bioinformatics. 2009;25(5):670–1. https://doi.org/10.1093/bioinformatics/btp024.
Bateman A, Coin L, Durbin R, Finn RD, Hollich V, Griffiths-Jones S, et al. The Pfam protein families database. Nucl Acids Res. 2004;32(suppl_1):D138–41.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, et al. The Pfam protein family database, towards a more sustainable future. Nucl Acids Res (database issue). 2016;44(D1):D279–85. https://doi.org/10.1093/nar/gkv1344.
Finn RD, Clements J, Eddy SR. HMMER web server, interactive sequence similarity searching. Nucl Acids Res. 2011;39(suppl 2):W29–37. https://doi.org/10.1093/nar/gkr367.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9(4):357–9. https://doi.org/10.1038/nmeth.1923.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014;15(12):550.
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-Seq data. Genome Biol. 2010;11(3):R25. https://doi.org/10.1186/gb-2010-11-3-r25.
Benjamini Y, Hochberg Y. Controlling the false discovery rate, a practical and powerful approach to multiple testing. J R Stat Soc B (Methodological). 1995;57(1):289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.
Fang H. dcGOR, An R package for analysing ontologies and protein domain annotations. PLoS Comput. Biol. 2014;10(10):e1003929.
Adedipe F, Grubbs N, Coates BS, Wiegmann B, Lorenzen M. Structural and functional insights into the Diabrotica virgifera virgifera ATP-binding cassette (ABC) transporters gene family. BMC Genomics. 2019;20(1):899. https://doi.org/10.1186/s12864-019-6218-8.
Attrill H, Falls K, Goodman JL, Millburn GH, Antonazzo G, Rey AJ, et al. FlyBase consortium. FlyBase, establishing a Gene Group resource for Drosophila melanogaster. Nucl Acids Res. 2016;44(D1):D786–92.
Kumar S, Doumanis J. The fly caspases. Cell Death Differ. 2000;7(11):1039–44. https://doi.org/10.1038/sj.cdd.4400756.
Thompson JD, Gibson TJ, Higgins DG. Multiple sequence alignment using ClustalW and ClustalX. Curr Prot Bioinformatics. 2003;1:2–3.
Kumar S, Nei M, Dudley J, Tamura K. MEGA, a biologist-centric software for evolutionary analysis of DNA and protein sequences. Brief Bioinformatics. 2008;9(4):299–306. https://doi.org/10.1093/bib/bbn017.
Marchler-Bauer A, Bo Y, Han L, He J, Lanczycki CJ, Lu S, et al. CDD/SPARCLE: functional classification of proteins via subfamily domain architectures. Nucl Acids Res. 2017;45(D1):D200–3.
Krogh A, Larsson B, von Heijne G, Sonnhammer ELL. Predicting transmembrane protein topology with a hidden Markov model, application to complete genomes. J Mol Biol. 2001;305(3):567–80. https://doi.org/10.1006/jmbi.2000.4315.
The findings and conclusions in this publication are those of the authors and should not be construed to represent any official USDA or U.S. Government determination or policy. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply a recommendation or endorsement by USDA. USDA ARS is an equal opportunity provider and employer. We thank Chad Nielson, USDA ARS North Central Agricultural Research Laboratory, Brookings, SD for work to create the inbred strain, Ped12 used in this study. We also thank Corinne Da Silva, Genoscope - Centre National de Séquençage: Evry, Île-de-France, for making the SRA and TSA submissions.
Sequences were obtained from the French Genoscope “Large-scale DNA sequencing” call of 2010/2011. Part of the research output reported here was a contribution from the United States Department of Agriculture (USDA), Agricultural Research Service (ARS) (CRIS Project 5030–22000-019-00D, Ecologically-based Management of Arthropods in the Maize Agroecosystem). Portions of the analyses were supported by resources provided by the SCINet project of the USDA ARS project number 0500–00093–001-00-D.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. Trimmed reads obtained from Illumina sequencing libraries.
Putative orthology of assembled transcripts in the Diabrotica virgifera virgifera reference transcriptome. Results based on BLASTx hits to protein models from Drosophila melanogaster (Dm), Tribolium castaneum (Tc), Dendroctonus ponderosae (Dp), and Anaplophora glabripennis (Ag). Number of unique D. v. virgifera derived proteins with shared evidence across all comparisons (n = 12,474) and all coleopteran species (n = 3883) are highlighted. A total of 34,684 assembled D. v. virgifera transcripts had ≥ 1 hit across all sets of reference protein models.
Number of trimmed Illumina single-end (SE) and paired-end (PE) reads that aligned to the Diabrotica virgifera virgifera reference transcriptome across replicates of each RNA-seq library. 1X = reads that aligned uniquely; > 1X = reads that aligned to greater than once location.
Transcripts with significant differences in read counts (expression) between Diabrotica virgifera virgifera larvae feeding on transgenic maize expressing the insecticidal B. thuringiensis (Bt) Cry3Bb1 toxin (T8; Table 2) compared to control non-Bt maize (T7). For each transcript, indication of differential expression also being shared when larvae were exposed to Heterorhabditis bacteriophora (Hb) and Metarhizium anisopliae (Ma) is indicated with a “1”. Standard output are shown for DESeq2 (baseMean = the average of the normalized counts taken over all samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange; stat = Wald statistic; pvalue = Wald test P-value; BH_padj = Benjamini-Hochberg adjusted p-value) and EdgeR (logFC = log2 fold change between the groups; logCPM = the average log2-counts-per-million; PValue = the two-sided P-value; BF_FDR = Bonferroni adjusted P-value). Lack of significance for any given transcript in DESeq2 or EdgeR results shown as “-”. The 1055 transcripts differentially expressed in both DESeq2 and EdgeR estimates are shown above the horizonal line. Transcript annotations include presence (Y) or absence (N) of signal_peptide probability (signalp_prob) threshold of > 0.700 shown only for “complete” proteins (match length = 1.0 to Drosophila melanogaster (Dm) or Trobiolium castaneum (Tc) protein models). PFAMSCAN_PfamA search results (format transcript query start, transcript query end: frame, strand (+ or -), range of hit to PFAM domain/PFAM domain name/ E-value/percent identity). Predicted protein information given (transcript start: transcript end: frame: strand (+ or -): amino acid sequence). BLASTx query results to indicated databases shown in standard output format with “/” indicating no values for queries receiving of hits. The 9 transcripts in blue italicized text are putatively derived from Hb or Ma.
Transcripts with significant differences in read counts (expression) between Diabrotica virgifera virgifera larvae feeding on transgenic maize expressing the insecticidal B. thuringiensis (Bt) Gpp34/Tpp35Ab1 toxin (T5; Table 2) compared to control non-Bt maize (T7). For each transcript, indication of differential expression also being shared when larvae were exposed to Heterorhabditis bacteriophora (Hb) and Metarhizium anisopliae (MA) is indicated with a “1”. Standard output are shown for DESeq2 (baseMean = the average of the normalized counts taken over all samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange; stat = Wald statistic; pvalue = Wald test P-value; BH_padj = Benjamini-Hochberg adjusted P-value) and EdgeR (logFC = log2 fold change between the groups; logCPM = the average log2-counts-per-million; PValue = the two-sided P-value; BF_FDR = Bonferroni adjusted P-value). Lack of significance for any given transcript in DESeq2 or EdgeR results shown as “-”. The 1374 transcripts differentially expressed in both DESeq2 and EdgeR estimates are shown above the horizonal line. Transcript annotations include presence (Y) or absence (N) of signal_peptide probability (signalp_prob) threshold of > 0.700 shown only for “complete” proteins (match length = 1.0 to Drosophila melanogaster (Dm) or Trobiolium castaneum (Tc) protein models). PFAMSCAN_PfamA search results (format transcript query start, transcript query end: frame, strand (+ or -), range of hit to PFAM domain/PFAM domain name/ E-value/percent identity). Predicted protein information given (transcript start: transcript end: frame: strand (+ or -): amino acid sequence). BLASTx query results to indicated databases shown in standard output format with “/” indicating no values for queries receiving of hits. The 13 transcripts in blue italicized text are putatively derived from Hb or Ma.
Transcripts with significant differences in read counts (expression) between Diabrotica virgifera virgifera larvae infected with Heterorhabditis bacteriophora (Hb) (T6; Table 2) compared to control non-Bt maize (T7). For each transcript, standard output are shown for DESeq2 (baseMean = the average of the normalized counts taken over all samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange; stat = Wald statistic; pvalue = Wald test P-value; BH_padj = Benjamini-Hochberg adjusted P-value) and EdgeR (logFC = log2 fold change between the groups; logCPM = the average log2-counts-per-million; PValue = the two-sided P-value; BF_FDR = Bonferroni adjusted P-value). Lack of significance for any given transcript in DESeq2 or EdgeR results shown as “-”. The 1562 transcripts differentially expressed in both DESeq2 and EdgeR estimates are shown above the horizonal line. Transcript annotations include presence (Y) or absence (N) of signal_peptide probability (signalp_prob) threshold of > 0.700 shown only for “complete” proteins (match length = 1.0 to Drosophila melanogaster (Dm) or Trobiolium castaneum (Tc) protein models). PFAMSCAN_PfamA search results (format transcript query start, transcript query end: frame, strand (+ or -), range of hit to PFAM domain/PFAM domain name/ E-value/percent identity). Predicted protein information given (transcript start: transcript end: frame: strand (+ or -): amino acid sequence). BLASTx query results to indicated databases shown in standard output format with “/” indicating no values for queries receiving of hits.
Transcripts with significant differences in read counts (expression) between Diabrotica virgifera virgifera larvae infected with Metarhizium anisopliae (Ma) (T3; Table 2) compared to control non-Bt maize (T7). For each transcript, standard output are shown for DESeq2 (baseMean = the average of the normalized counts taken over all samples; log2FoldChange = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange; stat = Wald statistic; pvalue = Wald test P-value; BH_padj = Benjamini-Hochberg adjusted P-value) and EdgeR (logFC = log2 fold change between the groups; logCPM = the average log2-counts-per-million; PValue = the two-sided P-value; BF_FDR = Bonferroni adjusted P-value). Lack of significance for any given transcript in DESeq2 or EdgeR results shown as “-”. The 1199 transcripts differentially expressed in both DESeq2 and EdgeR estimates are shown above the horizonal line. Transcript annotations include presence (Y) or absence (N) of signal_peptide probability (signalp_prob) threshold of > 0.700 shown only for “complete” proteins (match length = 1.0 to Drosophila melanogaster (Dm) or Trobiolium castaneum (Tc) protein models). PFAMSCAN_PfamA search results (format transcript query start, transcript query end: frame, strand (+ or -), range of hit to PFAM domain/PFAM domain name/ E-value/percent identity). Predicted protein information given (transcript start: transcript end: frame: strand (+ or -): amino acid sequence). BLASTx query results to indicated databases shown in standard output format with “/” indicating no values for queries receiving of hits.
Dispersion of raw and DESeq2 adjusted read counts about empirical mean for comparisons of triplicate RNA-seq data within A) control maize (Cn; Treatment T7) and Cry3Bb1 maize (VT3; T8), and C) control maize (Cn; T7) and Gpp34/Tpp35Ab1 maize (Hx; T5). MA-plots of Log2 transformed fold-change estimates and normalized mean read count are shown for B) control maize (Cn; Treatment T7) and Cry3Bb1 maize (VT3; T8), and D) control maize (Cn; T7) and Gpp34/Tpp35Ab1 maize (Hx; T5), with datapoints (change in transcript read counts) surpassing a Benjamini and Hochberg (1995) adjusted false discovery rate (FDR) of ≤0.05 shown in red.
Predicted functional PFAM domains encoded by transcripts differentially expressed by Diabrioca virgifera virgifera exposed to Heterorhabditis bacteriophora compared to controls. Counts provided for instances when ≥3 transcripts received an annotation within at least one of the treatments. InterPro identification (InterPro_ID) are also given, with NA indicating corresponding InterPro_ID not available.
Predicted functional PFAM domains encoded by transcripts differentially expressed by Diabrioca virgifera virgifera exposed to Metarhizium anisopliae compared to controls. Counts provided for instances when ≥3 transcripts received an annotation within at least one of the treatments. InterPro identification (InterPro_ID) are also given, with NA indicating corresponding InterPro_ID not available.
Gene Ontology (GO) terms enriched among transcripts differentially expressed in Heterorhabditis bacteriophora (Hb) and Metarhizium anisopliae (Ma) treatments that were shared with treatments A) Cry3Bb1 and B) Gpp34/Tpp35Ab1. Significantly overrepresented GO terms are shown within categories biological process (BP) and cellular component (CC) (FDR ≤ 1.0E− 5) and molecular function (MF) at level 2 (FDR ≤ 1.0E− 7; grey bars). Categories listed by GO ID and GO term. Number of transcripts encoding each PFAM domain within each functional category are indicated (black bars).
Transcripts with significant differences in read counts (expression) between Diabrotica virgifera virgifera larvae feeding on transgenic maize expressing the insecticidal B. thuringiensis (Bt) Cry3Bb1 toxin (T8; Table 2) and Gpp34/Tpp35Ab1 (T5) compared to control non-Bt maize (T7). For each transcript, indication of differential expression also predicted when larvae were exposed to Heterorhabditis bacteriophora (Hb) and Metarhizium anisopliae (Ma) is indicated with a “1”. Standard output are shown for DeSeq2 (baseMean = the average of the normalized counts taken over all samples; log2FC = log2 fold change between the groups; lfcSE = standard error of the log2FoldChange; stat = Wald statistic; pvalue = Wald test P-value; BH_padj = Benjamini-Hochberg adjusted P-value) and EdgeR (logFC = log2 fold change between the groups; logCPM = the average log2-counts-per-million; PValue = the two-sided P-value; BF_FDR = Bonferroni adjusted P-value). Lack of significance for any given transcript in DESeq2 or EdgeR results shown as “-”. Mean log2FC among DESeq2 and EdgeR estimates are provided for each transcript. Transcript annotations include presence (Y) or absence (N) of signal_peptide probability (signalp_prob) threshold of > 0.700 shown only for “complete” proteins (match length = 1.0 to Drosophila melanogaster (Dm) or Trobiolium castaneum (Tc) protein models). PFAMSCAN_PfamA search results (format transcript query start, transcript query end: frame, strand (+ or -), range of hit to PFAM domain/PFAM domain name/ E-value/percent identity). Predicted protein information given (transcript start: transcript end: frame: strand (+ or -): amino acid sequence). BLASTx query results to indicated databases shown in standard output format with “/” indicating no values for queries receiving of hits.
Alignment of partial enzymatic domain sequences from Drosophila melanogaster caspases DRONC, DRED, DAMM, STRICA, DCP-1, and DRICE (accessions in footnotes) with putative Diabrotica virgifera virgifera orthologs from the reference transcriptome assembly (DIAVI02NNNN) in this study and RefSeq gene models (XP_0281NNNNN) from the D. v. virgifera genome assembly Dvir_v2.0 (GCA_003013835.2; GenBank accession PXJM00000000.2). Conserved residues are highlighted grey, with those involved in binding and catalysis in black and yellow, respectively.
Structural annotation of inhibitor of apoptosis proteins (IAPs) encoded by transcripts differentially expressed following Diabrotica virgifera virgifera exposure to Cry3Bb1 or Gpp34/Tpp35Ab1. Annotations based on structure defined by Hay et al. (2000).
Multiple protein sequence alignment for B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor (BI) proteins (BI) and Lifeguard 4-like (LFG4) orthologs.
Orthology of stress-induced endoplasmic reticulum protein 2 (SERP2) encoded by the differentially expressed transcript DIAVI057195 in Diabrotica virgifera virgifera larvae exposed to Cry3Bb1 and Gpp34/Tpp35Ab1.
About this article
Cite this article
Coates, B.S., Deleury, E., Gassmann, A.J. et al. Up-regulation of apoptotic- and cell survival-related gene pathways following exposures of western corn rootworm to B. thuringiensis crystalline pesticidal proteins in transgenic maize roots. BMC Genomics 22, 639 (2021). https://doi.org/10.1186/s12864-021-07932-4
- Cell fate
- Stress response
- Pesticidal proteins
- Gene expression