Next Article in Journal
Phospholipase D from Loxosceles laeta Spider Venom Induces IL-6, IL-8, CXCL1/GRO-α, and CCL2/MCP-1 Production in Human Skin Fibroblasts and Stimulates Monocytes Migration
Next Article in Special Issue
Baseline Susceptibility of Field Populations of Helicoverpa armigera to Bacillus thuringiensis Vip3Aa Toxin and Lack of Cross-Resistance between Vip3Aa and Cry Toxins
Previous Article in Journal / Special Issue
Use of Carabids for the Post-Market Environmental Monitoring of Genetically Modified Crops
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Patterns of Gene Expression in Western Corn Rootworm (Diabrotica virgifera virgifera) Neonates, Challenged with Cry34Ab1, Cry35Ab1 and Cry34/35Ab1, Based on Next-Generation Sequencing

1
Department of Agronomy and Horticulture, University of Nebraska-Lincoln, Lincoln, NE 68583-0915, USA
2
Center for Biotechnology, School of Biological Sciences, UNL, Lincoln, NE 68583, USA
3
Dow AgroSciences, Indianapolis, IN 46268, USA
4
Monsanto, St. Louis, MO 63167, USA
5
Entomology and Nematology Department, University of Florida, Gainesville, FL 32611-0620, USA
*
Author to whom correspondence should be addressed.
Toxins 2017, 9(4), 124; https://doi.org/10.3390/toxins9040124
Submission received: 26 February 2017 / Revised: 22 March 2017 / Accepted: 27 March 2017 / Published: 30 March 2017
(This article belongs to the Special Issue The Insecticidal Bacterial Toxins in Modern Agriculture)

Abstract

:
With Next Generation Sequencing technologies, high-throughput RNA sequencing (RNAseq) was conducted to examine gene expression in neonates of Diabrotica virgifera virgifera (LeConte) (Western Corn Rootworm, WCR) challenged with individual proteins of the binary Bacillus thuringiensis insecticidal proteins, Cry34Ab1 and Cry35Ab1, and the combination of Cry34/Cry35Ab1, which together are active against rootworm larvae. Integrated results of three different statistical comparisons identified 114 and 1300 differentially expressed transcripts (DETs) in the Cry34Ab1 and Cry34/35Ab1 treatment, respectively, as compared to the control. No DETs were identified in the Cry35Ab1 treatment. Putative Bt binding receptors previously identified in other insect species were not identified in DETs in this study. The majority of DETs (75% with Cry34Ab1 and 68.3% with Cry34/35Ab1 treatments) had no significant hits in the NCBI nr database. In addition, 92 DETs were shared between Cry34Ab1 and Cry34/35Ab1 treatments. Further analysis revealed that the most abundant DETs in both Cry34Ab1 and Cry34/35Ab1 treatments were associated with binding and catalytic activity. Results from this study confirmed the nature of these binary toxins against WCR larvae and provide a fundamental profile of expression pattern of genes in response to challenge of the Cry34/35Ab1 toxin, which may provide insight into potential resistance mechanisms.

1. Introduction

The western corn rootworm (WCR), Diabrotica virgifera virgifera LeConte, is an important pest of field corn, Zea mays L. [1,2,3] both in terms of crop losses and costs associated with management practices. Managing corn rootworm populations to minimize risk of economic loss is extremely difficult, due in part to their unique capacity to evolve resistance to a variety of management practices including chemical insecticides [4,5,6], cultural control practices such as crop rotation [7,8], and transgenic corn hybrids expressing the Cry3Bb1 toxin [9,10]. The binary toxin, Cry34/Cry35Ab1 represents two proteins that are co-expressed in transgenic corn to control WCR and are commercialized both as single event hybrids and more recently as pyramided events with either Cry3Bb1 or mCry3A. However, with resistance to Cry3 toxins documented among WCR field populations [10,11], even the pyramided events may rely exclusively on the binary Cry34/35Ab1 toxin.
It has been suggested that Cry34Ab1 and Cry35Ab1 have specific binding sites on the brush border membrane of the rootworm midgut, and Cry34Ab1 enhances the specific binding of Cry35Ab1 [12]. Cry34Ab1 has limited toxicity by itself while Cry35Ab1 alone has no toxicity, and the binary toxin is necessary to achieve mortality of rootworm larvae in diet bioassays and to achieve root protection with transgenic maize plants [13]. The specific binding receptor(s) for these two Bt toxins have yet to be identified.
Next generation sequencing provides a simple and comprehensive approach to measure changes in expression of genes in response to environmental stressors such as insecticides [14] and Bt proteins [15,16,17,18]. The goal of this study was to identify genes responsive in western corn rootworms challenged with its individual components as well as the binary toxin, Cry34/35Ab1 which has been commercialized in transgenic corn since 2006 and which still performs effectively in the field although incomplete resistance to Cry34/Cry35Ab1 was recently reported [19]. The objective of the present research is to develop an overview of WCR genes responsive to Bt Cry34/Cry35Ab1 intoxication.

2. Results

2.1. Next Generation Sequencing

The sequencing run using the Illumina HiSeq2000 yielded a total of ~1345 million paired-end raw reads. The number of reads generated for each treatment is presented in Table S1. After removal of low quality reads (Q < 20), 287 million to 383 million (~97%) high-quality reads remained (Table S2) and were used for further analysis. All raw read data are available at the NCBI Sequence Read Archive (SRA) under Accession SRP037561 [20].

2.2. Mapping and Differential Expression Gene Analysis

Using the Bowtie program, nearly 70% of filtered reads for all four experimental conditions were aligned to the WCR reference transcriptome. The output (read counts) from alignment were applied directly in DESeq, and MA-plots were generated to depict a general view of the distribution of the differentially expressed transcripts (DETs) (p < 0.05) for all three treatments as compared to the control treatment (Figure 1). In general, transcripts with altered expression varied considerably with the source of Cry toxin challenge (Figure 1b,c). Among the three treatments, the Cry34/35Ab1 treatment produced almost 10,000 sequences defined as DETs (Figure 1c). In contrast, less than 1000 DETs were detected in the Cry34Ab1 treatment (Figure 1b) and none were detected in the treatment with Cry35 alone (Figure 1a).
The initial read count datasets were further processed by filtering transcripts with low read counts (cpm > 1) with edgeR, 26,218; 29,109 and 29,520 transcripts were remained in the Cry35Ab1, Cry34Ab1 and Cry34/35Ab1 compared to control comparisons, respectively, and further used as input in differential expression analysis.
As illustrated in Table 1, in total, 116, 132 and 135 DETs were identified with DESeq, edgeR and limma, respectively, in the Cry34Ab1treatment as compared to the control.
Among them, 114 DETs were commonly identified by these three methods. In the combined treatment with Cry34/35Ab1, 2215, 1673 and 2336 genes were classified as DETs by DESeq, edgeR and limma, respectively, and 1300 DETs were commonly identified, which is at least 10-fold greater than the number of DETs caused by exposure to Cry34Ab1 alone.
These 1300 DETs were used for subsequent analyses. As previously described, no DETs were detected by any of the three methods in the treatment with Cry35Ab1 alone.
Among the DETs identified in Cry34Ab1 (114) and Cry34/35 Ab1 (1300) treatment (Table 1), 92 DETs (31 up-regulated and 61 down-regulated) were found to be in common (Table 2). The remainder (22 and 1108 DETs) were assigned only to Cry34Ab1 or Cry34/35Ab1 treatment, respectively. Among the DETs common to both treatments, the average fold change was generally higher in the Cry34/35Ab1 treatment (Table 2).
For the unique DETs in each treatment, the average fold was similar (2.36 for Cry34Ab1 and 2.69 for Cry34/35Ab1 up-regulated genes and 3.16 and 3.3 for down-regulated genes). However, the range of fold change of DETs was greater in the Cry34/35Ab1 treatment (~6.46-fold) as compared to that in Cry34Ab treatment (3.57-fold) in down-regulated category.

Annotation of DETs

Nearly 70% of DETs in Cry34Ab1 treated WCR (114) and in Cry34/35Ab1 treated WCR (1300) treatment had no significant hits in the NCBI non-redundant (nr) database (Figure 2a,b). Among the DETs with hits, 35 DETs in the Cry34Ab1 treatments and 385 DETs in Cry34/35Ab1 treatments were well annotated with most hits to other Coleopterans including Tribolium castaneum and Dendroctonus ponderosae.
Gene ontology classification of the DETs is provided in Figure 3 (Cry34Ab1 compared to control) and Figure 4 (Cry34/35Ab1 compared to control).
For the Cry34Ab1 treatment, the largest number of DETs was assigned to molecular function (17 transcripts) and cellular component (17 transcripts). Of those assigned to molecular function, catalytic activity (11 transcripts, 64.7%) and binding (6 transcripts, 35.3%) accounted for the largest number of genes assigned while all 7 transcripts represented in biological process were associated with metabolic process. In contrast, the largest group of DETs for the Cry34/35Ab1 treatment was associated with biological process (1007 transcripts), in which metabolic process (217 transcripts, 21.5%), cellular process (175 transcripts, 17.4%), and single-organism process (159 transcripts, 15.8%) accounted for the largest categories. The second largest group involved molecular function (531 transcripts), in which binding (196 transcripts, 10.7%) and catalytic activity (292 transcripts, 16%) were most abundant.
In addition to the differences in number of transcripts that responded to Cry34/35Ab1 as compared to Cry34Ab1 alone, there were 18 more additional functional activities associated with these differences in the Cry34/35Ab1 treatment (Table S3). Activities included antioxidant activity, molecular transducer activity, enzyme regulatory activity, receptor activity, transporter activity and membrane-enclosed lumen.
Transcripts coding for putative Bt toxin receptors identified from other insect species, such as cadherin, aminopeptidase N and ATP-binding cassette transporter (ABC) and metalloprotease were not detected with either Cry34Ab1 or Cry34/35Ab exposure. However, two different alkaline phosphatases (Dv_137932_c0_seq1 and Dv_149197_c0_seq1), which have been associated with Bt toxin binding in Lepidoptera, were identified as differentially expressed. The Dv_ 137932_c0_seq1 was down-regulated 7.46- and 2-fold in both Cry34/35 and Cry34Ab1 treatment, respectively. The Dv_149197_c0_seq1was up-regulated 1.65- and 1.2-fold in Cry34Ab1 and Cry34/35Ab1 treatment, respectively, based on DESeq results only.

2.3. GO-Term Enrichment and Pathway Analysis

As illustrated in Table S4, five significantly overrepresented GO terms were associated with DETs in Cry34Ab1 challenge. Four of these were related to molecular function, in which two GO terms were up-regulated and correlated with zinc ion binding (GO:0008270) and transition metal ion binding (GO:0046914). The other two were down-regulated and were related to hydrolase activity (GO:0004553 and GO:0016798). The remaining GO term (GO:0005975) was down-regulated and associated with carbohydrate metabolism.
In the Cry34/35Ab1 treatment, a total of 168 GO terms were significantly enriched (Tables S5 and S6). Among them, 152 GO terms (Table S5) were up-regulated and 16 (Table S6) were down-regulated. Of the up-regulated GO terms, 35 (23%) were associated with molecular function and mostly related to binding, such as ATP binding (GO:0005542), cation binding (GO:0043169), anion binding (GO:0043168), ion binding (GO:0043167) and carbohydrate derivative binding (GO:0097367). The remaining 117 GO terms were identified as biological processes and the GO terms associated with regulation accounted for the largest group (23 GO terms, 19.6%), including regulation of signaling (GO:0023051), regulation of Ras protein signal transduction (GO:0046578), regulation of hydrolase activity (GO:0051336) and regulation of lipid catabolic processes. Interestingly, only two GO terms were related to cellular process (GO:0009987) and cellular metabolic process (GO:0044237) (Table S6) were found to be under-represented as down-regulated transcripts in the Cry34/35Ab1 treatment. Moreover, no cellular component-related GO term was found to be enriched in either the Cry34Ab1 or Cry34/35Ab1 treatment.
Due to the limited annotation of the reference transcriptome, KEGG analysis was also conducted to identify pathways associated with DETs to help identify higher-level functions. For the Cry34Ab1 treatment, 7 out of 114 (6%) DETs were assigned in two pathways: (1) amino sugar and nucleotide sugar metabolism and (2) starch and sucrose metabolism. In contrast, almost 323 out of 1300 (25%) DETs associated with the Cry34/35Ab1 treatment were associated with 42 different pathways (Table S7). Among them, the majority (10 in up-regulated DETs and 22 in down-regulated DETs) of the identified were related to “metabolic pathway”. The top two pathways with most designated DETs were pyrimidine metabolism with 7 DETs and glycan degradation with 6 DETs. In addition, 6 pathways were with 5 DEGs, 9 pathways with 4 DETs, 4 pathways with 3 DETs, 9 pathways with 2 DETs assigned and the rest 28 pathways with 1 DET assigned only. In addition, among these 42 pathways identified, two pathways related to detoxification, drug metabolism-cytochrome P450 and glutathione metabolism assigned with a transcript (Dv_138610_c0_seq1) were included with a 1.61-fold change in the down-regulated category in Cry34/34Ab1 treatment.

2.4. Validation with RT-qPCR

For all primers used in the validation experiment for four genes, a primer efficiency value between 92.1% and 104.3% at R2 (correlation coefficient) > 0.99 was obtained (Table S8).
As shown in Table 3, the qPCR results indicated that the expression of GH45 and ALP was repressed, whereas the transcripts corresponding to the GSC and PAT were enriched in Cry34/35Ab-treated neonates. The gene expression based qPCR was in agreement with digital results from RNAseq analysis for all four genes tested.

3. Discussion

Differences in expression after exposure to individual toxins and the combination of the two toxins support the binary nature of the Cry34/35Ab1 toxin [21,22] as co-expression of both components in transgenic corn hybrids is necessary for control of corn rootworms. It has been recently shown that the individual toxins exhibit different binding characteristics in the midgut of rootworm and that Cry34Ab1 serves to enhance the Cry35Ab1 specific binding. Cry35Ab1 alone exhibits very low binding capacity in the absence of Cry34Ab1 [12]. Susceptible neonates challenged by exposure to individual components of the binary toxin and with an effective level of the binary toxin resulted in changes in expression that are consistent with their respective toxicities [22,23]. The combined results indicate that 1208 unique DETs were altered in their expression when challenged with the Cry34/35Ab1 combination, which was ~54 times greater than that observed in treatment with Cry34Ab1 alone (22 unique DETs). Cry35Ab1, which is non-toxic to rootworm larvae, did not cause differences in gene expression. These results are consistent with the known toxicity pattern for the two individual toxins and their combination where Cry35Ab1 is non-toxic alone, Cry34Ab1 reduces growth with much less larval mortality and Cry34/35Ab1 causes significant growth inhibition and larval mortality. The combination of the two toxins is necessary to achieve high mortality and their combined expression in transgenic maize is critical to protect roots from damage.
GO term enrichment analysis has been previously employed to identify potential pathways that respond to environmental stressors such as Bt toxins [16,18,24] and insecticides [14,25] and to classify the functions of the predicted proteins in a number of different insects [17,26]. Like other stressors, the Bt toxins induced changes in WCR gene expression [17,27,28,29]. However, the number of genes altered and GO terms enriched vary greatly with the source of challenge and target organism [16,17,24,28,30], exposure time [18], and even within different populations of the same species. For example, in two Cry1Ac resistant Plutella xylostella populations originating from different collections, the number of enriched GO terms associated with each population was markedly different [17]. In the current study, some of the genes responsive to Cry34/35Ab1 exposure, such as those involving zinc ion binding (GO:0008270), lipase activity (GO:0016298), catalytic activity (GO:0050790), cell communication (GO:0010646) and Ras protein signal transduction (GO:0046578), were also reported in other Bt challenge studies [16,17,18,24], suggesting expression of these genes may be common to Bt toxin exposure.
A number of putative Bt protein receptors reported among other coleopterans include cadherin-like proteins [31], ADAM metalloprotease [32] and β-glucosidase [33], and in Lepidoptera include alkaline phosphatase, cadherin, aminopeptidase N, and ABC transporters [34]. In the current study, cadherin and amino peptidase N were not detected as DETs. However, one ADAM metalloprotease (Dv_149203_c0_seq1) and β-glucosidase precursors (Dv_139888_c0_seq1) were down-regulated in both Cry34Ab1 and Cry34/35Ab1 treatment, indicating that expression was repressed.
Although not reported from other coleopterans, alkaline phosphatase has been documented as Bt receptor and associated with Bt resistance in a number of lepidopterans [35,36]. In our study, two alkaline phosphatase transcripts (Dv_149197_c0_seq1 and Dv_137932_c0_seq1) were also identified as DETs in both Cry34Ab1 and Cry34/35Ab1 treatments. The Dv_137932_c0_seq1 was down-regulated almost 3-fold (p < 0.05), suggesting that this alkaline phosphatase (Dv_137932_c0_seq1) is associated with toxin response.
The over expression of glutathione S-transferase (GST) has been associated with insecticide resistance in many insect pests [37]. Recent studies have also shown that the GST expression was reduced in Cry3Aa intoxicated T. molitor larvae [24] and Cry1Ab resistant Asian corn borer (Ostrinia furnacalis) [29]. In our study, one transcript, Dv_138610_c0_seq1 assigned as glutathione transferase was also down-regulated 1.6-fold (p < 0.05) in the Cry34/35Ab1 treatment. Among the four metabolism pathways identified in the Cry34Ab1 treatment, three of them (purine metabolism, starch and sucrose metabolism and amino sugar and nucleotide sugar metabolism) were common with pathways identified in the Cry34/35Ab1 treatment although 55 more metabolism transcripts were triggered in neonates challenged by Cry34/35Ab1 treatment. The higher number of DETs in the Cry34/35Ab1 may be indicative of higher levels of stress levels imposed by the toxic combination of the two toxins. These results are in agreement with findings by Li et al. [12], in which Cry34Ab1 enhances the binding of Cry35Ab1.
Several other pathways identified in the Cry34/35Ab1 treatment, such as nicotinate and nicotinamide metabolism, tryptophan metabolism, pyruvate metabolism and starch and sucrose metabolism, have been associated with environmental [38] and insecticide induced stress [15] in insects, indicating these pathways might be common response in insects to stressors.

4. Conclusions

NGS provides an important tool to investigate changes in gene expression associated with environmental challenges [24,39,40], and the combination of different statistical methods in downstream analysis of expression effects (DESeq, edgeR and limma) improved the different gene-calling results. Multiple transcripts were detected as responsive to the challenges of Cry34Ab1 and Cry34/35Ab1 combined. However, Cry35Ab1 alone did not produce a response, which is consistent with the lack of toxicity for this toxin. No previously identified Bt receptor genes (except for an alkaline phosphatase) were identified as differentially expressed, suggesting the receptor for Cry34Ab1 and Cry35Ab1 might be unique or that expression of the specific receptor is unaffected by toxin exposure.
It is also possible that actual toxin receptor(s) are not responsive to toxin exposure and that the observed differences in expression are a function of cellular stress and subsequent repair processes. Further analysis is needed to assess not only the function of those genes significantly affected by exposure to Cry34/35Ab1, but also explore the relative changes of associated proteins especially given the large percentage of DEGs that had no annotation. The data obtained herein should facilitate a better understanding of the active response to Cry toxin challenge at a transcriptomic level and provide new insights into the interaction of WCR and the Cry34/35Ab1 binary toxin.

5. Materials and Methods

5.1. Insects

The susceptible western corn rootworm eggs from a non-diapause WCR strain, which has been reared continuously for more than 30 years in the absence of insecticide and any Bt toxins exposure, were purchased from Crop Characteristics, Inc. (Farmington, MN, USA) and incubated in a growth chamber at 26 ± 1 °C and 60% ± 10% relative humidity with a photoperiod of 12:12 h (L:D) until hatching occurred approximately two weeks later.

5.2. Bt Proteins

Full-length Cry35Ab1 protein (44 kDa) was digested with chymotrypsin to generate active protein core fragments. Briefly, full-length Cry35Ab1 was incubated with bovine pancreatic chymotrypsin (Sigma, St. Louis, MO, USA) at 50:1 (w/w ratio = Cry protein:enzyme) in 100 mM sodium citrate buffer, pH 3.0, at 4 °C with gentle shaking for 2–3 days. The resulting core fragment (40 kDa) was analyzed on a 12% SDS-PAGE gel as described by Crespo et al. [41]. The activated Cry35Ab1 and full-length Cry34Ab1were used for all experiments in this study.

5.3. Exposure

Exposure experiments were conducted in 24 well cell culture plates (Costa 3526, Corning Incorporated, NY, USA). One mL of Dow AgroSciences proprietary corn rootworm diet was dispensed into wells of each plate and the surface was coated with Bt protein(s) at 15 µg/cm2 for Cry34Ab1 or Cry35Ab1 alone and 15 µg/cm2 Cry34Ab1 + 15 µg/cm2 Cry35Ab1 combined. The Cry toxins were diluted in 20 mM sodium acetate solution at pH 3.5. Controls consisted of wells treated with sodium acetate solution only. Six replicates for each treatment, including controls, were prepared for a total of 24 samples. Approximately 32 neonates (<24 h after hatching) were transferred into each pre-coated well with a fine camel hair paint brush and were exposed to Bt protein(s) for 48 h at room temperature. All living neonates in each well were pooled, snap frozen in liquid nitrogen and stored at −80 °C until RNA extraction.

5.4. RNA Isolation

Total RNA was extracted from pooled samples using the RNeasy Mini Kit (Cat. 74104, Qiagen, Germantown, MD, USA) and treated with RNase-Free Dnase (Cat. 79254, Qiagen) to eliminate DNA contamination according to the manufacturer’s instructions. The quality of RNA samples was evaluated on 1% agrose gels and quantity was estimated on NanoDrop-1000 (Thermo Fisher Scientific, Bartlesville, OK, USA) before submission for RNAseq analysis.

5.5. Next Generation Sequencing

The RNA sample integrity of all twenty-four RNA samples was further assessed using an Agilent 2100 Bioanalyzer (Cat. G2940CA, Agilent Technologies, Richardson, TX, USA) at the Next Generation Sequencing Core Facility at Durham Research Center, University of Nebraska Medical Center and were processed for library construction and paired end sequencing on an Illumina HiSeq2000 system (San Diego, CA, USA). All samples were sequenced using 100 bp paired-end reads.

5.6. Read Mapping and Differential Expression Analysis

A stringent quality filter process was applied by removing reads that did not have a minimum Phred quality score (Q64) of 20 per base corresponding to a 1% expected error rate using Sickle/1.2 (version 1.2, San Francisco, USA, [42]) according to the manual instructions. To map the quality reads back to WCR reference transcriptome contigs [20], the Bowtie aligner (2013, version 1.0.0, Baltimore, MD, USA, [43]), Samtools (version 1.3, La Jolla, CA, USA, [44]) was used to retrieve read counts for all treatments and control. All data analyses were performed at Holland Computing Center (HCC) at the University of Nebraska.
To identify differentially expressed genes among the treatments, three commonly used statistical methods (DESeq [45], edgeR [46] and limma [47]) for detecting differential expression in RNA-seq studies were employed [48]. For each treatment, only differentially expressed transcripts (DETs) that were identified by all three methods were used in further analysis. The edgeR package was initially used to remove low read counts from all 24 samples at a threshold of >1 cpm (count per million). All filtered data were used as input for differential expression analysis at adjusted p < 0.05 (hereafter p < 0.05). For edgeR and limma, TMM normalization was used to identify differentially expressed genes. For DESeq, the size factor normalization was used.

5.7. Homology Searches and Gene Ontology Analysis

The resulting DETs from differential expression analysis described above were annotated using the BLASTx algorithm against the NCBI nonredundant (nr) database with a stringent e-value threshold cut-off of 10−25. For gene ontology analysis, BLAST2GO [49] was employed for further analysis at default settings. With the annotated WCR transcriptome database as reference, the GO enrichment analyses were conducted using the Fisher’s Exact Test implemented in the functional enrichment feature of BLAST2GO at default settings (FDR < 0.05) [49].

5.8. Validation of Gene Expression via qPCR

Two up-regulated DEGs, GSC (gut-specific chitinase, 12.55-fold) and PAT (proton-coupled amino acid transporter 1-like, 8.69-fold) and two down-regulated DEGs, GH45 (endo-beta glucansase 32-fold change) and ALP (alkaline phosphatase, 7.4-fold) identified in both Cry34Ab1 and Cry34/35Ab1 treatments were selected for validation by qRT-PCR analysis. The qPCR primers for target genes were designed with a web-based tool (Primer3plus, 2012, Singapore, Singapore [50]). The actin gene was used as endogenouse gene (housekeeping gene) and the primer efficiency tests were conducted as described by Rodrigues et al. [51]. Primers with efficiency between 90 and 110% were selected for qPCR (Table S8). The cDNA used in this validation experiment were synthesized with a different batch of RNA prepared from different samples collected in Exposure step previously. The RNAs from three biological samples from Cry34/35Ab1 treatment and control (buffer) were used for cDNA synthesis with QuantiTect Rev Transcription Kit (Qiagen, Cat. 205311) according to the manufacture’s instruction. The 2−ΔΔCT method was used to calculate the relative expression of target genes [52].

Supplementary Materials

The following are available online at www.mdpi.com/2072-6651/9/4/124/s1, Table S1: Reads generated on RNAseq for each treatment and control, Table S2: Summary of illumina generated read production and read for mapping after filtering, Table S3: 18 unique functional activities associated with 34/35Ab1 treatment as compared to 34Ab1 treatment, Table S4: Summary of GO term enrichment results of significantly regulated genes in treatment with Cry34Ab1 at FDR < 0.05, Table S5: Summary of GO term enrichment results of significantly up-regulated DEGs in treatment with Cry34/35Ab1 at FDR < 0.05, Table S6: Summary of GO term enrichment results of significantly down-regulated DEGs in treatment with Cry34/35Ab1 at FDR < 0.05, Table S7: KEGG pathways involving 323 DEGs in Cry34/35Ab1 treatment as compared to the control, Table S8: Primers used in validation of gene expression via qPCR.

Acknowledgments

The funding for this research is fully provided by the Dow AgroSciences. We also want to thank Natalie Matz in the toxicology laboratory, Department of entomology at UNL for RNA preparation. The RNAseq data analysis was carried out at Holland Computer Center (HCC) located at the University of Nebraska.

Author Contributions

H.W. and B.S. conceived and designed the experiments; H.W. performed the experiments; H.W. and K.A. analyzed the data; S.Y.T., E.M., H.L., M.D., C.K. and K.E.N. contributed reagents/materials/analysis tools; H.W., B.S., K.E.N., S.E. and P.G. wrote the paper. J.J. helped in sample collection.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Levine, E.; Oloumisadeghi, H. Management of diabroticite rootworms in corn. Annu. Rev. Entomol. 1991, 36, 229–255. [Google Scholar] [CrossRef]
  2. Sappington, T.W.; Siegfried, B.D.; Guillemaud, T. Coordinated diabrotica genetics research: Accelerating progress on an urgent insect pest problem. Am. Entomol. 2006, 52, 90–97. [Google Scholar] [CrossRef]
  3. Miller, N.J.; Guillemaud, T.; Giordano, R.; Siegfried, B.D.; Gray, M.E.; Meinke, L.J.; Sappington, T.W. Genes, gene flow and adaptation of diabrotica virgifera virgifera. Agric. For. Entomol. 2009, 11, 47–60. [Google Scholar] [CrossRef]
  4. Metcalf, R.L. Methods for the study of pest diabroticai. In Methods for the Study of Pest Diabrotica; Krysan, J.L., Miller, T.A., Eds.; Springer: New York, NY, USA, 1986; pp. 7–15. [Google Scholar]
  5. Meinke, L.J.; Siegfried, B.D.; Wright, R.J.; Chandler, L.D. Adult susceptibility of nebraska western corn rootworm (coleoptera: Chrysomelidae) populations to selected insecticides. J. Econ. Entomol. 1998, 91, 594–600. [Google Scholar] [CrossRef]
  6. Parimi, S.; Meinke, L.J.; French, B.W.; Chandler, L.D.; Siegfried, B.D. Stability and persistence of aldrin and methyl-parathion resistance in western corn rootworm populations (coleoptera: Chrysomelidae). Crop Prot. 2006, 25, 269–274. [Google Scholar] [CrossRef]
  7. Levine, E.; Spencer, J.L.; Isard, S.A.; Onstad, D.W.; Gray, M.E. Adaptation of the western corn rootworm to crop rotation: Evolution of a new strain in response to a management practice. Am. Entomol. 2002, 48, 94–107. [Google Scholar] [CrossRef]
  8. Gray, M.E.; Sappington, T.W.; Miller, N.J.; Moeser, J.; Bohn, M.O. Adaptation and invasiveness of western corn rootworm: Intensifying research on a worsening pest. Annu. Rev. Entomol. 2009, 54, 303–321. [Google Scholar] [CrossRef] [PubMed]
  9. Gassmann, A.J.; Petzold-Maxwell, J.L.; Keweshan, R.S.; Dunbar, M.W. Field-evolved resistance to bt maize by western corn rootworm. PLoS ONE 2011. [Google Scholar] [CrossRef] [PubMed]
  10. Gassmann, A.J.; Petzold-Maxwell, J.L.; Clifton, E.H.; Dunbar, M.W.; Hoffmann, A.M.; Ingber, D.A.; Keweshan, R.S. Field-evolved resistance by western corn rootworm to multiple bacillus thuringiensis toxins in transgenic maize. Proc. Natl. Acad. Sci. USA 2014, 111, 5141–5146. [Google Scholar] [CrossRef] [PubMed]
  11. Wangila, D.S.; Gassmann, A.J.; Petzold-Maxwell, J.L.; French, B.W.; Meinke, L.J. Susceptibility of nebraska western corn rootworm (coleoptera: Chrysomelidae) populations to bt corn events. J. Econ. Entomol. 2015, 108, 742–751. [Google Scholar] [CrossRef] [PubMed]
  12. Li, H.R.; Olson, M.; Lin, G.F.; Hey, T.; Tan, S.Y.; Narva, K.E. Bacillus thuringiensis Cry34Ab1/Cry35Ab1 interactions with western corn rootworm midgut membrane binding sites. PLoS ONE 2013. [Google Scholar] [CrossRef] [PubMed]
  13. Moellenbeck, D.J.; Peters, M.L.; Bing, J.W.; Rouse, J.R.; Higgins, L.S.; Sims, L.; Nevshemal, T.; Marshall, L.; Ellis, R.T.; Bystrak, P.G.; et al. Insecticidal proteins from bacillus thuringiensis protect corn from corn rootworms. Nat. Biotechnol. 2001, 19, 668–672. [Google Scholar] [CrossRef] [PubMed]
  14. Mamidala, P.; Wijeratne, A.J.; Wijeratne, S.; Kornacker, K.; Sudhamalla, B.; Rivera-Vega, L.J.; Hoelmer, A.; Meulia, T.; Jones, S.C.; Mittapalli, O. RNA-seq and molecular docking reveal multi-level pesticide resistance in the bed bug. BMC Genom. 2012. [Google Scholar] [CrossRef] [PubMed]
  15. Lin, Q.S.; Jin, F.L.; Hu, Z.D.; Chen, H.Y.; Yin, F.; Li, Z.Y.; Dong, X.L.; Zhang, D.Y.; Ren, S.X.; Feng, X. Transcriptome analysis of chlorantraniliprole resistance development in the diamondback moth Plutella xylostella. PLoS ONE 2013. [Google Scholar] [CrossRef] [PubMed]
  16. Sparks, M.E.; Blackburn, M.B.; Kuhar, D.; Gundersen-Rindal, D.E. Transcriptome of the lymantria dispar (gypsy moth) larval midgut in response to infection by bacillus thuringiensis. PLoS ONE 2013. [Google Scholar] [CrossRef] [PubMed]
  17. Lei, Y.Y.; Zhu, X.; Xie, W.; Wu, Q.J.; Wang, S.L.; Guo, Z.J.; Xu, B.Y.; Li, X.C.; Zhou, X.G.; Zhang, Y.J. Midgut transcriptome response to a cry toxin in the diamondback moth, Plutella xylostella (lepidoptera: Plutellidae). Gene 2014, 533, 180–187. [Google Scholar] [CrossRef] [PubMed]
  18. Canton, P.E.; Cancino-Rodezno, A.; Gill, S.S.; Soberon, M.; Bravo, A. Transcriptional cellular responses in midgut tissue of aedes aegypti larvae following intoxication with cry11aa toxin from Bacillus thuringiensis. BMC Genom. 2015. [Google Scholar] [CrossRef] [PubMed]
  19. Gassmann, A.J.; Shrestha, R.B.; Jakka, S.R.; Dunbar, M.W.; Clifton, E.H.; Paolino, A.R.; Ingber, D.A.; French, B.W.; Masloski, K.E.; Dounda, J.W.; et al. 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, 1872–1880. [Google Scholar] [CrossRef] [PubMed]
  20. Eyun, S.I.; Wang, H.C.; Pauchet, Y.; Ffrench-Constant, R.H.; Benson, A.K.; Valencia-Jimenez, A.; Moriyama, E.N.; Siegfried, B.D. Molecular evolution of glycoside hydrolase genes in the western corn rootworm (Diabrotica virgifera virgifera). PLoS ONE 2014. [Google Scholar] [CrossRef] [PubMed]
  21. Schnepf, H.E.; Lee, S.; Dojillo, J.; Burmeister, P.; Fencil, K.; Morera, L.; Nygaard, L.; Narva, K.E.; Wolt, J.D. Characterization of Cry34/Cry35 binary insecticidal proteins from diverse bacillus thuringiensis strain collections. Appl. Environ. Microbiol. 2005, 71, 1765–1774. [Google Scholar] [CrossRef] [PubMed]
  22. Ellis, R.T.; Stockhoff, B.A.; Stamp, L.; Schnepf, H.E.; Schwab, G.E.; Knuth, M.; Russell, J.; Cardineau, G.A.; Narva, K.E. Novel bacillus thuringiensis binary insecticidal crystal proteins active on western corn rootworm, Diabrotica virgifera virgifera leconte. Appl. Environ. Microbiol. 2002, 68, 1137–1145. [Google Scholar] [CrossRef] [PubMed]
  23. Herman, R.A.; Scherer, P.N.; Young, D.L.; Mihaliak, C.A.; Meade, T.; Woodsworth, A.T.; Stockhoff, B.A.; Narva, K.E. Binary insecticidal crystal protein from bacillus thuringiensis, strain PS149B1: Effects of individual protein components and mixtures in laboratory bioassays. J. Econ. Entomol. 2002, 95, 635–639. [Google Scholar] [CrossRef] [PubMed]
  24. Oppert, B.; Dowd, S.E.; Bouffard, P.; Li, L.; Conesa, A.; Lorenzen, M.D.; Toutges, M.; Marshall, J.; Huestis, D.L.; Fabrick, J.; et al. Transcriptome profiling of the intoxication response of tenebrio molitor larvae to bacillus thuringiensis Cry3Aa protoxin. PLoS ONE 2012. [Google Scholar] [CrossRef] [PubMed]
  25. Derecka, K.; Blythe, M.J.; Malla, S.; Genereux, D.P.; Guffanti, A.; Pavan, P.; Moles, A.; Snart, C.; Ryder, T.; Ortori, C.A.; et al. Transient exposure to low levels of insecticide affects metabolic networks of honeybee larvae. PLoS ONE 2013. [Google Scholar] [CrossRef] [PubMed]
  26. Dou, W.; Shen, G.M.; Niu, J.Z.; Ding, T.B.; Wei, D.D.; Wang, J.J. Mining genes involved in insecticide resistance of liposcelis bostrychophila badonnel by transcriptome and expression profile analysis. PLoS ONE 2013. [Google Scholar] [CrossRef] [PubMed]
  27. Vellichirammal, N.N.; Wang, H.C.; Eyun, S.; Moriyama, E.N.; Coates, B.S.; Miller, N.J.; Siegfried, S.D. Transcriptional analysis of susceptible and resistant european corn borer strains and their response to Cry1f protoxin. BMC Genom. 2015. [Google Scholar] [CrossRef] [PubMed]
  28. Yang, W.T.; Dierking, K.; Esser, D.; Tholey, A.; Leippe, M.; Rosenstiel, P.; Schulenburg, H. Overlapping and unique signatures in the proteomic and transcriptomic responses of the nematode caenorhabditis elegans toward pathogenic bacillus thuringiensis. Dev. Comp. Immunol. 2015, 52, 1–9. [Google Scholar] [CrossRef] [PubMed]
  29. Xu, L.N.; Wang, Y.Q.; Wang, Z.Y.; Hu, B.J.; Ling, Y.H.; He, K.L. Transcriptome differences between Cry1Ab resistant and susceptible strains of asian corn borer. BMC Genom. 2015. [Google Scholar] [CrossRef] [PubMed]
  30. Sayed, A.; Wiechman, B.; Struewing, I.; Smith, M.; French, W.; Nielsen, C.; Bagley, M. Isolation of transcripts from diabrotica virgifera virgifera leconte responsive to the bacillus thuringiensis toxin Cry3Bb1. Insect Mol. Biol. 2010, 19, 381–389. [Google Scholar] [CrossRef] [PubMed]
  31. Song, P.; Wang, Q.Y.; Nangong, Z.Y.; Su, J.P.; Ge, D.H. Identification of henosepilachna vigintioctomaculata (coleoptera: Coccinellidae) midgut putative receptor for bacillus thuringiensis insecticidal Cry7Ab3 toxin. J. Invertebr. Pathol. 2012, 109, 318–322. [Google Scholar] [CrossRef] [PubMed]
  32. Ochoa-Campuzano, C.; Real, M.D.; Martinez-Ramirez, A.C.; Bravo, A.; Rausell, C. An adam metalloprotease is a Cry3Aa bacillus thuringiensis toxin receptor. Biochem. Biophys. Res. Commun. 2007, 362, 437–442. [Google Scholar] [CrossRef] [PubMed]
  33. 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, 123–128. [Google Scholar] [CrossRef] [PubMed]
  34. Bravo, A.; Gill, S.S.; Soberon, M. Mode of action of bacillus thuringiensis Cry and Cyt toxins and their potential for insect control. Toxicon 2007, 49, 423–435. [Google Scholar] [CrossRef] [PubMed]
  35. Guo, Z.J.; Kang, S.; Chen, D.F.; Wu, Q.J.; Wang, S.L.; Xie, W.; Zhu, X.; Baxter, S.W.; Zhou, X.G.; Jurat-Fuentes, J.L.; et al. Mapk signaling pathway alters expression of midgut alp and abcc genes and causes resistance to bacillus thuringiensis cry1ac toxin in diamondback moth. PLoS Genet. 2015. [Google Scholar] [CrossRef] [PubMed]
  36. Pigott, C.R.; Ellar, D.J. Role of receptors in bacillus thuringiensis crystal toxin activity. Microbiol. Mol. Biol. Res. 2007, 71, 255–281. [Google Scholar] [CrossRef] [PubMed]
  37. Enayati, A.A.; Ranson, H.; Hemingway, J. Insect glutathione transferases and insecticide resistance. Insect Mol. Biol. 2005, 14, 3–8. [Google Scholar] [CrossRef] [PubMed]
  38. Teets, N.M.; Peyton, J.T.; Colinet, H.; Renault, D.; Kelley, J.L.; Kawarasaki, Y.; Lee, R.E.; Denlinger, D.L. Gene expression changes governing extreme dehydration tolerance in an antarctic insect. Proc. Natl. Acad. Sci. USA 2012, 109, 20744–20749. [Google Scholar] [CrossRef] [PubMed]
  39. Harrison, P.W.; Wright, A.E.; Mank, J.E. The evolution of gene expression and the transcriptome-phenotype relationship. Semin. Cell Dev. Biol. 2012, 23, 222–229. [Google Scholar] [CrossRef] [PubMed]
  40. Flagel, L.E.; Bansal, R.; Kerstetter, R.A.; Chen, M.; Carroll, M.; Flannagan, R.; Clark, T.; Goldman, B.S.; Michel, A.P. Western corn rootworm (diabrotica virgifera virgifera) transcriptome assembly and genomic analysis of population structure. BMC Genom. 2014. [Google Scholar] [CrossRef] [PubMed]
  41. Crespo, A.L.; Spencer, T.A.; Nekl, E.; Pusztai-Carey, M.; Moar, W.J.; Siegfried, B.D. Comparison and validation of methods to quantify Cry1Ab toxin from bacillus thuringiensis for standardization of insect bioassays. Appl. Environ. Microbiol. 2008, 74, 130–135. [Google Scholar] [CrossRef] [PubMed]
  42. Joshi, N.A.; Fass, J.N. Sickle: A Sliding-Window, Adaptive, Quality-Based Trimming Tool for Fastq Files. Available online: https://github.com/najoshi/sickle (accessed on 29 March 2017).
  43. Langmead, B.; Trapnell, C.; Pop, M.; Salzberg, S.L. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009. [Google Scholar] [CrossRef] [PubMed]
  44. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R.; 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and samtools. Bioinformatics 2009, 25, 2078–2079. [Google Scholar] [CrossRef] [PubMed]
  45. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Genome Biol. 2010. [Google Scholar] [CrossRef] [PubMed]
  46. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. Edger: A bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed]
  47. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Res. 2015. [Google Scholar] [CrossRef] [PubMed]
  48. Seyednasrollah, F.; Laiho, A.; Elo, L.L. Comparison of software packages for detecting differential expression in rna-seq studies. Brief. Bioinform. 2015, 16, 59–70. [Google Scholar] [CrossRef] [PubMed]
  49. Conesa, A.; Gotz, S.; Garcia-Gomez, J.M.; Terol, J.; Talon, M.; Robles, M. Blast2go: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21, 3674–3676. [Google Scholar] [CrossRef] [PubMed]
  50. Untergasser, A.; Nijveen, H.; Rao, X.; Bisseling, T.; Geurts, R.; Leunissen, J.A.M. Primer3plus, an enhanced web interface to primer3. Nucleic Acids Res. 2007, 35, W71–W74. [Google Scholar] [CrossRef] [PubMed]
  51. Rodrigues, T.B.; Khajuria, C.; Wang, H.C.; Matz, N.; Cardoso, D.C.; Valicente, F.H.; Zhou, X.G.; Siegfried, B. Validation of reference housekeeping genes for gene expression studies in western corn rootworm (Diabrotica virgifera virgifera). PLoS ONE 2014. [Google Scholar] [CrossRef] [PubMed]
  52. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T) (-delta delta C) method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
Figure 1. MA plots of differential expression in each treatment comparison with DESeq method. (a) buffer vs. Cry35Ab1; (b) buffer vs. Cry34Ab1; (c) buffer vs. Cry34/Cry35Ab1. Red dots represent the genes either up- or down-regulated at p (adj) < 0.05. MA plot is a way to display deferentially expressed genes versus expression strength (log2 fold change) between control and treatment.
Figure 1. MA plots of differential expression in each treatment comparison with DESeq method. (a) buffer vs. Cry35Ab1; (b) buffer vs. Cry34Ab1; (c) buffer vs. Cry34/Cry35Ab1. Red dots represent the genes either up- or down-regulated at p (adj) < 0.05. MA plot is a way to display deferentially expressed genes versus expression strength (log2 fold change) between control and treatment.
Toxins 09 00124 g001
Figure 2. Blast hit of DGEs. (a) top BLAST hits of DGEs in WCR neonates exposed to Cry34/35Ab1 combination; (b) top BLAST hits of DGEs in WCR neonates exposed to Cry34Ab1 alone.
Figure 2. Blast hit of DGEs. (a) top BLAST hits of DGEs in WCR neonates exposed to Cry34/35Ab1 combination; (b) top BLAST hits of DGEs in WCR neonates exposed to Cry34Ab1 alone.
Toxins 09 00124 g002
Figure 3. GO term categorization and distribution in Cry34Ab1 compared to control at level 2 under three main categories, i.e., biological process, molecular function and cellular component.
Figure 3. GO term categorization and distribution in Cry34Ab1 compared to control at level 2 under three main categories, i.e., biological process, molecular function and cellular component.
Toxins 09 00124 g003
Figure 4. GO term categorization and distribution in Cry34/35Ab1 vs. buffer control at level 2 under three main categories, i.e. biological process, molecular function and cellular component.
Figure 4. GO term categorization and distribution in Cry34/35Ab1 vs. buffer control at level 2 under three main categories, i.e. biological process, molecular function and cellular component.
Toxins 09 00124 g004
Table 1. RNAseq differential gene expression for WCR exposed to Cry34/35Ab1 toxin.
Table 1. RNAseq differential gene expression for WCR exposed to Cry34/35Ab1 toxin.
Analysis MethodCry34Ab1 vs. BufferCry35Ab1 vs. BufferCry34Ab1 + Cry35Ab1 vs. Buffer
updowntotalShared *updowntotalupdowntotalshared
DESeq4472116114000992122322151300
edgeR498313200064710261673
limma4887135000109312432336
* indicate the total DETs shared among three analysis method.
Table 2. Number of DETs shared and not shared between Cry34Ab and Cry34Ab1 + Cry35Ab1 treatment.
Table 2. Number of DETs shared and not shared between Cry34Ab and Cry34Ab1 + Cry35Ab1 treatment.
CategoryCategory# ContigsDigital Expression in Fold Change *
up-regulatedshared312.55 (2–4.66) a4.5 (2.3–12.5) b
unique12 a534 b2.36 (2.14–2.9) a2.69 (2–64) b
down-regulatedShared612.64 (2–5.6) a10.6 (2.49–194) b
unique10 a674 b3.16 (2–1) a3.3 (2–88) b
a challenged with Cry34Ab; b challenged with combination of Cry34/35Ab1; * all fold changes are from DESeq data at p < 0.05).
Table 3. Comparison of expression of four randomly selected genes with two analysis methods.
Table 3. Comparison of expression of four randomly selected genes with two analysis methods.
CategoryGeneFold in RNAseq AnalysisFold Change in qPCR
up-regulatedGSC12.55 *91.01 *
PAT8.69 *3.44 *
down-regulatedGH4532 *0.31 *
ALP7.4 *0.16 *
* The fold change either in RNAseq analysis or qPCR was from Cry34/35-treated samples and was compared to the expression of the same gene from the control; Abbreviation: GSC-gut-specific chitinase, PAT-proton-coupled amino acid transporter 1-like, GH45-endo-beta glucansase and ALP-alkaline phosphatase.

Share and Cite

MDPI and ACS Style

Wang, H.; Eyun, S.-i.; Arora, K.; Tan, S.Y.; Gandra, P.; Moriyama, E.; Khajuria, C.; Jurzenski, J.; Li, H.; Donahue, M.; 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, 124. https://doi.org/10.3390/toxins9040124

AMA Style

Wang H, Eyun S-i, Arora K, Tan SY, Gandra P, Moriyama E, Khajuria C, Jurzenski J, Li H, Donahue M, 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. https://doi.org/10.3390/toxins9040124

Chicago/Turabian Style

Wang, Haichuan, Seong-il Eyun, Kanika Arora, Sek Yee Tan, Premchand Gandra, Etsuko Moriyama, Chitvan Khajuria, Jessica Jurzenski, Huarong Li, Maia Donahue, and et al. 2017. "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 9, no. 4: 124. https://doi.org/10.3390/toxins9040124

APA Style

Wang, H., Eyun, S. -i., Arora, K., Tan, S. Y., Gandra, P., Moriyama, E., Khajuria, C., Jurzenski, J., Li, H., Donahue, M., Narva, K., & Siegfried, B. (2017). 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, 9(4), 124. https://doi.org/10.3390/toxins9040124

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop