# 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

Brad S Coates, Emeline Deleury, Aaron J Gassmann, Bruce E Hibbard, Lance J Meinke, Nicholas J Miller, Jennifer Petzold-Maxwell, B. Wade Wade French, Thomas W Sappington, Blair D Siegfried, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Brad S Coates, Emeline Deleury, Aaron J Gassmann, Bruce E Hibbard, Lance J Meinke, 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, 2021, 22 (1), pp.639. 10.1186/s12864-021-07932-4. hal-03482846

HAL Id: hal-03482846

https://hal.inrae.fr/hal-03482846

Submitted on 16 Dec 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# 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

Brad S. Coates¹*, Emeline Deleury³, Aaron J. Gassmann⁴, Bruce E. Hibbard⁵, Lance J. Meinke⁶, Nicholas J. Miller⁷, Jennifer Petzold-Maxwell⁴, B. Wade French⁸, Thomas W. Sappington¹, Blair D. Siegfried² and Thomas Guillemaud³

## Abstract

**Background:** 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.

**Results:** 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.

<footnote>* Correspondence: brad.coates@usda.gov
¹USDA-ARS, Corn Insects & Crop Genetics Research Unit, 103 Genetics Laboratory, Iowa State University, Ames, IA 50011, USA
Full list of author information is available at the end of the article</footnote>
**Conclusions:** 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.

**Keywords:** Diabrotica, Apoptosis, Cell fate, Stress response, Pesticidal proteins, Gene expression

**Background**

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–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) [4]. Transgenic maize hybrids that express pesticidal proteins are widely used by growers in the United States and several other countries worldwide [5]. 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 diapause 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 [11]. In some regions, adults are sprayed aerially with insecticides to reduce egg-laying and hence larval populations the following year [12]. After its first detection in 2009 [13], a high proportion of *D. v. virgifera* larvae in field populations in the United States now show resistance to Cry3Bb1 [14] along with cross-resistance to the structurally similar mCry3A [15] 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)) [4] maize is also documented in *D. v. virgifera* field populations, but these phenotypes show no cross-resistance to Cry3 proteins [18–20]. Analogous lack of Gpp34/Tpp35Ab1 cross-resistance with mCry3Aa was shown in laboratory selected *D. v. virgifera* [21]. This resistance occurred to transgenic hybrids that express a “low-dose” of Bt proteins [22] 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 [25]. 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 [26]. 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–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 $Ca^{2+}$ [30], 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 [31] which inserts into the enterocyte membrane following interactions with cell surface aminopeptidase N [32] or alkaline phosphatase [33]. ATP binding cassette (ABC) transporters are also implicated as receptors for Cry proteins [34–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–42], including the glycosyl hydrolase, α-amylase, from *Tenebrio molitor* [43], but their roles in Bt intoxication remains unknown.

An alternative model for Cry protein mode of action proposes a $Mg^{2+}$-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 requirements for pore formation [46].

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–50]. An amino acid change in a transmembrane domain of tetraspanin 1 is associated with *Helicoverpa armigera* resistance to Cry1Ac [51]. Cry1A resistance is also associated with reduced expression of one or more aminopeptidase N paralog in *Spodoptera exigua* [52], *Trichoplusia ni* [53]
and *Ostrinia nubilalis* [54], and alkaline phosphatase in *Heliothis virescens* [55]. The ABC transporter subfamily C member, *abcc2*, is linked or associated with lepidopteran resistance to Cry1Ac in *H. virescens* [34], *Plutella xylostella* [56], *Bombyx mori* [57], and *T. ni* [58], Cry2Ab2 in *Pectinophora possypiella* [59], and Cry1Fa in *O. nubilalis* [60] and *Spodoptera frugiperda* [61]. An ABC subfamily B member 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 [62]. An *abcb* gene is also located in proximity to a single QTL for Cry3Bb1 resistance in *D. v. virgifera* [63].

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 [66]. Conversely, transcripts encoding proteins involved in apoptosis (programmed cell death) are significantly up-regulated in *Manduca sexta* cells by Cry1Ac [67], and in the midgut of *S. exigua* by exposure to the Bt vegetative insecticidal protein (Vip) 3A protein [68]. Apoptosis was induced in *Choristoneura fumiferana* cells by a mechanism involving the mitogen activated protein kinase (MAPK) protein p38 following Cry1Ac exposure [69], and RNAi-mediated knockdown of MAPK p38 in *Chilo suppressalis* led to increased larval susceptibility to Cry1Ca [70].

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 [21]. Kinetic data demonstrate that Cry3Bb binds strongly to specific domains of the cadherin protein and enhances toxicity in *D. v. virgifera* [71] and other Coleoptera [72]. 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 [73]. Estimates of differential gene expression shows no significant induction of cadherin in susceptible compared to resistant larvae fed Cry3Bb1 [74], eCry3.1Ab [75], or between susceptible larvae exposed and not exposed to Cry3Bb1 [74] or Gpp34/Tpp35Ab1 [76]. 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 [74] and eCry3.1Ab [75], 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 [74].

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.

## Results

### 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
<table><caption>**Table 1** Developmental stages and conditions (single replicates) in the normalized reference transcriptome from *Diabrotica virgifera* inbred line Ped12</caption><thead><tr><th>ID</th><th>Stage</th><th>Stage or condition</th><th>Description</th></tr></thead><tbody><tr><td>C01</td><td>Adult</td><td>Heat (30 min at 37 °C)</td><td>Heat shock</td></tr><tr><td>C02</td><td>Adult</td><td>Starved 72 h</td><td>Starvation</td></tr><tr><td>C03</td><td>Adult</td><td>Mated females</td><td>Post-mating</td></tr><tr><td>C04</td><td>Adult</td><td>Mated males</td><td>Post-mating</td></tr><tr><td>C05</td><td>Adult</td><td>Unmated females</td><td>Virgin</td></tr><tr><td>C06</td><td>Adult</td><td>Unmated males</td><td>Virgin</td></tr><tr><td>C07</td><td>1st instar</td><td>Larval development; early-stage</td><td>Larval instar</td></tr><tr><td>C08</td><td>2nd instar</td><td>Larval development; mid-stage</td><td>Larval instar</td></tr><tr><td>C09</td><td>3rd instar</td><td>Larval development; late-stage</td><td>Larval instar</td></tr><tr><td>C10</td><td>Adult</td><td>Unmated adults; mixed sex</td><td>Adult (early; virgin)</td></tr><tr><td>C11</td><td>Egg (&lt; 1 day)</td><td>Embryonic; early development</td><td>Embryonic (early)</td></tr><tr><td>C12</td><td>Egg (12–16 day)</td><td>Embryonic; late development</td><td>Embryonic (late)</td></tr><tr><td>C13</td><td>2nd and 3rd instar</td><td>Exposed to mCry3A maize Event MIR604 <sup>A</sup></td><td>Insecticidal protein toxin</td></tr><tr><td>C14</td><td>1st and 2nd instar</td><td>Exposed to Cry3Bb1 maize Event MON88017 <sup>B</sup></td><td>Insecticidal protein toxin</td></tr><tr><td>C15</td><td>1st and 2nd instar</td><td>Exposed to Gpp34/Tpp35Ab1 maize Event DAS-59122-7 <sup>C</sup></td><td>Insecticidal protein toxin</td></tr><tr><td>C16</td><td>1st and 2nd instar</td><td>Exposed to drought conditions</td><td>Weather stress</td></tr><tr><td>C17</td><td>2nd and 3rd instar</td><td>Exposed to <i>Heterorhabditis bacteriophora</i></td><td>Entomopathogen (nematode)</td></tr><tr><td>C18</td><td>2nd and 3rd instar</td><td>Exposed to <i>Metarhizium anisopliae</i></td><td>Entomopathogen (fungus)</td></tr><tr><td>C19</td><td>Adult</td><td>Exposed to Thiomethoxam insecticide</td><td>Chemical insecticide</td></tr><tr><td>C20</td><td>Adult</td><td>Exposed to Curcurbitacin feeding stimulant</td><td>Semiochemical attractant</td></tr></tbody></table>

A. Syngenta, Basel, Switzerland;

B. Bayer Crop Sciences/Monanto Company, St. Louis, MO, USA

C. Corteva/Pioneer, Indianapolis, IN, USA

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-

<table><caption>**Table 2** Treatments performed in triplicate on *Diabrotica virgifera virgifera* inbred line Ped12 prior to RNA-seq library construction and sequencing</caption><thead><tr><th>ID</th><th>Treatment</th><th>Treatment (exposure) description</th><th>Libraries (SRA Accessions)</th></tr></thead><tbody><tr><td>T1</td><td>DiapEgg</td><td>Eggs diapause conditions.</td><td>LibA, LibB, LibC (ERR2791371–27913713)</td></tr><tr><td>T2</td><td>1st-Instar</td><td>1st instars.</td><td>LibD, LibE, LibF (ERR2791374–2791376)</td></tr><tr><td>T3</td><td>Ma</td><td>2nd and 3rd instars exposed to <i>Metarhizium anisopliae</i> (Ma).</td><td>LibG, LibH, LibI (ERR2791378–2791380)</td></tr><tr><td>T4</td><td>Gut</td><td>3rd instars exposed to the non-Bt hybrid Pioneer 38B85; fore-, mid- and hindgut</td><td>LibJ, LibK, LibL (ERR2791381–2791383)</td></tr><tr><td>T5</td><td>Hx</td><td>3rd instars exposed to maize Event DAS-59122-7; Herculex® XTRA Mycogen hybrid 2 T789 (Gpp34/Tpp35Ab1) for 2 days.</td><td>LibM, LibN, LibO (ERR-2791384-2,791,386)</td></tr><tr><td>T6</td><td>Hb</td><td>2nd and 3rd instars infected with <i>Heterorhabditis bacteriophora</i> (Hb).</td><td>LibP, LibQ, LibR (ERR2791387–2791389)</td></tr><tr><td>T7</td><td>Cn</td><td>3rd instars exposed to non-Bt hybrid 38B85 (control).</td><td>LibS, LibT, LibU (ERR2791390–2791392)</td></tr><tr><td>T8</td><td>VT3</td><td>3rd instars exposed to Event MON88017; DeKalb VT TriplePro® hybrid (Cry3Bb1) for 2 days.</td><td>LibV, LibW, LibX (ERR2791393–2791395)</td></tr></tbody></table>

Library read data used for estimates of differential gene expression are highlighted gray. All experimental accessions are located under BioSample ERS2715388; BioProject PRJEB28633
<table><caption>**Table 3** Representation of the de novo assembled *Diabrotica virgifera virgifera* reference transcriptome</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td colspan="2"><strong>A) Assembly and predictions</strong></td></tr><tr><td>Number of trimmed reads</td><td>769,834,956</td></tr><tr><td>Paired end reads</td><td>598,089,648</td></tr><tr><td>Single end reads</td><td>171,745,308</td></tr><tr><td>Number of in silico normalized reads</td><td>82,501,989</td></tr><tr><td>Paired end reads</td><td>51,579,260</td></tr><tr><td>Single end reads</td><td>30,922,729</td></tr><tr><td>Total assembled transcripts (≥ 200 bp)</td><td>228,885</td></tr><tr><td>Non-redundant assembled transcripts*</td><td>116,070</td></tr><tr><td>Minimum transcript length (nt)</td><td>200</td></tr><tr><td>Maximum transcript length (nt)</td><td>33,331</td></tr><tr><td>Total assembly size (bases)</td><td>101,481,135</td></tr><tr><td>N50</td><td>1197</td></tr><tr><td>Mean transcript size (nt)</td><td>875</td></tr><tr><td>Transcripts with predicted peptide</td><td>56,656</td></tr><tr><td>Predicted peptides with PFAM annotation</td><td>23,728</td></tr><tr><td colspan="2"><strong>B) Comparative analyses</strong></td></tr><tr><td>Transcripts with a BLASTx hit**</td><td></td></tr><tr><td>SwissProt database (v2015-02-04)</td><td>20,961</td></tr><tr><td>Drosophila melanogaster peptides v5.46</td><td>17,322</td></tr><tr><td>Tribolium castaneum proteins v3.0</td><td>24,543</td></tr><tr><td>Dendroctonus ponderosae proteins</td><td>26,965</td></tr><tr><td>A. glabripennis proteins</td><td>25,340</td></tr><tr><td>Overall<sup>+</sup> covering ≥80% of subject length</td><td><strong>14,574</strong></td></tr><tr><td>Overall<sup>+</sup> putative full-length</td><td><strong>7568</strong></td></tr></tbody></table>

* after redundancy suppression;

** E-value threshold = 10⁻⁷ and HSP cut-off = 25 (AA);

+ independent and non-overlapping set

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 estimated 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 Log₂(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).
<figure><img src="image_4.png" /></figure>

There was a strong correlation between Log₂(fold-change) estimates between DESeq2 and EdgeR methods (R² = 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 Log₂(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 Log₂(fold-change) were highly correlated between EdgeR and DESeq2 (R² = 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. 1 b). 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
<figure><img src="image_5.png" /></figure>

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
<table><caption>**Table 4** Count of differentially expressed transcripts between treatment pairs (FDR ≤ 0.05 in DESeq2 and EdgeR)</caption><thead><tr><th rowspan="2"></th><th colspan="4">ID/Treatment (as in Table 2)</th></tr><tr><th>T7-T8</th><th>T7-T5</th><th>T7-T3</th><th>T7-T6</th></tr></thead><tbody><tr><td>Transcript partitions</td><td>Cn v. VT3</td><td>Cn s. Hx</td><td>Cn v. Ma</td><td>Cn v. Hb</td></tr><tr><td>Total differentially expressed</td><td>1064</td><td>1387</td><td>1199</td><td>1562</td></tr><tr><td>total up-regulated</td><td>950</td><td>1296</td><td>638</td><td>1269</td></tr><tr><td>total down-regulated</td><td>114</td><td>91</td><td>563</td><td>293</td></tr><tr><td>total with PFAM domains</td><td>773</td><td>1031</td><td>658</td><td>820</td></tr><tr><td>total with signal peptide</td><td>293</td><td>366</td><td>212</td><td>204</td></tr><tr><td>Endogenous* (total/Bt only)<sup>#</sup></td><td>1055/ 798</td><td>1374/1079</td><td>NA</td><td>NA</td></tr><tr><td>up-regulated (total/Bt only)<sup>#</sup></td><td>942/ 771</td><td>1291/1059</td><td>636/NA</td><td>1269/NA</td></tr><tr><td>down-regulated (total/Bt only)<sup>#</sup></td><td>113/ 27</td><td>83/ 20</td><td>563/NA</td><td>293/NA</td></tr></tbody></table>

(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 sub-families 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 lec tin 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 Log₂(fold-change) estimates for these transcripts were highly correlated between DESeq2 and EdgeR analyses (R² ≥ 0.8524; Fig. 4a). Following removal of 133 transcripts that were also differentially expressed in Hb

<footnote>* no significant BLASTx hits to Wolbachia, *Heterorhabditis bacteriophora* (Hb), or *Metarhizium anisopliae* (Ma) protein databases (E-values > 1.0⁻⁷)</footnote>

<footnote># Bt only values obtained by filtering differentially expressed transcripts from Cry3Bb1 or Gpp34/Tpp35Ab1 datasets which were also differentially expressed in Hb and Ma treatments</footnote>

<footnote>NA Not applicable</footnote>
<table><caption>**Table 5** Predicted PFAM domains encoded by transcripts differentially-expressed in susceptible *Diabrotica virgifera virgifera* Ped12</caption><thead><tr><th rowspan="2">PFAM domain name</th><th rowspan="2">PFAM domain description</th><th rowspan="2">PFAM_ID</th><th rowspan="2">InterPro_ID</th><th colspan="3">Comparisons</th></tr><tr><th>T7-T5</th><th>T7-T8</th><th>Both</th></tr></thead><tbody><tr><td>Sugar_tr</td><td>Sugar transporter</td><td>PF00083</td><td>IPR005828</td><td>Hx v Cn</td><td>26</td><td>17</td><td>14</td></tr><tr><td>p450</td><td>Cytochrome P450</td><td>PF00067</td><td>IPR001128</td><td></td><td>22</td><td>17</td><td>15</td></tr><tr><td>Lectin_C</td><td>Lectin C-type domain</td><td>PF00059</td><td>IPR001304</td><td></td><td>21</td><td>12</td><td>11</td></tr><tr><td>Inhibitor_l29</td><td>Cathepsin propeptide inhibitor domain</td><td>PF08246</td><td>IPR013201</td><td></td><td>19</td><td>22</td><td>14</td></tr><tr><td>Trypsin</td><td>Trypsin</td><td>PF00089</td><td>IPR001254</td><td></td><td>19</td><td>20</td><td>9</td></tr><tr><td>Peptidase_C1</td><td>Papain family cysteine protease</td><td>PF00112</td><td>IPR000668</td><td></td><td>18</td><td>21</td><td>13</td></tr><tr><td>CBM_14</td><td>Carbohydrate binding module</td><td>PF01607</td><td>IPR002557</td><td></td><td>17</td><td>15</td><td>12</td></tr><tr><td>LRR_8</td><td>Leucine-rich repeat</td><td>PF13855</td><td>IPR001611</td><td></td><td>15</td><td>11</td><td>10</td></tr><tr><td>COesterase</td><td>Carboxyesterase, type B</td><td>PF00135</td><td>IPR002018</td><td></td><td>14</td><td>18</td><td>10</td></tr><tr><td>adh_short</td><td>Short chain dehydrogenase</td><td>PF00106</td><td>IPR002347</td><td></td><td>13</td><td>5</td><td>3</td></tr><tr><td>EcKinase</td><td>Ecdysteroid kinase</td><td>PF02958</td><td>IPR004119</td><td></td><td>11</td><td>12</td><td>8</td></tr><tr><td>UDPGT</td><td>UDP glycosyltransferases</td><td>PF00201</td><td>IPR002213</td><td></td><td>11</td><td>11</td><td>8</td></tr><tr><td>Kunitz_BPTI</td><td>Kunitz/Bovine pancreatic</td><td>PF00014</td><td>IPR002223</td><td></td><td>10</td><td>3</td><td>3</td></tr><tr><td>MFS_1</td><td>Major facilitator superfamily</td><td>PF07690</td><td>IPR011701</td><td></td><td>10</td><td>10</td><td>7</td></tr><tr><td>DUF1397</td><td>Unknown function</td><td>PF11901</td><td>IPR024518</td><td></td><td>9</td><td>3</td><td>3</td></tr><tr><td>Serpin</td><td>Serpin serine protease inhibitor</td><td>PF00079</td><td>IPR026796</td><td></td><td>9</td><td>3</td><td>3</td></tr><tr><td>Abhydrolase_1</td><td>Alpha/beta hydrolase fold</td><td>PF00561</td><td>IPR000073</td><td></td><td>8</td><td>11</td><td>6</td></tr><tr><td>Glyco_hydro_18</td><td>Glycosyl hydrolase family 18</td><td>PF00704</td><td>IPR001223</td><td></td><td>7</td><td>5</td><td>4</td></tr><tr><td>LRR_5</td><td>BspA type Leucine rich repeat region</td><td>PF13306</td><td>IPR026906</td><td></td><td>7</td><td>4</td><td>4</td></tr><tr><td>EF-hand_7</td><td>EF-hand domain pair</td><td>PF13499</td><td>IPR002048</td><td></td><td>6</td><td>3</td><td>3</td></tr><tr><td>Fibrinogen_C</td><td>Fibronectin</td><td>PF00147</td><td>IPR002181</td><td></td><td>6</td><td>7</td><td>4</td></tr><tr><td>GST_C</td><td>Glutathione S-transferase, C-terminal</td><td>PF00043</td><td>IPR004046</td><td></td><td>6</td><td>7</td><td>5</td></tr><tr><td>PBP_GOBP</td><td>PBP/GOBP family</td><td>PF01395</td><td>IPR006170</td><td></td><td>6</td><td>2</td><td>2</td></tr><tr><td>ADH_zinc_N</td><td>Zinc-binding dehydrogenase</td><td>PF00107</td><td>IPR013149</td><td></td><td>5</td><td>3</td><td>2</td></tr><tr><td>Aldo_ket_red</td><td>Aldo/keto reductase</td><td>PF00248</td><td>IPR023210</td><td></td><td>5</td><td>6</td><td>3</td></tr><tr><td>Ank_2</td><td>Ankyrin repeats</td><td>PF12796</td><td>IPR020683</td><td></td><td>5</td><td>4</td><td>3</td></tr><tr><td>Far-17a_AIG1</td><td>Far-17a_AIG1-like protein</td><td>PF04750</td><td>IPR006838</td><td></td><td>5</td><td>4</td><td>4</td></tr><tr><td>GILT</td><td>γ-interferon-inducible lysosomal thiol</td><td>PF03227</td><td>IPR004911</td><td></td><td>5</td><td>5</td><td>4</td></tr><tr><td>Glyco_hydro_1</td><td>Glycosyl hydrolase family 1</td><td>PF00232</td><td>IPR001360</td><td></td><td>5</td><td>14</td><td>5</td></tr><tr><td>Ig_2</td><td>Immunoglobulin domain</td><td>PF13895</td><td>IPR007110</td><td></td><td>5</td><td>0</td><td>0</td></tr><tr><td>JHBP</td><td>Juvenile hormone-binding protein</td><td>PF06585</td><td>IPR010562</td><td></td><td>5</td><td>10</td><td>4</td></tr><tr><td>Lipase</td><td>Lipase</td><td>PF00151</td><td>IPR013818</td><td></td><td>5</td><td>4</td><td>3</td></tr><tr><td>Ras</td><td>Ras family</td><td>PF00071</td><td>IPR001806</td><td></td><td>5</td><td>3</td><td>3</td></tr><tr><td>7tm_2</td><td>7 transmembrane receptor</td><td>PF00002</td><td>IPR000832</td><td></td><td>4</td><td>4</td><td>4</td></tr><tr><td>ABC2_membrane</td><td>ABC-2 type transporter</td><td>PF01061</td><td>IPR013525</td><td></td><td>4</td><td>2</td><td>1</td></tr><tr><td>ABC_tran</td><td>ABC transporter</td><td>PF00005</td><td>IPR003439</td><td></td><td>4</td><td>2</td><td>1</td></tr><tr><td>Abhydro_lipase</td><td>Alpha/beta hydrolase lipase region</td><td>PF04083</td><td>IPR006693</td><td></td><td>4</td><td>6</td><td>2</td></tr><tr><td>Asp</td><td>Eukaryotic aspartyl protease</td><td>PF00026</td><td>IPR033121</td><td></td><td>4</td><td>3</td><td>2</td></tr><tr><td>DUF3421</td><td>Protein of unknown function</td><td>PF11901</td><td>IPR024518</td><td></td><td>4</td><td>4</td><td>3</td></tr><tr><td>GST_N</td><td>Glutathione S-transferase, N-terminal</td><td>PF02798</td><td>IPR004045</td><td></td><td>4</td><td>5</td><td>3</td></tr><tr><td>I-set</td><td>Immunoglobulin I-set domain</td><td>PF07679</td><td>IPR013098</td><td></td><td>4</td><td>4</td><td>3</td></tr><tr><td>Lys</td><td>C-type lysosome/alpha-lactalbumin</td><td>PF00062</td><td>IPR001916</td><td></td><td>4</td><td>0</td><td>0</td></tr></tbody></table>
<table><caption>**Table 5** Predicted PFAM domains encoded by transcripts differentially-expressed in susceptible *Diabrotica virgifera virgifera* Ped12 (Continued)</caption><thead><tr><th></th><th></th><th></th><th colspan="3">Comparisons</th></tr></thead><tbody><tr><td>Sina</td><td>Seven in abstensia protein family</td><td>PF03145</td><td>IPR018121</td><td>4</td><td>4</td><td>3</td></tr><tr><td>SVWC</td><td>von Willebrand factor type C</td><td>PF15430</td><td>IPR029277</td><td>4</td><td>4</td><td>4</td></tr><tr><td>Tetraspannin</td><td>Tetraspannin family</td><td>PF00335</td><td>IPR018499</td><td>4</td><td>2</td><td>2</td></tr></tbody></table>

Counts predicted within translated products for differentially-expressed transcripts for Bt susceptible *D. v. virgifera* Ped12 larval exposures to transgenic maize roots that express Gpp34/Tpp35Ab1 (T5; Herculex® (Hx) hybrid) and Cry3Bb1 (T8; VT TriplePro® (VT3) hybrid) compared to control (T7). The intersection of differentially-expressed transcripts in both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments also indicated. All counts representative of transcripts after in silico subtraction of those also differentially-expressed in exposures to the entomopathogens *Heterorhabditis bacteriophora* (Hb) and *Metarhizium anisopliae* (Ma). Values provided for instances when ≥4 transcripts received an annotation within at least one of the treatments. InterPro identification (InterPro_ID) are also given

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

<table><caption>**Table 6** Candidate *Bacillus thuringiensis* (Bt) insecticidal protein receptors differentially in susceptible *Diabrotica virgifera virgifera* Ped12 larvae</caption><thead><tr><th rowspan="2">Transcript</th><th rowspan="2">Dvir_v2.0 Protein model(s)#</th><th rowspan="2">Putative annotation</th><th colspan="4">Log<sub>2</sub>(fold-change)</th><th rowspan="2">Tribolium castaneum ortholog</th></tr><tr><th>T7-T8</th><th>T7-T5</th><th>T7-T6</th><th>T7-T3</th></tr></thead><tbody><tr><td>DIAVI009217</td><td>XP_028139024.1 &amp; 25.1</td><td>ABC subfamily G</td><td>+ 2.42</td><td>NS</td><td>NS</td><td>NS</td><td>XP_973526.1; 07047; TcABCG-4H; 53%/ 646</td></tr><tr><td>DIAVI012990</td><td>XP_028134948.1 to 52.1</td><td>ABC subfamily C</td><td>+ 2.41</td><td>+ 2.28</td><td>NS</td><td>NS</td><td>XP_969781.1; 14,775; TcABCC-ST; 57%/ 439</td></tr><tr><td>DIAVI001032</td><td>XP_028145606.1</td><td>ABC subfamily C</td><td>+ 2.07</td><td>NS</td><td>NS</td><td>NS</td><td>XP_8193834.1; 14,403; TcABCC-SR; 44%/1374</td></tr><tr><td>DIAVI006431</td><td>XP_028139807.1</td><td>ABC subfamily B</td><td>NS</td><td>+ 1.98</td><td>NS</td><td>NS</td><td>XP_974441.2; 15,192; TcABCB-6A; 77%/ 612</td></tr><tr><td>DIAVI025188</td><td>XP_028135097.1</td><td>ABC subfamily G</td><td>NS</td><td>+ 1.96</td><td>NS</td><td>NS</td><td>XP_008198312.1; ND; ND; 36%/ 156</td></tr><tr><td>DIAVI008609</td><td>XP_028129515.1 &amp; 16.1</td><td>ABC subfamily E</td><td>NS</td><td>+ 1.58</td><td>NS</td><td>NS</td><td>XP_968009.1; 10,519; TcABCE-3A; 91%/ 557</td></tr><tr><td>DIAVI001096</td><td>XP_028128982.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>+ 0.81</td><td>NS</td><td>XP_972534.1; 09892; TcABCC-7B; 63%/1280</td></tr><tr><td>DIAVI001196</td><td>XP_028139572.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>+ 2.39</td><td>+ 2.13</td><td>XP_971802.2; 14,383; TcABCC-5 N; 44%/1260</td></tr><tr><td>DIAVI001308</td><td>XP_028155748.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>+ 1.57</td><td>NS</td><td>XP_971908.1; 14,385; TcABCC-5P; 46%/1266</td></tr><tr><td>DIAVI001589</td><td>XP_028142713.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>+ 1.16</td><td>+ 1.03</td><td>XP_971802.2; 14,383; TcABCC-5 N; 45%/1280</td></tr><tr><td>DIAVI001763</td><td>XP_028131534.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>+ 1.11</td><td>+ 1.22</td><td>XP_971802.2; 14,383; TcABCC-5 N; 48%/1258</td></tr><tr><td>DIAVI001938</td><td>XP_028155747.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>-1.38</td><td>-2.05</td><td>XP_971908.1; 14,385; TcABCC-5P; 47%/1240</td></tr><tr><td>DIAVI012237</td><td>XP_028137191.1 to 93.1</td><td>ABC subfamily C</td><td>NS</td><td>NS</td><td>NS</td><td>-2.17</td><td>XP_971802.2; 14,383; TcABCC-5 N; 47%/ 587</td></tr><tr><td>DIAVI021979</td><td>XP_028140422.1</td><td>Tetraspanin</td><td>+ 2.48</td><td>+ 2.34</td><td>NS</td><td>NS</td><td>XP_969428.111703; CD63 antigen;53%/ 271</td></tr><tr><td>DIAVI004770</td><td>XP_028137187.1 &amp; 88.1</td><td>Tetraspanin</td><td>NS</td><td>+ 1.84</td><td>NS</td><td>NS</td><td>XP_966752.212071; Tetraspanin33; 91%; 291</td></tr></tbody></table>

Fold-change and direction (+ or -) of differential expression among transcripts is indicated within exposures to roots of maize hybrids that express Bt pesticidal proteins Cry3Bb1 (T8; VT TriplePro® (VT3) hybrid) and Gpp34/Tpp35Ab1 (T5; Herculex® (Hx) hybrid), or entomopathogens *Heterorhabditis bacteriophora* (T6) and *Metarhizium anisopliae* (T3), compared to the non-Bt maize control treatment (T7)

# Protein models in RefSeq accession GCF_003013835.1 annotated from assembly accession GCA_003013835.2 under BioProject PRJNA432972; NS: not significant with FDR > 0.05; ND: not determined; Accn: GenBank accession; GLEAN: *Tribolium castaneum* GLEAN protein nomenclature; Ortholog: *T. castaneum*. ABC transporter classification according to Broehan et al. (2013) and putative annotation as defined by Adepipe et al. (2019); ID: identity from BLASTp alignments; Aligned L: length of BLASTp alignment
<table><caption>**Table 7** Differential expressed apoptotic and autophagy-related protein-encoding transcripts</caption><thead><tr><th rowspan="2">Transcript</th><th rowspan="2">Putative annotation</th><th colspan="4">Log<sub>2</sub>(fold-change)</th></tr><tr><th>T7-T8</th><th>T7-T5</th><th>T7-T6</th><th>T7-T3</th></tr></thead><tbody><tr><td>DIAVI022989</td><td>Initiator caspase-like (orthologs: STRICA or DAMM)</td><td>+ 3.21</td><td>+ 3.05</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI027204</td><td>Effector caspase-like (ortholog: DECAY)</td><td>+ 3.00</td><td>+ 2.47</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI057195</td><td>Stress-associated endoplasmic reticulum protein 2</td><td>+ 2.26</td><td>+ 2.13</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI029891</td><td>Lifeguard 4-like (BAX inhibitor domain protein)</td><td>+ 2.02</td><td>+ 1.76</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI013501</td><td>Autophagy-related gene 13 (atg13)-encoded protein</td><td>+ 2.20</td><td>+ 1.71</td><td>-1.42</td><td>-1.52</td></tr><tr><td>DIAVI011561</td><td>Programmed cell death protein 4</td><td>ns</td><td>+ 1.67</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI026079</td><td>BAX inhibitor 1 (BI-1)</td><td>ns</td><td>+ 1.52</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI011972</td><td>Inhibitor of apoptosis 1 (IAP1) isoform X1</td><td>ns</td><td>+ 2.22</td><td>ns</td><td>ns</td></tr><tr><td>DIAVI007715</td><td>Inhibitor of apoptosis 1 (IAP1) isoform X5</td><td>+ 4.36</td><td>ns</td><td>ns</td><td>ns</td></tr></tbody></table>

Response among susceptible *Diabrotica virgifera virgifera* Ped12 larvae exposed to *Bacillus thuringiensis* (Bt) maize roots and entomopathogens. Fold-change and direction (+ or -) of differential expression for transcripts is indicated for treatments with exposures to maize roots that express Bt pesticidal proteins Cry3Bb1 (T8; VT TriplePro® (VT3) hybrid) and Gpp34/Tpp35Ab1 (T5; Herculex® (Hx) hybrid), and to the entomopathogens *Heterorhabditis bacteriophora* (T6) and *Metarhizium anisopliae* (T3), compared to non-Bt control (T7). Transcripts differentially expressed in the intersection of both Bt maize treatments are highlighted in gray. # Putative homologs from protein models in RefSeq accession GCF_003013835.1 annotated from genome assembly accession GCA_003013835.2 under the BioProject PRJNA432972: DIAVI022989 = XP_028132548.1; DIAVI027204 = XP_028144295.1; DIAVI057195 = XP_028135864.1; DIAVI029891 = XP_028155171.1; DIAVI013501 = XP_028130453.1; DIAVI011561 = XP_028128855.1 and XP_028128856.1; DIAVI026079 = XP_028141533.1; DIAVI011972 and DIAVI007715 = XP_028141650.1 to 58.1

<figure><img src="image_6.png" /></figure>

**Fig. 3** Enriched Gene Ontology (GO) terms among transcripts differentially expressed in susceptible *Diabrotica virgifera virgifera* larvae exposed to maize roots expressing *Bacillus thuringiensis* pesticidal proteins. For independent **A) Cry3Bb1**\* and **B) Gpp34/Tpp35Ab1** treatments\*, significantly over-represented GO terms are shown within the categories of biological process (BP) and cellular component (CC) (FDR ≤ 1.0E⁻⁵), and molecular function (MF) at level 2 (FDR ≤ 1.0E⁻⁶; gray bars). Categories listed by GO ID and GO term in order of highest significance (bottom) to lowest significance (top). Number of transcripts encoding each PFAM domain within each functional category are indicated (black bars). \* Filtered differentially expressed transcripts following in silico subtraction that removed those shared with transcripts also differentially expressed within entomopathogen treatments (Hb and Ma)
<figure><img src="image_7.png" /></figure>

putatively most enriched for those localized in the extra-
cellular space, secretory vesicles, and microbodies, and
having coenzyme binding function and are involved in
small molecule catabolism (Fig. 3). These same categor-
ies are also most significantly enriched among differen-
tially 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 upregu-
lated in both Cry3Bb1 and Gpp34/Tpp35Ab1 treatments
(Table 7).
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Enriched Gene Ontology (GO) terms among transcripts within the intersection of differentially expressed in susceptible *Diatrotica virgifera* larvae independently exposed to maize roots expressing different *Bacillus thuringiensis* (Bt) pesticidal proteins*. Significantly over-represented GO terms within the categories of biological process (BP), cellular component (CC) (FDR ≤ 1.0E⁻⁵), and molecular function (MF) at level 2 (FDR ≤ 1.0E⁻⁶; gray bars). Categories listed by GO ID and GO term in order of highest significance (bottom) to lowest significance (top). Number of transcripts encoding each PFAM domain within each functional category are indicated (black bars). * Set of 563 transcripts with estimated significant levels of differential expression between *D. v. virgifera* larvae exposed to hybrid maize roots that express Bt Cry3BB1 (T8; Hx) and Gpp34/Tpp35Ab1 (T5; Hx) compared to non-Bt control maize roots (T7; Cn) within independent treatments (intersection). These shared transcripts were further filtered by in silico subtraction to remove differentially expressed transcripts that were also differentially expressed within entomopathogen treatments (*Heterorhabditis bacteriophora* (Hb) and *Metarhizium anisopliae* (Ma))</figcaption></figure>

## 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
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Phylogenetic relationship and orthology of putative *Diabrotica virgifera virgifera* caspases. Maximum-likelihood (ML) analysis of a partial enzymatic domain alignment between putative *D. v. virgifera* caspases with seven known caspases from *Drosophila melanogaster* (Supplementary Fig. 5). The *D. v. virgifera* transcripts (DIAVI0NNNN) and protein models (XP_0281NNNN.1) from RefSeq accession GCF_003013835.1 for the draft genome assembly accession GCF_003013835.1 Dvir_v2 are integrated for homolog estimation. Caspase-encoding transcripts DIAVI027204 and DIAVI022989 up-regulated in both Cry3Bb1 and Cry34/35Ab1 exposed larvae (Table 7) are enclosed in boxes. Clades corresponding to effector and initiator classes are indicated by brackets. The ML of the constructed tree using the LG model of protein sequence evolution (Le and Gascuel 2008) [77] with a gamma shape parameter (G) = 2.7553 and 6.77% of sites defined as evolutionarily invariable () was minimized at a log likelihood of -2224.28, and resulted in a tree with a total branch length of 10.157. Node support obtained using 1000 bootstrap pseudo-replications of the aligned protein sequence data</figcaption></figure>

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.

# Discussion

## 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 [27]. 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–76, 79–82]. This phenomenon is also observed as part of insect responses to chemical insecticide exposures [83] and pathogen infections [84]. Such outcomes not only impose challenges for interpreting results, but also in formulating hypotheses to guide future research.

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 [51]. 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 [51], 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* [87], 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 [75]. 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 [76], or where transcripts were not detectable in gut tissues [74]. 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–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 [95], 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* [75], 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 [74], 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 [74] 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 [98], including regulation of insect ecdysone and juvenile hormone pathways [99], cuticle formation [100], and xenobiotics detoxification [101]. Uncoupling of P450 oxygenation reactions results in production of reactive oxygen species (ROS), hydrogen peroxide and superoxide [102]. Cellular homeostasis can become disrupted during times of cell stress when excess ROS is produced due to high P450 activity. ROS can trigger apoptosis [103], or act as second messengers that modulate other cellular processes [104]. Stress responses triggered by ROS during high metabolic states are intimately tied to increased ABC transporter activities [105], 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 [106]. 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 [67] and in midgut tissues of *S. exigua* following a sub-lethal exposure to the Bt Vip3A protein [68].

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 [107]. 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 [108]. 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 [110], suggesting their functional conservation. IAPs are up-regulated in response to cell stress conditions [111], 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 [112], 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) [113], 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 [114]. 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 [115], but evidence suggests it has a limited role in triggering apoptosis [116]. TMBIM proteins modulate stress through several intracellular mechanisms [117]. 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 [120]. 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 [121]. Therefore, BI-1 promotes cell survival pathways by suppressing factors that would otherwise promote apoptotic signaling [122]. SERP also functions within the UPR by enhancing protein stability [123, 124]. Somewhat analogously, LFG4 remediates stress in the Golgi and ER [117], 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* [125]. BI-1 and LFG4 structure and function are remarkably conserved, where ectopic expression of viral orthologs can rescue knockdown phenotypes in mammalian cells [113], 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 in-
dividual contributions and impacts on cellular and
organismal outcomes.

Conclusions

A greater understanding of cellular responses among
susceptible insects to pesticidal protein exposure can in-
form future research into mechanisms that lead to re-
sistance. To date mutations in cell surface receptors
have mainly been implicated in pesticidal protein resist-
ance among insects, although some evidence suggests al-
teration 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 pat-
terns 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 de-
termine cell fate. Because we used whole larvae, we can-
not predict any possible tissue- or cell type-specific
responses. In general, exposure level has a role in cell re-
sponse to pore forming proteins suggesting a “high-
dose” leads to death by oncosis (cell swelling and bleb-
bing due to osmotic imbalance) as opposed to a “low-
dose” that tends to trigger apoptosis [126]. Therefore,
the responses of D. v. virgifera to the “low-dose” pesti-
cidal protein exposures characteristic of current trans-
genic 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.

Methods

Samples, treatments, and collections

A Bt susceptible non-dia pausing strain of D. v. virgifera
[127] maintained at the United States Department of
Agriculture, Agricultural Research Service, North Cen-
tral 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 F₁ to F₅; (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 main-
tained thereafter as a closed colony of approximately
1000 individuals for 8 generations prior to use in this
study.

Ped12 individuals were sampled at different develop-
mental 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 [128] at
USDA-NCARL. For conditions C13 to C16 (Table 1)
maize seedling mats were germinated from seeds of hy-
brids expressing modified mCry3A (Event MIR604; Syn-
genta, 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. virgi-
fera 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 de-
scribed 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
[129]. 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) re-
ceived no water (Table 1).

Additionally, Ped12 larvae were separately exposed to
the entomopathogenic nematode *Heterorhabditis bacter-
iophora* (Hb) strain BU (Becker Underwood, Ames, IA,
USA; C17 in Table 1), and the entomopathogenic fun-
gus, *Metarhizium anisopliae* (Ma) strain F52 (provided
by Stefan Jaronski, USDA-ARS; C18 in Table 1). In brief,
300 µl of a 510 Hb ml⁻¹ suspension was aliquoted onto
filter paper in 10 cm Petri plates, and six 2nd and 3rd in-
stars were added and exposed for 48 h. Analogously,
300 µl of 1.5 × 10⁷ Ma conidia ml⁻¹ 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 expos-
ure assay [130]. For this, 20 ml vials were coated with
sublethal levels of the neonicotinoid insecticide, thia-
methoxam (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 [128] 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 [131]. 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. 1 a). 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 [132] with parameter *--maxIPCov* 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 [133] and OASES 0.2.06 [134]. 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 [135] 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 [136] and NRCL tools with a progressive reduction in the stringency parameters [137]. 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) [138]. Transcripts were mapped to 248 CEGs, and the set of 1066 universal single copy orthologs from Arthropoda obtained from OrthoDB v 9 [139] as identified using BUSCO v 3 [140] (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 [141] was used for querying, and the results filtered for those with an *E*-value ≤ 10<sup>-7</sup> 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 [142]. Putative protein family domains were assigned by searches of the PFAM A database v. 27.0 [143, 144] using the program HMMSEARCH [145] 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. 1 b). 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 [146] 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 [147] and EdgeR v3.8.6 [148]. 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 [149]. 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/Tpp35Ab1 exposure

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 [150]. Subsequent enrichment analyses by dcGOR 1.0.6 at GO level 2 applied significance thresholds of *E*-values ≤ 1.0<sup>-5</sup> for CC and MF terms, and ≤ 1.0<sup>-6</sup> 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 [151]. 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 [152]. This was repeated for *H. armigera* TSPAN1.

## Phylogenetics and structural annotations

Protein sequences of *Drosophila melanogaster* caspases were downloaded from FlyBase [152] by gene name search: DRONC, DRED, DAMM, STRICA, DECAY, DCP-1, and DRICE, [153]. 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 [141] (*E*-value cutoffs ≤ 10<sup>-20</sup>). 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 [154] within the MEGA8.0 alignment utility [155] (default parameters). The LG + G + I model of protein sequence evolution [77] 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 [155].

Sequences from up-regulated transcripts encoding putative inhibitor of apoptosis proteins (IAPs; DIAV I011972 and DIAVI007715), B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor (BI) proteins BI-1 (DIAVI026079) and Lifeguard 4-like (LFG4; DIAV I029891), 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 [140] (*E*-value cutoffs ≤ 10<sup>-60</sup>). 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) [156] 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 [109]. 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 [157].

## Abbreviations

Bt: *B. thuringiensis*; BI-1: B-cell-lymphoma protein 2 (Bcl-2)-associated X (BAX) inhibitor-1; Cry: crystalline protein; Gpp: AeGerolysin-like pesticidal protein; LFG4: Lifeguard 4; SERP2: stress-associated ER protein 2; Tpp: toxin-10-like pesticidal protein

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12864-021-07932-4.

**Additional file 1: Supplementary Table S1.** Trimmed reads obtained from Illumina sequencing libraries.

**Additional file 2: Supplementary Fig. S1.** 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.

**Additional file 3: Supplementary Table S2.** 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.

**Additional file 4: Supplementary Table S3.** 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; ifcSE = standard error of the log2Fold-Change; 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 horizontal line. Transcript annotations include presence (Y) or absence (N) of signal肽 probability (signalp Prob) threshold of > 0.700 shown only for "complete" proteins (match length = 1.0 to *Drosophila melanogaster* (Dm) or *Tribolium 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.

**Additional file 5: Supplementary Table S4.** 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; log2Fold-Change = log2 fold change between the groups; IfcSE = 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 horizontal 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 *Troiolium 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.

**Additional file 6: Supplementary Table S5.** 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; IfcSE = standard error of the log2Fold-Change; 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 horizontal 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 *Troiolium 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.

**Additional file 7: Supplementary Table S6.** 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; IfcSE = standard error of the log2Fold-Change; 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 horizontal 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 *Troiolium 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.

(signalp PROB) threshold of > 0.700 shown only for "complete" proteins (match length = 1.0 to *Drosophila melanogaster* (Dm) or *Troiolium 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.

**Additional file 8: Supplementary Fig. S2.** 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.

**Additional file 9: Supplementary Table S7.** Predicted functional PFAM domains encoded by transcripts differentially expressed by *Diabrotica 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.

**Additional file 10: Supplementary Table S8.** Predicted functional PFAM domains encoded by transcripts differentially expressed by *Diabrotica 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.

**Additional file 11: Supplementary Fig. S3.** 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).

**Additional file 12: Supplementary Table S9.** 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; IfcSE = 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 *Troiolium 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.
databases shown in standard output format with "/" indicating no values for queries receiving of hits.

**Additional file 13: Supplementary Fig. S4.** 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_0281NNNN) 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.

**Additional file 14: Supplementary Fig. S5.** 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).

**Additional file 15: Supplementary Fig. S6.** 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.

**Additional file 16: Supplementary Fig. S7.** 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.

## Acknowledgements

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.

## Authors' contributions

BSC conceived and performed parts of the study, analyzed data, interpreted results and wrote the manuscript. ED assembled, processed RNA-seq data, generated annotations, and wrote part of the manuscript. BWF, BEH, AJG, LJM, NJM, JPM provided materials and conceived and/or performed experiments. TWS, BDS, and TG conceived parts of the study, advised, and/or provided funding. All co-authors read and approved the content of this research.

## Funding

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 SCI-Net project of the USDA ARS project number 0500-00093-001-00-D.

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

## Declarations

### Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that they have no competing interests.

## Author details

¹USDA-ARS, Corn Insects & Crop Genetics Research Unit, 103 Genetics Laboratory, Iowa State University, Ames, IA 50011, USA. ²University of Florida, Gainesville, FL, USA. ³INRAE, CNRS, Université Côte d'Azur, ISA, France.

⁴Department of Entomology, Iowa State University, Ames, IA, USA.

⁵USDA-ARS, Plant Genetics Research Unit, Columbia, MO, USA. ⁶Department of Entomology, University of Nebraska, Lincoln, NE, USA. ⁷Illinois Institute of Technology, Chicago, IL, USA. ⁸USDA-ARS, North Central Agricultural Research Laboratory, Brookings, SD, USA.

Received: 7 April 2021 Accepted: 4 August 2021

Published online: 04 September 2021

## References

1. Roush RT, McKenzie JA. Ecological genetics of insecticide and acaricide resistance. Annu Rev Entomol. 1987;32(1):361-80. https://doi.org/10.1146/a.nnurev.32.010187.002045.

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

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

4. Crickmore N, Berry C, Panneerselvam S, Mishra R, Connor TR, Bonning BC. A structure-based nomenclature for *bacillus thuringiensis* and other bacteri-derived pesticidal proteins. J Invertebr Pathol. 2021. https://doi.org/10.1016/j.jip.2020.107438.

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

6. Peralta C, Palma L. Is the insect world overcoming the efficacy of *bacillus thuringiensis* Toxins. 2017;1:39.

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

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

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

10. Urias-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.

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

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

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

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

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

16. 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.
17. 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.

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

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

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

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

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

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

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

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

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

27. Pigott CR, Ellar DJ. Role of receptors in *bacillus thuringiensis* crystal toxin activity. Microbiol Mol Biol Rev. 2007;71(2):255–81.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

45. Zhang X, Candas M, Griko NB, Taussig R, Bulla LA. A mechanism of cell death involving an adenyllyl 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.

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

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

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

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

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

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

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

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

54. 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.
55. 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.

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

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

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

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

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

61. 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 Cry1A105 insecticidal proteins. Sci Reports. 2018;8(1):1–11.

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

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

64. 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-A RCH6>3.0.CO;2-6.

65. 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/0958315992811.

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

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

68. Hernández-Martínez P, Gomis-Cebolla J, Ferré J, Escrige 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.

69. Portugal L, Muñoz-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.

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

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

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

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

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

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

76. Wang H, Eyun SL, 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.

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

78. Adang MJ, Crickmore N, Jurat-Fuentes JL. Diversity of *bacillus thuringiensis* crystal toxins and mechanism of action. Adv Insect Physiol. 2014;47:39–87.

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

80. Bel Y, Jakubowska AK, Costa J, Herrero S, Escrige 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.

81. Sparks ME, Blackburn MB, Kuhar D, Gundersen-Rindal DE. Transcriptome of the *Lymannia dispar* (gypsy moth) larval midgut in response to infection by *bacillus thuringiensis*. PLoS One. 2013;8(5):e61190.

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

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

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

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

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

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

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

89. 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.1 8868.

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

91. Krishnamurthy P, Schuetz JD. The ABC transporter Abcg2/Bcrp, role in hypoxia mediated survival. Biometals. 2005;18(4):349–58.

92. Leslie EM. Arseniceglutathione conjugate transport by the human multidrug resistance proteins (MRPs/ABCCs). J Inorg Biochem. 2012;108:141–9.

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

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

95. Dermauw W, Van Leeuwen T. The ABC gene family in arthropods, comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Cell. 2014;45:89–110. https://doi.org/10.1016/j.ibmb.2013.11.001.

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

97. Eming SA, Wynn TA, Martin P. Inflammation and metabolism in tissue repair and regenstation. Science. 2017;356(6342):1026–30. https://doi.org/10.1126/science.aam7928.

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

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

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

101. Feyereisen R. Insect P450 enzymes. Annu Rev Entomol. 1999;44(1):507–33. https://doi.org/10.1146/annurev.ento.44.1.507.

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

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

104. Khan AU, Wilson T. Reactive oxygen species as cellular messengers. Chem Biol. 1995;2(7):437–45.

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

106. Cooper DM, Granville DJ, Lowenberger C. The insect caspases. Apoptosis. 2009;14(3):247–56. https://doi.org/10.1007/s10495-009-0322-1.

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

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

109. Hay BA. Understanding IAP function and regulation, a view from *Drosophila*. Cell Death Differ 2000;7(11):1045–1056.

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

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

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

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

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

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

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

117. Rojas-Rivera D, Hetz C. TMBIM protein family, ancestral regulators of cell death. Oncogene. 2015;34(3):269–80.

118. Xu Q, Reed JC. Bax inhibitor-1, a mammalian apoptosis suppressor identified by functional screening in yeast. Mol Cell. 1998;1(3):337–46.

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

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

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

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

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

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

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

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

127. Branson TF. The selection of a non-diapause strain of *Diabrotica virgifera* (Coleoptera, Chrysomelidae). Entomol Exper et Appl. 1976;19(2):148–54.

128. Branson TF, Jackson JJ, Sutter GR. Improved method for rearing *Diabrotica virgifera virgifera* (Coleoptera, Chrysomelidae). J Econ Entomol. 1988;81:410–4.

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

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

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

132. Haas BJ, Papanicolou 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.

133. Zerbino DR, Birney E. Velvet, algorithms for de novo short read assembly using de Brujin graphs. Genome Res. 2008;18(5):821–9. https://doi.org/10.1110/gr.074492.107.

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

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

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

137. Duhoux A, Carrère S, Gouzy J, Bonin L, Délye C. RNA-Seq analysis of ryegrass 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.

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

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

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

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

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

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

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

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

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

147. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014;15(12):S50.

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

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

150. Fang H. dcGOR, An R package for analysing ontologies and protein domain annotations. PLoS Comput. Biol. 2014;10(10):e1003929.

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

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

153. Kumar S, Doumanis J. The fly caspases. Cell Death Differ. 2000;7(11):1039-44. https://doi.org/10.1038/sj.cdd.4400756.

154. Thompson JD, Gibson TJ, Higgins DG. Multiple sequence alignment using ClustalW and ClustalX. Curr Prot Bioinformatics. 2003;1:2-3.

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

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

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

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_10.png" /></figure>