Differential gene expression in response to eCry3.1Ab ingestion in an unselected and eCry3.1Ab-selected western corn rootworm (Diabrotica virgifera virgifera LeConte) population

Diabrotica virgifera virgifera LeConte, the western corn rootworm (WCR) is one of the most destructive pests in the U.S. Corn Belt. Transgenic maize lines expressing various Cry toxins from Bacillus thuringiensis have been adopted as a management strategy. However, resistance to many Bt toxins has occurred. To investigate the mechanisms of Bt resistance we carried out RNA-seq using Illumina sequencing technology on resistant, eCry3.1Ab-selected and susceptible, unselected, whole WCR neonates which fed on seedling maize with and without eCry3.1Ab for 12 and 24 hours. In a parallel experiment RNA-seq experiments were conducted when only the midgut of neonate WCR was evaluated from the same treatments. After de novo transcriptome assembly we identified differentially expressed genes (DEGs). Results from the assemblies and annotation indicate that WCR neonates from the eCry3.1Ab-selected resistant colony expressed a small number of up and down-regulated genes following Bt intoxication. In contrast, unselected susceptible WCR neonates expressed a large number of up and down-regulated transcripts in response to intoxication. Annotation and pathway analysis of DEGs between susceptible and resistant whole WCR and their midgut tissue revealed genes associated with cell membrane, immune response, detoxification, and potential Bt receptors which are likely related to eCry3.1Ab resistance. This research provides a framework to study the toxicology of Bt toxins and mechanism of resistance in WCR, an economically important coleopteran pest species.

The Bt (Bacillus thuringiensis) Cry3 δ-endotoxin superfamily is known for its specificity to coleopteran species 1,2 . Transgenic maize hybrids expressing Cry3Bb1, mCry3A, eCry3.1Ab as well as Cry34/35Ab1 have been introduced to the market to manage western corn rootworm (WCR, Diabrotica virgifera virgifera LeConte). As a highly adaptive species, WCR has developed resistance to broadcast soil insecticides 3 , aerial spray insecticides for reducing adult number 4,5 , crop rotation 6 , and Cry3 Bt toxins. Laboratory selection experiments have developed WCR colonies resistant to Cry3Bb1, mCry3A, Cry34/35Ab1 and eCry3.1Ab by continuously rearing WCR on transgenic maize lines [7][8][9][10] . Resistance to Cry3 Bt toxins has become a practical issue since field resistant WCR populations have been reported in many locations [11][12][13][14] . Although eCry3.1Ab is the most recent Bt toxin 15 to the market without reported control failure in the field, laboratory selection 9 and cross-resistance experiments 13 indicate that resistance to eCry3.1Ab is likely in the field due to the resistance to other Cry3 proteins 16 .
Three Cry toxins (Cry3Bb1, mCry3A and eCry3.1Ab) for WCR control are derived from the Cry3 superfamily. In lepidopteran systems, the consumption of Bt protoxin by larvae is followed by activation via cleavage by midgut

Results and Discussion
Transcriptome assembly and annotation of WCR transcriptome. The transcriptome of whole larvae and midgut were separately de novo assembled. Reads from either whole larvae or midgut were individually pooled to increase the assembly coverage. Low quality and mitochondrial sequences were removed prior to assembly. Transcriptome assembly resulted in a whole larval transcriptome with 204,842 contigs from 57 Gb of reads and midgut transcriptome with 226,115 contigs from 137 Gb of reads. The two transcriptomes had comparable average sequence length (measured as N50, a weighted median statistic that more than half of the nucleotides of a transcriptome belong to the contigs of this size N50 or longer), GC content and length distribution ( Table 1, Fig. 1). After removing duplicate contigs with more than 95% sequence similarity, we obtained 187,570 and 209,167 contigs respectively from larval and midgut transcriptome (Table 1). Herein, we refer to the two reduced redundancy sets of contigs as "unigene" sets. BUSCO 37 was applied to evaluate the completeness of each unigene set. All of the transcriptome and unigene sets cover over 95% of insecta single-copy orthologs. Both unigene sets maintain the coverage and integrity while duplication is reduced, indicating the acceptable qualities for functional analysis (Table 2).
Transcriptomes and their unigene sets were aligned to the NCBI non-redundant protein database (NR) using BLASTX. Only 42% of larval and 38.2% of midgut unigenes had significant BLASTX results. The species distribution of BLASTX top hits indicated that a predominant number of unigenes could be annotated by the Coleopteran model species Tribolium castaneum (Herbst) ( Table 3, Fig. 2). We also identified 372 and 400 unigenes from larval and midgut transcriptome, respectively, that were directly aligned to 13 known genes from Diabrotica species (Tables 3 and 4).

Differential expression analysis of WCR transcriptome.
To identify the genes involved in eCry3.1Ab response, we analyzed the differentially expressed unigenes from whole larval and midgut transcriptomes between eCry3.1Ab-feeding and non-Bt isoline feeding WCR at both 12-and 24-hour time points. To understand the expression differences between eCry3.1Ab-selected and susceptible, control WCR from the same original population, we also compared the expression differences in whole larvae and midgut between two populations when being fed with the same kind of maize root. After the alignment and filtering out unigenes with extremely low expression levels, only 31,875 of larval and 22,954 of midgut transcriptome unigenes were proceeded to edgeR analysis 38 . The patterns of differentially expressed genes are shown in Fig. 3. Regardless of exposure time or tissue type, susceptible WCR differentially expressed a much larger number of genes in response to eCry3.1Ab intoxication. In contrast, eCry3.1-selected WCR had dramatically fewer differentially expressed genes (DEGs). The selected and unselected WCR shares many DEGs in both whole larvae and midgut transcriptome, while some unique DEGs were colony specific, especially in unselected colony (Fig. 4).
Albeit smaller in size, we identified DEGs between eCry3.1Ab-selected and unselected WCR, especially when both groups of neonates were fed on eCry3.1Ab maize root. The function of those DEGs and the pathways in which they were involved may reveal the physiological differences and the mechanism of eCry3.1Ab resistance in the selected population.
www.nature.com/scientificreports www.nature.com/scientificreports/ GO annotation and pathway analysis on eCry3.1Ab feeding WCR midgut. To further investigate the molecular and physiological adaptation to intoxication on eCry3.1Ab-selected WCR we applied gene ontology (GO) analysis to differentially expressed genes using Blast2GO. We compared the top 20 level-2 GO terms called from the unigene set from midgut DEGs between selected and unselected WCR feeding on eCry3.1Ab maize for 24 h. The GO terms were ranked based on number of unigenes of each GO term. The terms "metabolic  Table 1. Summary statistics of WCR whole larvae and midgut transcriptome assemblies and their unigene sets.   www.nature.com/scientificreports www.nature.com/scientificreports/ process", "catalytic activity" and "membrane" accounted for the majority of each of the three ontologies [Biological Process (BP), Molecular Function (MF) and Cellular Component (CC), respectively], suggesting the primary location and functions of eCry3.1Ab resistance. Other than that, the BP term "cellular process" and MF term "binding" also implied the binding and processing of toxin may have a role in eCry3.1Ab sensitivity.  Table 3. Summary of BLASTX of WCR whole larval and midgut transcriptome assemblies and their unigene sets.   www.nature.com/scientificreports www.nature.com/scientificreports/ We compared the distribution of level-2 GO terms between the DEG unigene set described above and midgut transcriptome unigene set ( Table 5). The BP terms "Multiple-organism process", "immune system process", "detoxification" and "cell killing" and the CC terms "extracellular region" and "extracellular part" were among the most over represented GOs in DEG unigene sets, while the BP term "reproduction" and MF term "structural molecule activity" were the most under represented. However gene set enrichment analysis (GSEA) showed that no specific GO term was significantly enriched in the DEG unigene set.
We further predicted enzyme codes (EC) from GO terms for unigenes and used the EC to map differentially expressed unigenes to pathways from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database 39 . The differentially expressed midgut unigenes were involved in 49 pathways, including purine metabolism (KEGG ID: 00230), glutathione metabolism (KEGG ID: 00480), fatty acid synthesis (KEGG ID: 00061), glycerophospholipid metabolism (KEGG ID: 00564) and drug metabolism (KEGG ID: 00983). The fatty acid synthesis pathway was also identified in whole larval differentially expressed unigenes. These results suggest that differences in eCry3.1Ab tolerances might arise from the alternation of genes related to detoxification, membrane functions and metabolism.
Expression of potential and novel eCry3.1Ab resistant genes. Research has revealed that cadherins 26 , ABC transporters 40 , and aminopeptidase N (APN) 41 are Bt receptors in lepidopteran insects. Both cadherins and ABC transporters have been found in WCR 22,23 . Upon BLASTX, we found 102 unigenes with cadherin function, 209 unigenes with ABC transporter function, and 50 unigenes with APN function from the midgut unigene set. In whole larval unigene set, there were 78 unigenes with cadherin function, 191 unigenes with ABC transporter function, and 62 unigenes with APN function.
The midgut gene expression profile showed that cadherin was not significantly differentially expressed, while some ABC transporters (multidrug resistance-associated protein) and APNs were differentially expressed following Cry3.1Ab ingestion, especially in susceptible WCR (Table 6, Fig. 5). The expression level of one APN midgut unigene (comp127589_c1_seq. 1) was significantly higher in susceptible WCR midgut, when both selected and susceptible-WCR were given eCry3.1Ab for 24 hours. Two ABC transporter unigenes (comp121889_c0_seq. 1, comp126268_c0_seq. 1) showed the same pattern, but the increased expression in susceptible WCR was not significant. The protein structure of eCry3.1Ab toxin contains partial Cry1Ab sequence at its C-terminus 15,42 . In multiple lepidopteran species APNs serve as binding receptors of Cry1A toxins and mutations or reduced expression of APN result in resistance to Bt toxins 41 . Our results suggest that APN is a potential eCry3.1Ab target and the reduced expression level under intoxication might contribute to resistance in WCR.
Proteases have been associated with Bt resistance either by increasing digestion of toxins 43,44 , or by decreasing the proteolytic activation of Bt pro-toxins 45 . We identified 112 unigenes with metalloprotease functions and 199 unigenes with cathepsin functions in the larval unigene set. In the midgut unigene set, there were 131 unigenes with metalloproteases functions and 234 unigenes with cathepsin functions. Expression levels of some proteases were regulated by Bt intoxication either in selected or susceptible WCR, while none of those digestive proteases was up regulated in selected WCR midgut compared to the susceptible one when both insects were fed with eCry3.1Ab maize. www.nature.com/scientificreports www.nature.com/scientificreports/ Considering the expression patterns, GO annotation, and pathway analysis, we infer that at least two novel genes are likely involved in the resistance to eCry3.1Ab. Esterase has been reported involved in Cry1Ac resistance in Helicoverpa armigera 46 . In WCR, we observed a novel esterase (comp36305_c0_seq. 1) expressed in the midgut of selected colony was higher than unselected one after 24 hours of intoxication. The second gene is a dynein heavy chain-like protein (comp127369_c3_seq. 14). It was constitutively up regulated in selected WCR regardless of diet and time. Since dynein is a cytoskeletal motor protein involved in intracellular transportation and the movement of chromosomes, we propose that selected-WCR many have a stronger activity in either endocytosis or cell mitosis to remove the attached eCry3.1Ab molecules 47 , or to repair damaged epithelial cells 48 .

Conclusions
After a comprehensive analysis of an RNA-seq experiment with eCry3.1Ab-selected and susceptible, control WCR, the transcriptome, unigene sets and reads provided numerous resources for studying the interaction between Cry3 and this coleopteran species. This study is the first step approaching the delineation of Bt resistance mechanisms in WCR. We propose more than one potential mechanism of resistance to eCry3.1Ab -a dual action Bt toxin. With the recently published WCR genome sequence, future research will detect the genomic-wide genetic variations associated with Bt resistance. Studies to explore how Bt toxins affect gene alternative splicing and whether the alternatively sliced genes are related to Bt resistance in WCR. We are also developing cellular and molecular methods i.e. cultured cell expression, RNAi gene silencing, individual genotyping to further study the detailed mechanism of resistance to eCry3.1Ab as well as other Bt toxins. Continuous discoveries in this field will lead to improving strategies for insect resistant management and the developing of novel entomotoxins.

Material and Methods
Insects and bioassay. The eCry3.1Ab-selected resistant WCR colony was initially selected and reared on nonelite non-commercial eCry3.1Ab-expressing transgenic maize (event 5130) under laboratory conditions 9 . Both the selected and susceptible control colonies were developed from a single population and had been maintained on eCry3.1Ab-expressing transgenic maize (material ID 12MG00345) and its near-isoline (material ID 12MG001181), respectively, for more than 30 generations. For the current experimental design (Supplementary Table 1, column 1), both eCry3.1Ab-expressing and isoline maize seeds were surface sterilized and germinated in Petri dishes with moistened filter paper at 23 °C for 4-6 days without illumination. Approximately 30 neonates hatched within 24 hours were transferred to a Petri dish containing 3-4 seedlings of each line. After 12 or 24 hours feeding, the living first-instar larvae were recovered. In a separate identical experiment (Supplementary Table 1, column 2), the midgut was dissected from 30-40 recovered larvae of each Petri dish. Both whole larvae and midgut were flash frozen in liquid nitrogen and stored at −80 °C until processing. Both the whole larvae and midgut bioassays were replicated independently three times as full biological replications.  www.nature.com/scientificreports www.nature.com/scientificreports/ manufacturer's instructions. For whole larvae and midgut respectively, 24 libraries (2 insect colonies × 2 corn lines × 2 times × 3 biological replications) were normalized, pooled and sequenced on two lanes of Illumina HiSeq2000 sequencer using 100-nucleotide pair-end protocol (Global Biologics LLC, Columbia, MO, USA).
De novo assembly of transcriptome. Adapters and low quality reads were trimmed using FastqMcf (version 1.04.803, https://github.com/ExpressionAnalysis/ea-utils/blob/wiki/FastqMcf.md) and Trimmomatic (version 0.36) 49 , respectively. To remove the mtDNA, the reads were aligned to the WCR mitochondrial genome 50 using Bowtie 2 (version 4.7.7) 51 with default settings. Paired reads concordantly aligned with no mismatch were considered as mitochondrial reads and discarded. These steps resulted in "clean reads" for assembly and differential expression analysis.
Trinity (version r2013-11-10) 52 was used to de novo assemble larval and midgut transcriptomes using cleaned reads pooled from libraries of each sample type with default setting. The unigene set of each transcriptome was obtained by removing sequence with over 95% of similarity in Blast2GO 53 (v4.1.9). BUSCO analyses (version v3 with Insecta odb9 dataset) 37 were performed on transcriptome assemblies and unigene sets to evaluate their quality and completeness. Differential expression analysis. Cleaned reads were aligned to the corresponding transcriptome unigene set using Bowtie 2 with default pair-end settings 51 . The output SAM files were converted to BAM format using SAMtools (version 0.1.20) 54 . The differential expression analysis was conducted in R (version 3.4.1) 55 . The read counts were called in R using the GenomicAlignment (version 1.6.3) and GenomicRanges (1.22.4) packages 56 . We counted only concordant alignment pairs while accepting multiple mapping reads due to the potential existence of isoforms in the unigene sets. The low expression unigenes were removed by applying filters with at least 2 count per million (CPM) over 3 samples. The differentially expressed contigs were assessed using the edgeR-robust algorithm of the edgeR package (version 3.12.1) 38,57 with the trimmed mean of M-values (TMM) normalization method 58,59 . False discovery rate (FDR) was controlled at 0.05 by the edgeR package and was used to determine the significance of differentially expressed genes (DEGs).  Table 6. Description of treatment groups for WCR larvae shown in Fig. 5. Figure 5. Expression levels of 11 potential eCry3.1Ab resistance related genes in 8 WCR midgut treatment group (see Table 6). Expression level s are quantified by count per million read (cpm). The candidate genes are categorized into: 1: potential Bt receptors; 2: digestive proteases; 3: detoxification enzymes; 4: enzymes involve in drug metabolism pathway; 5: enzymes involve in membrane related pathways; 6: other candidates.