Comparative Transcriptome Analysis of Bt Resistant and Susceptible Strains in Ostrinia furnacalis (Guenée) (Lepidoptera: Crambidae)

The evolution of target pest population resistance to Bt toxins is the most relevant threat to the sustainability of Bt technology, thus it is necessary to clarify insect resistance mechanisms. Firstly, the resistance level of Asian corn borer was determined by bioassay. After 28 generations selection in the lab, the Cry1Ie-resistant strain (ACB-IeR) developed more than 862-fold resistance to Cry1Ie, and the Cry1F-resistant strain (ACB-FR) developed 961-fold resistance to Cry1F. The results show that long-term exposure to Bt toxins can lead to resistance. Then, we compared the differential expression genes (DEGs) of ACB-FR and ACB-IeR with susceptible strain (ACB-BtS), and analyzed GO function and KEGG pathway through transcriptome sequencing. The comparison showed that in Bt-resistant strains, many genes have a significant down-regulated trend. Several Bt-resistance candidate genes were differentially expressed in both resistant strains. Furthermore, the DEGs were verified by RT-qPCR and showed similar trend. These results provide candidate genes for further research on the Bt resistance mechanism.


Introduction
Corn (Zea mays L.), as one of the most important food crops, is a high-quality feed crop and a high-yield cash crop, as well as an industrial raw material [1]. Asian corn borer (ACB),

Bioassay and Resistance Ratio
Three O. furnacalis laboratory strains were raised in State Key Laboratory for the Biology of the Plant Diseases and Insect Pests, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China. Cry1Ie-resistant strains (ACB-IeR) are exposed to Cry1Ie toxins with a concentration of 16.0 µg/g and Cry1F-resistant strains (ACB-FR) are exposed to Cry1F toxins with a concentration of 3.0 µg/g to maintain the intensity of selection, and susceptible strain (ACB-BtS) is kept on an artificial diet without any toxins. The Cry1F toxin protein used in this experiment was produced by Marianne P. Carey, Case Western Reserve University, USA. Cry1Ie toxin protein was purchased from ABZYMO Biosciences, Beijing, China.
The protein-food feeding method was used to determine the toxicological reaction of each strain to Bt toxin [22]. ACB-BtS was exposed to six concentrations of Cry1Ie toxin (0, 0.2, 0.4, 1, 2, 5 µg/g) and seven concentrations of Cry1F toxin (0, 0.04, 0.08, 0.20, 0.40, 1.00, 2.00 5.00 µg/g); ACB-IeR was exposed to six concentrations of Cry1Ie toxin (0, 40, 80, 200, 400, 1000 µg/g); ACB-FR was exposed to six concentrations of Cry1F toxin (0, 10, 20, 50, 100, 250, 500 µg/g), and the corresponding concentration of artificial feed was distributed into 48-well plates. A single starvation-treated 24-h-old neonate was put into each hole; a total of 144 larvae were tested at each concentration. ddH 2 O replaced the toxins to mix with an artificial diet that was used as a control. The larvae were reared at 27 ± 1 • C, 70-80% relative humidity (RH) and a photoperiod of 16:8 h (L: D). Investigating the results after seven days, the number of survivals and total body weight were recorded, and larvae were considered dead if they had died or weighed less than 0.1 mg. The Probit statistical model in SPSS 19.0 was used to calculate the lethal concentration of each strain and each toxin that causes 50% of the individual's mortality (LC 50 ), as well as the LC 50 ratio of the susceptible strain and the resistant strain when exposed to the same toxin as the resistance ratio (RR).

Midgut Dissection and RNA Extraction
Fourth instar larvae were put on ice to make them comatose, the midgut was dissected and put into a PBS buffer to remove the food content, then transferred into water to wash it again, the water absorbed with filter paper and then put into 200 µL Trizol. The midguts of five individuals were collected as one biological replicate, and three biological replicates were collected from ACB-BtS, ACB-FR and ACB-IeR strains.
Total RNA was extracted from the midgut of ACB-BtS, ACB-FR and ACB-IeR strains following the manufacturer's instructions for Trizol reagent. The quantity and quality of RNA were evaluated by 1% agarose gel electrophoresis and spectrophotometry on Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA). All samples were stored at -80 • C until they were used for first-strand cDNA synthesis.

Library Preparation for Transcriptome Sequencing
Nine single-ended cDNA libraries were constructed using total RNA samples from the midgut of ACB-BtS, ACB-FR and ACB-IeR. Total RNA was extracted using Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's protocol. RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase-free agarose gel electrophoresis. After total RNA was extracted, mRNA of each sample was enriched by Oligo (dT) magnetic beads. Then, the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcription into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP and buffer. Then, the cDNA fragments were purified with QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, single A base added, and ligated to Illumina sequencing adapters. The ligation products were size-selected by agarose gel electrophoresis, PCR amplified and sequenced using Illumina Novaseq 6000 by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China).

Assembly and Functional Gene Annotation
By filtering raw data, lower-quality reads were filtered. The linker sequence was removed to obtain high-quality clean data by fastp (version 0.18.0) [23]. After screening the raw data, the Trinity program was used to assemble high-quality clean data [24]. To minimize the sequence's redundancy, the CD-HIT program was used to obtain the longest transcript contigs [25]. After the library quality has been evaluated to meet the requirements of bioinformatics analysis, expression level and gene structure analysis can be carried out. DESeq2 was used to calculate the expression level of candidate genes in different treatments [26]. BLASTX was used to compare unigene sequences with NCBI protein databases: set e-value ≤ 10 −5 , and using the annotation information with the highest similarity in the comparison as the final comment information. DEGs were calculated as FPKM (fragments per kilobase pair of exon model per million fragments mapped), which were used to compare the expression of up-regulated or down-regulated transcripts with ACB-BtS strain.

Selection of Differentially Expressed Genes among ACB-BtS, ACB-FR and ACB-IeR
Blast2GO was used to predict the classification of the function of unigenes in different treatments by EuKaryotic of orthologous groups (KOG) database [27]. KEGG pathway enrichment was conducted to analysis by the online KEGG Automatic Annotation Server (KAAS) [28,29]. In order to understand the differential gene, we not only performed GO function and KEGG pathway enrichment analysis for all the differential genes in the transcriptome, but also used the false discovery rate (FDR) method to determine the threshold. The FDR ≤ 0.05 and the absolute value of log2FC ≥ 1 were used as evaluation criteria to compare the significance level of gene expression between different ACB strains. In addition, the genes annotated as candidate Bt resistance genes were verified by RT-qPCR in O. furnacalis midgut.

RT-qPCR Analysis
The genes associated with Bt resistance were screened out from the transcriptome and summarized pathway annotation candidate genes in the KEGG pathway with p-value < 0.05, and the part of the screened resistance-related genes that was duplicated with genes involved in the KEGG pathway was deleted, then RT-qPCR was used to verify these genes. Total RNA was prepared from the midgut of ACB-FR, ACB-IeR and ACB-BtS strains with three biological replicates. The TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, China) was used to synthesize cDNAs. The primers used for RT-qPCR were designed with Primer 3PLUS (Supplementary Materials: Table S1).
Reactions were set up in triplicate using 20 µL reaction volumes containing 10 µL of 2× PerfectStart Green qPCR SuperMix (TransGen Biotech, China), 0.8 µL of PCR forward primer (10 µM), 0.8 µL of PCR reverse primer (10 µM), 0.4 µL of passive reference dye (50×), 1.0 µL of cDNA as template and 7.0 µL sterilized ultra-pure grade water. The ABI 7500 fast real-time PCR system was used for amplification at 95 • C (30 s), followed by 40 cycles at 95 • C (10 s) and 60 • C (20 s). Then, melting curve analysis was conducted to evaluate the specificity of amplification. Before RT-qPCR, the amplification efficiency of all the primers was verified. The actin gene of ACB was used to normalize the RT-qPCR data. Relative expression levels for target genes in different strains were calculated via the 2 −∆∆CT method [30]. The expression difference among ACB-FR, ACB-IeR and ACB-BtS strains was analyzed by one-way analyses of variance (ANOVAs) statistical method. Differences in the three strains were compared using Turkey test.

Insect Resistance Levels
After 28 generations of selection in the lab, ACB-IeR had developed more than 862-fold resistance to Cry1Ie (Table 1). Analogously, ACB-FR had developed 961-fold resistance to Cry1F (Table 1).

RNA-Seq and Sequence Assembly
Quality control of data was obtained from transcriptome sequencing, and the RNA sequencing results of ACB-FR, ACB-IeR and ACB-BtS ranged from 39,687,801 to 50,484,709 ( Table 2). The clean sequence of each library ranges from 39,127,741 to 49,893,262 reads. In addition, the GC content ranges from 48.99% to 49.95%. The total number of reads that can be mapped to the genome ranges from 31,285,151 to 41,890,268 and the proportion of effective reads was from 79.96% to 85.34%. The raw reads of ACB-BtS, ACB-FR and ACB-IeR strain were submitted to the SRA database with accession number SRR16556445, SRR16571663 and SRR16572026.

Differentially Expressed Genes between Cry1F and Cry1Ie-Resistant Strains Compared with Susceptible Strains of O. furnacalis
Differentially expression genes (DEGs) were screened with FDR ≤ 0.05 and log2FC absolute value > 1 as a significantly different gene (Supplementary Materials: Table S5). The results showed that a total of 283 down-regulated and 172 up-regulated genes were differentially expressed ( Figure 1A) in ACB-FR compared with ACB-BtS strains. Moreover, a total of 229 down-regulated and 133 up-regulated genes were selected in ACB-IeR compared with ACB-BtS strains ( Figure 1B). This comparison showed that a larger proportion of genes were significantly down-regulated in the Bt-resistant strain. Based on former reports and DEGs in the transcriptome analysis, some well-known Bt receptors such as cadherin, ABC transporter, ALP and APN were screened out and listed in Supplementary Table S2. Other genes include UDP-glucuronosyltransferase gene, trypsin, heat shock protein (HSP) 68-like gene, glutathione S-transferase 1-like gene (GSTs), etc. Several of these Bt-resistance candidate genes were differentially regulated in both of the resistant strains. Most of the heat shock proteins (Supplementary Materials: Table S6) were down-regulated in ACB-FR compared with ACB-BtS strains, but it was expressed higher in ACB-IeR. UDP-glucuronosyltransferase (Supplementary Materials: Table S7) and trypsin, especially alkaline C-like (Supplementary Materials: Table S8), were up-regulated in ACB-IeR and ACB-FR compared with ACB-BtS strains. The differential expression of GSTs in Cry1F-and Cry1Ie-resistant strains was shown in Table S9. From the list, five GSTs were up-regulated in ACB-FR and ACB-IeR, only ncbi_114360773 significantly up-regulated with log2FC value of 2.75 in ACB-FR strain and MSTRG.22972 up-regulated with log2FC value of 1.91 in the ACB-IeR strain. Three GSTs were down-regulated in ACB-FR and ACB-IeR strains, but only ncbi_114365341 was significantly down-regulated in ACB-IeR strains with log2FC value −1.66.
up-regulated in ACB-IeR and ACB-FR compared with ACB-BtS strains. The differential expression of GSTs in Cry1F-and Cry1Ie-resistant strains was shown in Table S9. From the list, five GSTs were up-regulated in ACB-FR and ACB-IeR, only ncbi_114360773 significantly up-regulated with log2FC value of 2.75 in ACB-FR strain and MSTRG.22972 up-regulated with log2FC value of 1.91 in the ACB-IeR strain. Three GSTs were downregulated in ACB-FR and ACB-IeR strains, but only ncbi_114365341 was significantly down-regulated in ACB-IeR strains with log2FC value -1.66.

Gene Ontology and Pathway Enrichment
In order to enrich the identified differential genes into specific functional categories and explore the biological significance behind transcriptomics data, ACB-FR or ACB-IeR was compared with ACB-BtS unigene for GO analysis, according to the Biological Process, Cellular Component and Molecular Function for classification and annotation. The three classification categories were subdivided into forty-five secondary GO-terms ( Figure 2).
In terms of Biological Process, unigenes were allocated to 22 aspects such as cellular process, single-organism process and metabolic process. Concerning cellular components, the cell part and cell and organelle represent most of the genes. In the Molecular Function category, catalytic activity and binding were the most abundant. ACB-FR and ACB-IeR, compared with ACB-BtS, have the largest number of genes annotated in the secondary classification of catalytic activity, single-organism process and cellular process.

Gene Ontology and Pathway Enrichment
In order to enrich the identified differential genes into specific functional categories and explore the biological significance behind transcriptomics data, ACB-FR or ACB-IeR was compared with ACB-BtS unigene for GO analysis, according to the Biological Process, Cellular Component and Molecular Function for classification and annotation. The three classification categories were subdivided into forty-five secondary GO-terms ( Figure 2).
In terms of Biological Process, unigenes were allocated to 22 aspects such as cellular process, single-organism process and metabolic process. Concerning cellular components, the cell part and cell and organelle represent most of the genes. In the Molecular Function category, catalytic activity and binding were the most abundant. ACB-FR and ACB-IeR, compared with ACB-BtS, have the largest number of genes annotated in the secondary classification of catalytic activity, single-organism process and cellular process. Comparing ACB-FR with ACB-BtS (Figure 2A Table S3).
A total of 130 candidate genes with pathway annotation were allocated to reference pathways in KEGG when ACB-FR was compared with ACB-BtS. In the result, 15 pathways with p ≤ 0.05 (Supplementary Materials: Table S4) were enriched, including "Longevity regulating pathway-multiple species" and "Protein processing in endoplasmic reticulum", which had the lowest p-value. The top 20 most-enriched pathways were shown in Figure 3A. However, comparing ACB-IeR with ACB-BtS, 94 candidate genes with pathway annotations were assigned to the reference pathways in KEGG, and 18 pathways were enriched with p ≤ 0.05 (Supplementary Materials: Table S4), including those with the lowest p-value "Pancreatic secretion" and "Protein digestion and absorption". The top 20 most-enriched pathways were shown in Figure 3B. A total of 130 candidate genes with pathway annotation were allocated to refe pathways in KEGG when ACB-FR was compared with ACB-BtS. In the resu pathways with p ≤ 0.05 (Supplementary Materials: Table S4) were enriched, incl "Longevity regulating pathway-multiple species" and "Protein processin endoplasmic reticulum", which had the lowest p-value. The top 20 most-enr pathways were shown in Figure 3A. However, comparing ACB-IeR with ACB-B candidate genes with pathway annotations were assigned to the reference pathwa KEGG, and 18 pathways were enriched with p ≤ 0.05 (Supplementary Materials: Tab including those with the lowest p-value "Pancreatic secretion" and "Protein digestio absorption". The top 20 most-enriched pathways were shown in Figure 3B.

RT-qPCR Verification of Differentially Expressed Genes
According to the fold change of RT-qPCR analysis, the results support differential expression at the gene level. Compared with susceptible strains, some differentially expressed genes have the same trend in resistant strains, but several of these Bt-resistance candidate genes were differentially regulated in both resistant strains. Additionally, the expression levels of glutathione S-transferase 1-like, UDP-glucuronosyltransferase 2B14like and ABCG1 were all up-regulated in both resistant strains (Figure 4). However, membrane-bound alkaline phosphatase-like (MBALP), heat shock protein 68-like and ABCA1 were down-regulated in both resistant strains (Figure 4). Moreover, pancreatic triacylglycerol lipase-like (PNLIP) was up-regulated in ACB-IeR and down-regulated in ACB-FR. The vertical axis is the pathway name, and the horizontal axis is the rich factor (the ratio of differentially expressed genes to all annotated genes in the pathway). The size of q-value is represented by the color of the point. Additionally, the smaller the q-value, the deeper of the red color. The size of the bubble indicates the number of differential genes contained in each pathway.

RT-qPCR Verification of Differentially Expressed Genes
According to the fold change of RT-qPCR analysis, the results support differential expression at the gene level. Compared with susceptible strains, some differentially expressed genes have the same trend in resistant strains, but several of these Bt-resistance candidate genes were differentially regulated in both resistant strains. Additionally, the expression levels of glutathione S-transferase 1-like, UDP-glucuronosyltransferase 2B14-like and ABCG1 were all up-regulated in both resistant strains (Figure 4). However, membranebound alkaline phosphatase-like (MBALP), heat shock protein 68-like and ABCA1 were down-regulated in both resistant strains (Figure 4). Moreover, pancreatic triacylglycerol lipase-like (PNLIP) was up-regulated in ACB-IeR and down-regulated in ACB-FR. Agriculture 2021, 11, x FOR PEER REVIEW 9 of 15 The significance threshold is set to p ≤ 0.05, "ns" indicates that the difference is not significant, "*" indicates the difference is significant, and the number of "*" indicates the degree of significance.

Discussion
The evolution of resistance of target pest populations to Bt toxins is the most relevant threat to the sustainability of Bt technology in the future [7], and clarifying insect resistance mechanisms will help improve Bt technology [31]. In the present study, the ACB-IeR strain has a higher estimated Cry1Ie resistance ratio, more than 862-fold increase compared to ACB-BtS, after 28 generations of selection (Table 1). In another study, the ACB-IeR strain also had more than 850-fold resistance to Cry1Ie, and Cry1F-resistant strains have not shown cross-resistance to Cry1Ie [32]. When the target pest is resistant to one Bt toxin, it will produce cross-resistance to other Bt toxins with the same or similar action sites and binding receptors, while it has no or low cross-resistance to Bt toxins with different types of action sites [33,34]. As previously reported, the high cross-resistance of European corn borer Cry1Ab-resistant strains to Cry1Ac is due to a common binding site [35,36]. In addition, Cry1F-resistant European corn borer does not have cross-resistance with Cry1Ab and Cry9C, because Cry1F has no common binding site with Cry1Ab and Figure 4. RT-qPCR results of differentially expressed genes. Three biological replicates of each strain of ACB-BtS, ACB-FR and ACB-IeR were performed in three technical replicates. The height of the bar graph represents the average of 2 −∆∆CT values. The significance threshold is set to p ≤ 0.05, "ns" indicates that the difference is not significant, "*" indicates the difference is significant, and the number of "*" indicates the degree of significance.

Discussion
The evolution of resistance of target pest populations to Bt toxins is the most relevant threat to the sustainability of Bt technology in the future [7], and clarifying insect resistance mechanisms will help improve Bt technology [31]. In the present study, the ACB-IeR strain has a higher estimated Cry1Ie resistance ratio, more than 862-fold increase compared to ACB-BtS, after 28 generations of selection (Table 1). In another study, the ACB-IeR strain also had more than 850-fold resistance to Cry1Ie, and Cry1F-resistant strains have not shown cross-resistance to Cry1Ie [32]. When the target pest is resistant to one Bt toxin, it will produce cross-resistance to other Bt toxins with the same or similar action sites and binding receptors, while it has no or low cross-resistance to Bt toxins with different types of action sites [33,34]. As previously reported, the high cross-resistance of European corn borer Cry1Ab-resistant strains to Cry1Ac is due to a common binding site [35,36]. In addition, Cry1F-resistant European corn borer does not have cross-resistance with Cry1Ab and Cry9C, because Cry1F has no common binding site with Cry1Ab and Cry9C [34,35]. To understand the mechanism of Bt resistance, it is necessary to study the binding sites of Bt toxins in O. furnacalis. This has guiding significance for the development and utilization of bivalent transgenic insect-resistant maize and different transgenic insectresistant maize rotations.
Furthermore, a comparative analysis of transcriptome data shows that more genes in ACB-FR and ACB-IeR strains are down-regulated compared to susceptible strains (Figure 1). In the current study, the transcriptome analysis revealed 62-63% of genes were downregulated in the Bt-resistant strains compared to the susceptible strain. Similarly, an analysis of chlorantraniliprole-susceptible and -resistant strains of Plutella xylostella revealed similar global gene down-regulation, where including many receptors and transporters [37]. Conversely, up-regulated genes related to immunity and detoxification genes were found in the current study and in Cry1A-resistant P. xylostella [38]. Thus, the down-regulation of Bt targets and transporters and simultaneous induction of detoxification mechanisms could lead to increased resistance to insecticides. In our transcriptome analysis, most genes within the ABCG family were down-regulated in ACB-IeR and ACB-FR compared with ACB-BtS; a similar result was reported in the ACB-AcR strain [20], which indicates that ABCG genes are putative receptors of Cry1Ac, Cry1Ie and Cry1F, and need more data to support this conclusion. Similarly, trypsin-like genes were up-regulated in ACB-IeR and ACB-FR compared with ACB-BtS (Table S2), but the functional relationship with Bt resistance is still unclear.
As the well-known receptor of Bt toxins, cadherin, APN and ALP was listed in Table S2. In total, 45 cadherins, 9 ALPs and 30 APNs were differentially expressed in ACB-IeR or -FR strains. Cadherins acted as receptors of Cry proteins are reported in many insect orders [39,40]. Binding with Cry1 toxins can form critical toxic oligomer that can effectively kill the insect [41]. Much research indicates that mutation of cadherin can generally link with Cry1A proteins resistance in Pectinophora gossypiella, Heliothis virescens and Helicoverpa armigera [14,42,43]. Knocking down by RNAi and knocking out by CRISPR/Cas 9 of cadherin can decrease the insect's susceptibility to Bt proteins [44]. Synergism between cadherin and ABCC2 were reported by in vitro and in vivo test [45,46].
In particular, Bt-resistance mechanisms have involved several ABC transporters from subgroups A [47], C [14,48] and G [49]; as reported previously [20], members of the three subgroups are down-regulated in Cry1Ab-and Cry1Ac-resistant O. furnacalis. The results of RT-qPCR showed that the abundance of Ofabca2 was reduced in the midgut of Cry1F-and Cry1Ie-resistant larvae (Figure 4), and interestingly, Ofabcg1 increased both in Cry1F-and Cry1Ie-resistant larval midgut transcripts and RT-qPCR results ( Figure 4). This may indicate that although Cry toxins have some of the same receptors, the molecular mechanisms involved in developing resistance characteristics of different resistant O. furnacalis strains may be different.
Less differential expression ALPs (nine ALPs) were found in the two resistant strains. GPI-tethered or soluble ALPs involve dephosphorylation. Decreasing expression of ALPs were found in resistant strains of Ostrinia furnacalis, Plutella xylostella, Spodoptera frugiperda and Helicoverpa armigera [50,51], but in a resistant strain of Diatraea saccharalis, ALP transcript levels were not reduced [52]. In our result, only ncbi_114355045 was significantly upregulated in ACB-FR and -IeR strains, ncbi_114362710 was significantly down-regulated in ACB-IeR strains and ncbi_114355386 was significantly down-regulated in ACB-FR strains. Further experiments could unearth the functions of these genes.
In addition, more differential expression APNs (30 APNs) were found in ACB-FR and -IeR strain. Not all the APNs could interaction with Cry toxins; some of them were reported associated with resistance [53]. In H. armigera Cry1Ac resistant strain, there was a 22 aa deletion of APN, and in Spodoptera exigua, APN1 was absent in Cry1Ca resistant strain midgut [54]. Further SNP or deletion can be analyzed between resistant and susceptible strains of O. furnacalis. The knocking down of APNs could decrease the susceptibility to Cry1A toxins in Chilo suppressalis [55], the gene down-regulated in ACB-FR or ACB-IeR, probably the receptors of Cry1F or Cry1Ie.
Moreover, previous studies have shown that most of the genes annotated as GSTs and trypsin in ACB-AhR have low expression levels [21]. Interestingly, we screened data related to GSTs from the transcript; most of the GSTs in ACB-FR showed an up-regulation trend, but it was down-regulated in ACB-IeR (Supplementary Materials: Table S9). The different expression GSTs in ACB-FR may be involved in Bt resistance, and could be partially responsible for Cry1F-sepecific resistance. The difference was that the genes annotated as trypsin alkaline C-like are mostly up-regulated in Cry1F-and Cry1Ie-resistant strains (Supplementary Materials: Table S8). According to previous reports, GSTs were related to resistance and participated in the detoxification of exogenous substances [56]. Trypsin is also considered the main protease for detoxification [57]. Thus, further research is needed to clarify the role of genes in the mechanism of Bt resistance.
In addition, GO and KEGG pathway analysis was performed on the DEGs of ACB-FR, ACB-IeR and ACB-BtS to find more genes related to Bt resistance in ACB. It can be seen that a large part of genes is related to catalytic activity in terms of molecular function (Supplementary Materials: Table S3). In the previous report, the molecular function part of the resistant strain of P. xylostella involved the most genes in the category of catalytic activity. These findings indicate that the resistance mechanism may also be related to catalytic activity [38,58]. Furthermore, KEGG pathway enrichment analysis showed that ACB-FR and ACB-IeR was significantly enriched in a series of genes in the same KEGG pathways (Supplementary Materials: Table S4), such as retinol metabolism, ABC transporters, metabolism of xenobiotics by cytochrome P450, protein digestion and absorption, pancreatic secretion, chemical carcinogenesis, fat digestion and absorption and drug metabolism-cytochrome P450; these pathways may be related to Bt toxins resistance. ABC transporters play a key role in the mode of action of the Bt toxin. Previous reports have shown that mutations in the ABC protein can confer high levels of resistance [59]. The two KEGG pathways about P450 are very important in the metabolism of xenobiotics and drugs. CYP450 includes many subfamilies, but only six exist in insects, five are insect specific, CYP6/CYP9/CYP12/CYP18 and CYP28, and one, CYP4, includes sequences from vertebrates [60,61], especially for CYP4 and CYP6 which are reported to relate with insecticide resistance [62,63]. Bt toxins also act as xenobiotics or drugs for insects, and enrichment on the two pathways in the comparison of resistant and susceptible strains indicates P450s are related to Bt resistance.

Conclusions
Cry1Ie-and Cry1F-resistant strains were selected in our lab and compared with susceptible strains. ACB-IeR strains evolved more than 862.07-fold greater resistance to Cry1Ie than susceptible strain ACB-BtS, and ACB-FR strains evolved more than 961-fold greater resistance to Cry1F than ACB-BtS. Further analysis of the transcriptomes of the three strains showed that a larger proportion of genes were significantly down-regulated in the Bt-resistant strain. This study is based on the transcripts of Cry1Ie-and Cry1F-resistant strains, focusing on finding the differences between the two strains at the molecular level, which is believed to provide ideas for delaying Bt resistance.