A Novel Reference for Bt-Resistance Mechanism in Plutella xylostella Based on Analysis of the Midgut Transcriptomes

Simple Summary Plutella xylostella is a very serious pest to cruciferous vegetables. At present, the control methods used are mainly traditional insecticides and the cultivation of Bt crops. However, with the long-term and large-scale use of insecticides, the diamondback moth has developed strong resistance to many kinds of insecticides and Bt crops. The Cry1S1000 strain of P. xylostella used here is a strain with more than 8000 times resistance to Bt Cry1Ac protoxin. In this paper, we used transcriptome sequencing to determine the midgut transcriptome of the G88-susceptible strain, Cry1S1000-resistant strain and its corresponding toxin-induced strains to find more genes related to Bt resistance. Our results will provide a reference for optimizing the control strategy of diamondback moth resistance and improving the control efficiency of biopesticides and Bt crops. Abstract The diamondback moth, Plutella xylostella, is a lepidopteran insect that mainly harms cruciferous vegetables, with strong resistance to a variety of agrochemicals, including Bacillus thuringiensis (Bt) toxins. This study intended to screen genes associated with Bt resistance in P. xylostella by comparing the midgut transcriptome of Cry1Ac-susceptible and -resistant strains together with two toxin-treated strains 24 h before sampling. A total of 12 samples were analyzed by BGISEQ-500, and each sample obtained an average of 6.35 Gb data. Additionally, 3284 differentially expressed genes (DEGs) were identified in susceptible and resistant strains. Among them, five DEGs for cadherin, 14 for aminopeptidase, zero for alkaline phosphatase, 14 for ATP binding cassette transport, and five heat shock proteins were potentially involved in resistance to Cry1Ac in P. xylostella. Furthermore, DEGs associated with “binding”, “catalytic activity”, “cellular process”, “metabolic process”, and “cellular anatomical entity” were more likely to be responsible for resistance to Bt toxin. Thus, together with other omics data, our results will offer prospective genes for the development of Bt resistance, thereby providing a brand new reference for revealing the resistance mechanism to Bt of P. xylostella.


Introduction
The diamondback moth, Plutella xylostella, is an important pest of cruciferous crops worldwide [1]. The total economic cost of its damage and management worldwide is more than US $4-5 billion per year [1], primarily due to strong resistance to multifarious Many studies have been conducted to examine transcriptome of P. xylostella, such as tissue-specific transcriptome analysis of a Bt-resistant strain [49]. Comparative studies were also reported on transcriptomic differences between Bt-susceptible and -resistant strains, although the focus is on different Bt toxins [43,[45][46][47]50,51]. Moreover, transcriptome analyses of different insecticides and resistant strains have been reported [52][53][54][55]. All this research provides valuable candidate gene resources for studying Bt resistance in insects. For the G88 and Cry1S10000 strains used in this research, similar studies have been carried out in the past. Lei et al. [43] used transcriptome analysis to compare the difference in midgut transcriptome between two resistant strains and one susceptible strain. One resistant strain used in this study had the same source as the resistant strain used in that study. However, the difference was that the resistant strain used in their study had been screened for toxin Cry1Ac in each generation, while the resistant strain used in this study remained highly resistant to Bt Cry1Ac toxin after only one screening with Cry1Ac toxin. Later studies [33,56] on the same genes (PxABCC2 and PxABCC3) also produced different results, so we believe that the resistance strains used in these two studies should have great differences and that our study results provide rich candidate gene resources for studying Bt resistance.
In this study, the transcriptomes of G88-susceptible and Cry1S1000-resistant strains were compared by RNA-Seq. The protoxin-treated strains were also compared to eliminate the expression differences of gene transcription levels caused by toxin induction. In addition, DEGs were further identified and validated by quantitative real-time PCR (qRT-PCR) by comparing both strains [33]. The results of this study provide candidate genes for the evolution of resistance to Bt and a reference for revealing the resistance mechanism to Bt in P. xylostella.

Insect Rearing and Sample Collection
The two P. xylostella strains (G88 and Cry1Ac-R) were provided by Dr. Anthony M. Shelton in 2016 [57]. Among them, the Cry1Ac-R strain was screened by 1000 µg/mL Cry1Ac protoxin and was never exposed to insecticide again. Here, we used Cry1S1000 to name the Bt resistant strain, as previously reported [33]. Larvae were fed an artificial diet at 26 ± 1 • C, 60 ± 5% RH (relative humidity), and 16:8 h (light:dark) photoperiod. During the adult mating stage, 10% honey water was used for supplemental nutrition. In this study, the resistance ratio of Cry1S1000 to Cry1Ac was >8000-fold compared to the G88 strain. The fourth-instar larvae of both strains were dissected to obtain the midgut tissues. In addition, 24 h before sampling, both strains were fed an artificial diet containing a lower concentration (0.006 µg/mL, LC 30 of G88 strain) of Cry1Ac protoxin to eliminate toxin-induced changes in transcription levels and ensure that most larvae were alive before sampling. These strains were also dissected to obtain midgut tissues.

Cry1Ac Protoxin Purification
Cry1Ac protoxin was produced by Btk strain HD73. The Bt culture was cultured in 1/2 LB (1 L:2.5 g yeast extract; 5 g tryptone; 5 g NaCl) at 230 rpm for 84 h [58,59]. Bacterial cells were suspended in 1 M NaCl, then centrifuged and washed twice with distilled water to prepare protoxin crystals. Then, the protoxin crystals were solubilized in the lysate (50 mM Na 2 CO 3 ; 50 mM EDTA; PH 9.5) with 5% β-mercaptoethanol. After centrifugation, the dissolved Cry1Ac protoxin was collected from the supernatant, and then 1/7 volume of 4 M sodium acetate (PH 4.5) was added to the precipitate. After two washes with distilled water, the precipitated protoxin was centrifuged and resuspended in 50 mM Na 2 CO 3 (PH 9.5). The Bradford method was used to determine the protoxin concentration by using BSA as a standard, and the quality of soluble protoxin was analyzed by 12% SDS-PAGE [60].

RNA Extraction, Library Construction, and Illumina Sequencing
Total RNA was extracted from the susceptible (G88) and resistant (Cry1S1000) strains, as well as from two protoxin-treated groups, before sampling using Trizol reagent (Invitrogen, Carlsbad, CA, USA), according to the manual's instructions. RNA samples were only used if they had a 260/280 ratio from 1.9 to 2.1 and an RNA integrity number (RIN) higher than 8. The quantity and quality of the total RNA were assessed by a Nano Drop2000 (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent 2100 bioanalyzer (Agilent, Palo Alto, CA, USA). Then, oligo (dT)-attached magnetic beads were used to purify Poly (A) mRNA. The purified mRNAs were segmented with fragment buffer at an appropriate temperature. Then, the first-strand cDNA was generated by reverse transcription using random hexamer-primes, and the second-strand cDNA was synthesized. Afterward, the short fragments were then connected to the sequencing adapters. After agarose gel electrophoresis, appropriate fragments were selected as templates for PCR amplification. Finally, Illumina HiSeq TM 2000 (BGI, Shenzhen, China) was used to sequence the library.

Gene Function Annotation and Characterization
The BLASTx search was carried out in protein databases, including KOG (eukaryotic Orthologous Group database), Nr (non-redundant) protein database, SwissPort, and KEGG (Kyoto Encyclopedia of Genes and Genomes protein database), to determine the functional annotation. Furthermore, a BLASTn search was performed using the Nt database. The InterPro annotation and GO (Gene Ontology) of unigenes were acquired using the Blast2GO and InterProScan5 program with Nr annotation, respectively [63]. Then, GO classification was performed to elucidate the distribution of DEG functions, including biological process, cellular component, and molecular function [64]. According to the annotation results of GO, KEGG, and official classification, functional classification of the differential genes was performed, and Phyper R in R (v.3.6.1) software was used for enrichment analysis. The p-value calculation method is as follows: Then, the false discovery rate (FDR) correction was performed based on the p-value, and the function with a Q-value of ≤0.05 was generally considered to be significantly The TPM method can eliminate the influence of different gene lengths and sequencing levels on gene expression calculation, and it can make the total expression level in different samples consistent. The TPM was calculated to represent the transcript-level expression [65]. The DESeq2 (v1.4.5) method is based on the negative binomial distribution principle, and differential expression analysis was performed using DESeq2 (v1.4.5) with Q value (Adjusted p value) ≤ 0.05, while the other parameters were the default values [66].

Real-Time Quantitative PCR Analysis of Gene Expression
qPCR was used to verify transcriptome sequencing results. Firstly, total RNA was extracted as described above from the fourth-instar larvae midguts. The genomic DNA was then removed using DNase I, and the concentration of total RNA was measured by NanoDrop 2000. The first-strand cDNA was synthesized from 2 ug RNA using "FastKing gDNA Dispelling RT SuperMix" (TianGen, Beijing, China) kit, following the manufacturer's instructions. In addition, primer pairs (Additional file 3) were designed using Primer 5. PCR was conducted using "Eastep qPCR Master Mix" (Promega, Shanghai, China) on a CFX96 (BioRad, Hercules, CA, USA) sequence detection system. The qPCR mixture included 10 µL qPCR master Mix, 0.4 µL of each primer (10 µM), 7.2 µL of RNase-free Water, and 2 µL cDNA. The qPCR program started at 95 • C for 2 min, followed by a total of 40 cycles of 95 • C for 15 s and 60 • C for 50 s. In the melting curve analysis, an automatic dissociation step cycle was added to determine the specificity of primers. Relative quantification of genes was performed using the 2 −∆∆Ct method [67], and the expression level was normalized to the ribosomal L32 gene (GenBank: AB180441) [68]. For each treatment, we used three biological replicates and three technical replicates. Moreover, one-way ANOVA with Holm-Sidak's test (p ≤ 0.05) was used to determine whether differences between treatments were statistically significant.

Illumina Sequencing Analysis
Illumina high-throughput sequencing was performed on cDNA samples from midgut tissues of G88 and Cry1S1000 strains of P. xylostella. We tested a total of 12 samples using BGISEQ-500, and every sample produced an average of 6.35 Gb data. Upon mapping these reads to the P. xylostella DBM V2, we obtained a mapping rate from 51.94% to 53.66% of the compared genome and 61.56% to 65.02% of the compared genes. We defined low-quality reads as those with a mass value of less than 15 bases that accounted for more than 20% of the total number of bases in the reads (Table 1). Moreover, the predicted new genes were 3641, and a total of 18,042 expressed genes were detected including 14,474 known genes and 3568 predicted new genes. In addition, 19,415 new transcripts were detected, of which 10,525 belonged to novel isoforms, 3641 were new protein-coding gene transcripts, and the remaining 5249 belonged to long noncoding RNAs ( Table 2 and Table S1).

Samples
Raw   represents G88-susceptible strain; DBMB represents G88-susceptible strain with toxin treatment; DBMC represents Cry1S1000-resistant strain; DBMD represents Cry1S1000-resistant strain with toxin treatment. GO analysis showed that the distribution patterns of DEG function categories in each group were basically consistent, except for "biological adhesion", "detoxification", "locomotion", "pigmentation", and "rhythmic process" (red arrows), between susceptible and resistant strains.

Functional Classification by KEGG
A KEGG pathway analysis was performed to obtain further insights into the global network. A total of 33 pathways were obtained between susceptible (G88) and resistant (Cry1S1000) strains, and they were distributed to five major groups: cellular components (789 DEGs), genetic information processing (543 DEGs), environmental information pro- Figure 1. Histogram presentation of gene ontology (GO) classification between any two groups from susceptible, resistant, and two toxin-treated strains. The GO categories shown on the X-axis are divided into three main ontologies: molecular function, biological process, and cellular component. The Y-axis indicates the number of DEGs in each category. DBMA represents G88-susceptible strain; DBMB represents G88-susceptible strain with toxin treatment; DBMC represents Cry1S1000-resistant strain; DBMD represents Cry1S1000-resistant strain with toxin treatment. GO analysis showed that the distribution patterns of DEG function categories in each group were basically consistent, except for "biological adhesion", "detoxification", "locomotion", "pigmentation", and "rhythmic process" (red arrows), between susceptible and resistant strains.

Functional Classification by KEGG
A KEGG pathway analysis was performed to obtain further insights into the global network. A total of 33 pathways were obtained between susceptible (G88) and resistant (Cry1S1000) strains, and they were distributed to five major groups: cellular components (789 DEGs), genetic information processing (543 DEGs), environmental information processing (721 DEGs), metabolism (1731 DEGs), and organismal systems (1479 DEGs). Among them, the majority of these DEGs were assigned to "global and overview maps" (619 DEGs), "signal transduction" (466 DEGs), and "transport and catabolism" (359 DEGs) ( Figure 2). These KEGG annotations provide a meaningful resource for studying the functions, specific processes, and pathways in the resistance research of P. xylostella.

Differentially Expressed Genes (DEGs) in Midgut Transcripts
Before the analysis of DEGs among four groups, the PCA analysis was performed on the midgut samples collected by transcriptome sequencing using princomp in R software (v.3.6.1), which showed that the three samples from each group belonged to a separate cluster. This result suggests that the data were sufficient and reproducible ( Figure S1).
Moreover, Pearson correlation analysis among samples showed that the R 2 values of biological replicates in each group were all greater than 0.98, which was much higher than 0.92 under ideal sampling and experimental conditions. Results showed that the similarity between the biological duplicates was very high, which met the ideal sampling require-

Differentially Expressed Genes (DEGs) in Midgut Transcripts
Before the analysis of DEGs among four groups, the PCA analysis was performed on the midgut samples collected by transcriptome sequencing using princomp in R software (v.3.6.1), which showed that the three samples from each group belonged to a separate cluster. This result suggests that the data were sufficient and reproducible ( Figure S1).
Moreover, Pearson correlation analysis among samples showed that the R 2 values of biological replicates in each group were all greater than 0.98, which was much higher than 0.92 under ideal sampling and experimental conditions. Results showed that the similarity between the biological duplicates was very high, which met the ideal sampling requirements for the subsequent analysis of DEGs ( Figure 3). DEGs were then analyzed among the four treatment groups by comparing the midgut transcriptome. The results revealed 599, 3284, 4313, and 63 significant DEGs between DBMA-DBMB, DBMA-DBMC, DBMB-DBMD, and DBMC-DBMD in the midgut transcriptome of P. xylostella (Padj ≤ 0.05) ( Figure S2).
Among these, 1528 genes were up-regulated, and 1756 were down-regulated between G88 and Cry1S1000 strains. Among all the DEGs, 4 were down-regulated, 2 were up-regulated by more than 10 times, 827 were down-regulated, and 631 were up-regulated by 2 to 10 times (Figure 4).  Figure S2).
Among these, 1528 genes were up-regulated, and 1756 were down-regulated between G88 and Cry1S1000 strains. Among all the DEGs, 4 were down-regulated, 2 were upregulated by more than 10 times, 827 were down-regulated, and 631 were up-regulated by 2 to 10 times (Figure 4).

Differentially Expressed Genes (DEGs) Involved in Insecticide Resistance
DEGs encoding Cry toxin receptors or participating in the Cry toxin pathway were further studied, including cadherin, aminopeptidase-N/P, alkaline phosphatase, ATPbinding cassette transporter family, heat shock protein, and V-type proton ATPase catalytic subunit A. Among them, only 4 DEGs for cadherin, 10 for aminopeptidase, 3 for alkaline phosphatase, 24 for ATP binding cassette transporter, and 7 for heat shock protein between G88 and Cry1S1000 strains were detected, all of which may be related to the resistance in P. xylostella (Table 3 and Table S2). gut transcriptome. The results revealed 599, 3284, 4313, and 63 significant DEGs between DBMA-DBMB, DBMA-DBMC, DBMB-DBMD, and DBMC-DBMD in the midgut transcriptome of P. xylostella (Padj ≤ 0.05) ( Figure S2).
Among these, 1528 genes were up-regulated, and 1756 were down-regulated between G88 and Cry1S1000 strains. Among all the DEGs, 4 were down-regulated, 2 were up-regulated by more than 10 times, 827 were down-regulated, and 631 were up-regulated by 2 to 10 times (Figure 4).   In this study, we examined not only DEGs that encode potential Cry toxin receptors but also DEGs that encode detoxification enzymes, such as carboxylesterase (CarEs), cytochrome P450 (P450s), glutathione S-transferase (GSTs), and some other genes related to metabolic insecticide resistance and insecticide targets. There were 24 DEGs encoding P450s, 15 of which were up-regulated, and 9 were down-regulated. In addition, 3 DEGs encoding CarEs were up-regulated, and the other 11 DEGs were down-regulated. The three DEGs encoding GST were all down-regulated (Table 3). Immune genes associated with hemolymph melanization, such as serpin protease and serpin protease inhibitor, were also detected. However, no serpin protease was found in the present study. Meanwhile, serpin protease inhibitor-related genes were searched, in which two of them were upregulated and the other two were down-regulated compared with the G88-susceptible strain (Table 3).

Validation of Differentially Expressed Genes by qRT-PCR
To eliminate the expression difference caused by toxin induction, we subtracted the DEGs between G88 and Cry1S1000 strains and the DEGs between these two strains and their toxin induction treatment group successively (DBMB, DBMD). Finally, we obtained a total of 3029 differentially expressed genes ( Figure 5).

Insects 2021, 12, x FOR PEER REVIEW 11 of 21
To eliminate the expression difference caused by toxin induction, we subtracted the DEGs between G88 and Cry1S1000 strains and the DEGs between these two strains and their toxin induction treatment group successively (DBMB, DBMD). Finally, we obtained a total of 3029 differentially expressed genes ( Figure 5). Firstly, we excluded the new genes predicted by BGI and the DEGs for which the TPM value was less than 10. Thus, after the exclusion step, 376 DEGs were left. Then, we carried out subsequent analysis combined with the log2 fold-change, in which we retained the DEGs having an absolute value of greater than 2 log2 fold-change. Thus, we obtained a total of 16 differential candidate genes, of which 6 genes were down-regulated, and the log2 fold-change was less than or equal to −2, while 10 genes were up-regulated and the log2 fold-change was greater than or equal to 2 (Table 4, Figure 6).  Firstly, we excluded the new genes predicted by BGI and the DEGs for which the TPM value was less than 10. Thus, after the exclusion step, 376 DEGs were left. Then, we carried out subsequent analysis combined with the log 2 fold-change, in which we retained the DEGs having an absolute value of greater than 2 log 2 fold-change. Thus, we obtained a total of 16 differential candidate genes, of which 6 genes were down-regulated, and the log 2 fold-change was less than or equal to −2, while 10 genes were up-regulated and the log 2 fold-change was greater than or equal to 2 (Table 4, Figure 6).  To verify the expression patterns of genes, we selected 16 DEGs (including 6 downregulated DEGs and 10 up-regulated DEGs) for qRT-PCR. In addition, the ribosomal L32 gene (GenBank: AB180441) was set as the candidate reference gene for qRT-PCR normalization. The results showed that the transcriptome expression pattern was consistent with  To verify the expression patterns of genes, we selected 16 DEGs (including 6 downregulated DEGs and 10 up-regulated DEGs) for qRT-PCR. In addition, the ribosomal L32 gene (GenBank: AB180441) was set as the candidate reference gene for qRT-PCR normalization. The results showed that the transcriptome expression pattern was consistent with the qRT-PCR results (Figure 7).

Discussion
P. xylostella is a highly destructive vegetable pest, and its resistance has badly been affecting the continuable use of insecticides. However, the resistance mechanism of P. xylostella to Cry1Ac remains to be improved. In this study, transcriptome analysis was per-

Discussion
P. xylostella is a highly destructive vegetable pest, and its resistance has badly been affecting the continuable use of insecticides. However, the resistance mechanism of P. xylostella to Cry1Ac remains to be improved. In this study, transcriptome analysis was performed on Cry1Ac-susceptible (DBMA) and -resistant (DBMC) strains and two toxin-treated (DBMB, DBMD) strains 24 h before sampling. A total of 12 samples were examined using the BGISEQ-500 platform, and every sample produced an average of 6.35 Gb of clean data. In addition, 3284 DEGs were detected between susceptible (DBMA) and resistant (DBMC) strains, 599 DEGs between susceptible (DBMA) and its toxin-treated strain (DBMB), and 63 DEGs between the resistant strain (DBMC) and its toxin-treated strain (DBMD). Among them, DEGs encoding Bt Cry toxin receptors cadherin, GPI-anchored alkaline phosphatase (ALP), glycosylphosphatidylinositol (GPI)-anchored aminopeptidase N (APN), and ATP-binding cassette transports (ABC) were investigated in detail. Furthermore, DEGs associated with metabolic insecticide resistance and insecticide targets and immune genes related to hemolymph melanization were also detected. Presumably, all the significantly differentially expressed genes detected above may be associated with the resistance to Bt toxin Cry1Ac.
Among all the DEGs, five were annotated as cadherin: four of them were up-regulated in the resistant strain, and only one was down-regulated. In addition, no DEGs were upregulated in the toxin-treated susceptible (DBMB) strain compared with the G88-susceptible strain (DBMA), which is consistent with the previous finding that the cadherin-like gene is not associated with resistance to Bt toxin Cry1Ac in P. xylostella [69]. Baxter et al. [70] concluded that the resistance of P. xylostella to Cry1Ac was not related to the mutation of cadherin but through other mechanisms. This finding contradicts the correlation between cadherin mutation and Cry1Ac resistance in Heliothis virescens and Pectinophora gossypiella. After analyzing the genetic mapping of resistance genes already reported in other species on P. xylostella, they also [71] found that these genes were not on the chromosome where the Cry1Ac resistance locus was located. In addition, the bioassay results of the knockout G88 susceptible strain of P. xylostella after the mutation of cadherin showed no significant difference in resistance ratio (unpublished). Unlike other lepidopteran insects, P. xylostella itself is complex and changeable, with many SNPs, which increases the difficulty of studying the resistance mechanism of P. xylostella to Bt by analogy. Therefore, it remains to be further investigated whether the orthologous genes of cadherin will play a role in resistance.
Aminopeptidase is a peptide chain end-hydrolase and is usually divided into four categories: N (AP-N; EC3.4.11.7), A (AP-A; EC3.4.11.2), P (AP-P; EC3.4.11.9), and W (AP-W; EC3.4.11.16) [72]. Aminopeptidase can catalyze the splitting of the amino-terminal residues of many proteins, and it is widely present in animals and plants. Among them, APN is the first reported receptor protein of Bt toxin. Since then, the APN of many lepidopteran insects has been identified as a binding protein of Bt toxin, which plays an essential role in insect resistance. Knight et al. [73] cloned the Cry1Ac toxin receptor protein gene in the midgut of the Manduca sexta and found that the protein was encoded by APN. Similarly, Gill et al. [74] used the same method and successfully found the receptor protein APN for the Bt toxin on the Heliothis virescen. In addition, the toxin binding test found that APN is the Cry1Ac toxin-binding receptor for P. xylostella [75], Lymantria dispa [76,77], Bombyx mori [78], and Trichoplusia ni [79]. Tiewsiri and Wang [14] found that Bt Cry1Ac resistance was related to the down-regulation of APN1 but not with the up-regulation of APN6 in cabbage loopers. Qiu et al. [80] demonstrated that RNA interference knockdown of APN genes decreases the susceptibility of Chilo suppressalis larvae to Cry1Ab/Cry1Ac and Cry1Ca-expressing transgenic rice. A recent study also found that knockdown of the APN genes decreases the susceptibility of C. suppressalis larvae to Cry1Ab/Cry1Ac and Cry1Ca [16]. In this study, 14 DEGs were annotated as aminopeptidases between susceptible and resistant strains, among which only two were up-regulated. Twelve DEGs encoding aminopeptidases were down-regulated compared with the G88 strain.
Another GPI-anchored Cry receptor involved in Bt resistance is alkaline phosphatase (ALP). In this study, we did not find any differentially expressed genes annotated as ALP. Another important group of receptors involved in Cry1Ac resistance is the ATP-binding cassette transporters. Many previous studies have proven that resistance to Cry1Ac is related to ABC transports. Among them, the structural variation or expression changes of ABCA, ABCB, ABCC, ABCG, and ABCH are all related to Bt resistance [29,31,33,34,40,[81][82][83][84][85][86]. The results from Yang [34] confirmed that ABCA2 is essential for the toxicity of Cry2Ab in T. ni, whereby mutation of ABCA2 confers resistance to Cry2Ab in the resistant T. ni strain derived from a Bt-resistant greenhouse population. In addition, Niu et al. [81] demonstrated that heterologous expression of DvABCB1 in Sf9 and HEK293 cells conferred sensitivity to the cytotoxic effects of Cry3A toxins in Western corn rootworm (WCR), Diabrotica virgifera verifier. Tian et al. [82] found that elevated expression of PxPgp1 was observed in P. xylostella after they were exposed to abamectin treatment. Zhou et al. [83] also found that reduced expression of the P-glycoprotein gene PxABCB1 is linked to resistance to Bt Cry1Ac toxin in P. xylostella (L.). Tanaka et al. [29] demonstrated Cry toxin receptor functionality for ABCC2 and highlighted the crucial role of this protein and cadherin. Meanwhile, bioassays of CRISPR-based mutant strains demonstrated that the deletion of PxABCC2 or PxABCC3 alone provided <4-fold tolerance to Cry1Ac, while disruption of both genes together conferred >8000-fold resistance to Cry1Ac, suggesting the redundant or complementary roles of PxABCC2 and PxABCC3 in P. xylostella [33]. Furthermore, RNA interference (RNAi)mediated suppression of Pxwhite gene expression significantly reduced larval susceptibility to Cry1Ac toxin. Genetic linkage analysis confirmed that down-regulation of the Pxwhite gene is tightly linked to Cry1Ac resistance in P. xylostella. Additionally, silencing PxABCH1 by a relatively high dose of dsRNA dramatically reduced its expression and resulted in larval and pupal lethal phenotypes in both susceptible and Cry1Ac-resistant P. xylostella strains [31]. In addition, Zhou et al. [83] found that P-glycoprotein gene PxABCB1 was significantly down-regulated in two resistant strains through preliminary transcriptome analysis. More importantly, knockdown of this gene in susceptible strain DBM1AC-S led to a significant reduction in sensitivity. Thus, they concluded that the decrease in PxABCB1 expression was closely related to Cry1Ac resistance. Xu et al. [87] constructed a stable PxABCC2 cell line and found that cell lines with stable PxABCC2 were significantly less sensitive to avermectin and chlorfenapyr compared to the control group. Their study showed that the upregulation of PxABCC2 gene was associated with insecticide resistance. This study also provides new insights into the cross-resistance between Bt toxins and chemical insecticides.
We also counted the number of different genes related to insecticide resistance in each group, as shown in Table 3. In the follow-up study, we conducted a functional verification of whether the CYP gene is involved in Bt resistance mechanism. The results provide new information on the interaction between Bt resistance and chemical insecticides. This study found 14 differentially expressed genes belonging to the ABC family between susceptible and resistant strains. Among them, five were up-regulated, and nine were down-regulated. Whether these genes are also involved in Bt resistance remains to be further studied. Our study provides rich candidate gene resources for the study of transporter involvement in Bt resistance.
The change of midgut protease-induced protoxin activation is one of the proposed mechanisms associated with Bt resistance, with some controversial evidence reported previously. For example, Gong et al. [88] found dramatically decreased activities of casein lysase and trypsin in the resistant P. xylostella strain SZ-R, due to significantly down-regulated expression of PxTryp_SPc1. A comparable finding, significantly lowered expression of PxTryp_CFT1 (Px007616) in the resistant P. xylostella strain CryS1000, was also identified in this study. However, Wei et al. [89] demonstrated that the LF120 strain (resistant strain) of H. armigera developed strong resistance to both activated toxin and protoxin, and the use of protease inhibitors did not change the LC 50 of the resistant strain to Cry1Ac pro-toxin. The role of midgut protease-induced protoxin activation, therefore, is well worth further investigation.
Our results also revealed DEGs encoding heat-shock proteins. A previous study showed that PxHsp90 assists Cry1A proteins by enhancing their binding to the receptor and protecting Cry protoxins from gut protease degradation [39]. In this study, five DEGs encoding heat-shock proteins were found between the susceptible and resistant strains. Among them, one was up-regulated, and four were down-regulated. Meanwhile, genes encoding serpin proteases and serpin protease inhibitors were involved in the immune function to regulate hemolymph melanization. However, we only found one DEG between G88 and Cry1S1000 strains.
Genes encoding the detoxification enzymes of carboxylesterase (CarE), cytochrome P450 (P450), and glutathione S-transferase (GST) were also identified. Among them, P450 family members play a significant role in the degradation of chemical insecticides. In this study, 24 DEGs were identified as P450s, 4 were CarEs, and 4 were GSTs. A previous study showed that RNA interference-mediated knockdown of CYP6BG1 from P. xylostella reduced larval resistance to permethrin [90]. In addition, Chen [91] demonstrated that the overexpression of three CYP genes, CYP6CY14, CYP6CY22, and CYP6UN1, contributed to dinotefuran resistance in Aphis gossypii. Furthermore, knockdown of CYP4PR1 increased the susceptibility to deltamethrin in pyrethroid-resistant insects [92]. A dual-Luciferase Reporter assay, a yeast one-hybrid (Y1H) assay, and RNA interference confirmed that the mutation of the PxABCG1 promoter in the resistant strain resulted in Bt being unable to combine with the toxin, resulting in strong resistance [93]. Qin et al. [94] showed that MAPK-activated PxJun inhibited the expression of PxABCB1, leading to the resistance of P. xylostella to Cry1Ac. Their study is the first to identify transcription factors involved in the transcriptional regulation of Cry receptor genes in the gut of Bt-resistant insects. Similar results were also found in several CYP genes found in our transcriptome data. Furthermore, preliminary sequence verification found that there were significant variations in the promoter regions of several CYP genes between susceptible and resistant strains. Notably, it is widely reported that the P450 gene is mainly related to insecticide resistance. However, whether or not the P450 gene is involved in the detoxification metabolism of Bt toxin, and what role it plays in this process, needs further validation.
To find other resistance-related genes, we analyzed the GO function and KEGG function classifications for DEGs. The analysis of the GO function classification showed that 1932 genes were assigned to molecular functions, among which "binding", "catalytic activity", and "transporter activity" had the largest number. Furthermore, 1461 genes were assigned to biological processes and 1138 genes to cellular components. Specifically, "cellular process" and "metabolic process" categories in the biological process domain were represented by 574 and 468 DEGs between G88 and Cry1S1000 strains, respectively. In the cellular component domain, the top three categories were "cellular anatomical entity", "intracellular", and "protein-containing complex", and the combined total number of DEGs was more than 1100. Therefore, future work should be focused on DEGs under these categories, specifically searching genes related to Bt resistance. In our study, a large number of genes were significantly enriched in the "hydrolase activity", "catalytic activity", "metabolic process", "primary metabolic process", and "organic substance metabolic process"-related GO categories ( Figure S3).
Meanwhile, the KEGG pathway analysis indicated that 466 DEGs were annotated as "signal transduction" under the environmental information processing. In the metabolism domain, 619 DEGs were annotated as "global and overview maps", representing the largest number of DEGs between G88 and Cry1S1000 strains in the KEGG function classification. This result suggests that there might be many genes related to Bt resistance in this category. In addition, DEGs annotated as "transport and catabolism", "digestives system", and "endocrine system" also need further investigation in future research. KEGG category analysis revealed a significant range of enrichment genes in "fanconi anemia pathway"related KEGG pathways ( Figure S4).
After continuous laboratory screening with activated Cry1Ac toxin for 20 generations, the SZ strain that originally had only 20 times resistance to Cry1Ac reached a 1200 times resistance level. Moreover, studies showed that resistance of the SZBT strain to Cry1Ac was controlled by a single, autosomal, and incomplete recessive gene. This study provides good material for studying the mechanism of Bt resistance at single loci [95]. At the same time, it also shows that the change of resistance may be relatively simple under laboratory selection methods. In particular, we can find Bt resistance genes with stronger specificity through this method to analyze the mechanism of Bt resistance. The G88 susceptible strain used in this paper has been raised for more than 110 generations in the laboratory, and the Cry1S1000 strain has been raised for more than 90 generations. Stable indoor breeding strains provide good experimental materials for subsequent studies, and they also provide more possibilities for the search of Bt resistance genes.
In the previous study [43], a comparative transcriptomic analysis was conducted between a susceptible strain (MM) and two resistant strains (GK and MK) of P. xylostella, of which the sources of MM and MK were consistent with ours. However, Lei et al. continuously exposed resistant strains to Bt toxin to maintain their resistance multiples. The resistant strain (Cry1S1000) used in this study [33] was screened with Cry1Ac toxin only once when they first arrived in our laboratory, and then it was bred and passed on in an environment without toxin exposure until now. Due to the selective induction of Bt toxin resistance genes, we can actually understand these two resistance genes as different resistant lines; that is, when we compare them with a susceptible strain at the same time, the DEGs should be different. In addition, when we verified the function of the same gene (PxABCC2 and PxABCC3), we obtained completely different results [33,56], so it is reasonable for us to compare two completely different resistance strains. As shown in Table 3, we obtained a large number of CYP genes. However, in subsequent functional verification experiments (unpublished), we found sequence differences in coding regions and promoter regions of multiple CYP genes between susceptible and resistant strains, which may be related to resistance differences between the two strains. Unlike the previous studies [43], before sampling we performed a short 24 h induction process for transcriptome sequencing using low concentrations of Bt toxin in both resistant and susceptible strains. This process will eliminate some differences in gene transcription levels due to toxin induction, thereby gaining more constitutive types of DEGs between both strains.

Conclusions
In conclusion, this study reveals many DEGs between G88-susceptible and Cry1S1000resistant strains. In view of our results, several factors may participate in P. xylostella resistance to Cry1Ac toxins. The widely reported cadherin and aminopeptidase N, alkaline phosphatase, and ATP-binding cassette transporter may be involved in Cry1Ac resistance. Among them, we have shown 14 DEGs encoding ABC transporters between G88 and the Cry1S1000 strains, suggesting that genes from the ABC family are more likely to play an important role in Bt resistance. Moreover, DEGs related to detoxification metabolism may account for P. xylostella resistance to the Cry1Ac toxins, such as P450s (24 DEGs), CarEs (4 DEGs), and GSTs (4 DEGs). These important genetic resources, coupled with ongoing antibody analysis resources, are more conducive to screening and obtaining genes related to Cry1Ac resistance in P. xylostella. RNAi or CRISPR/Cas9 technologies could be used to verify relationships between these candidate genes and Bt resistance.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/insects12121091/s1, Figure S1. PCA analysis of midgut samples of P. xylostella; Figure S2. The number of differentially expressed genes between groups; Figure S3. GO enrichment of differentially expressed genes between G88 susceptible (DBMA) and Cry1S1000 resistant (DBMC) strains; Figure S4. KEGG enrichment of differentially expressed genes between G88 susceptible (DBMA) and Cry1S1000 resistant (DBMC) strains; Table S1. Gene ID and annotation information of the novel transcript; Table S2. Gene ID of differentially expressed genes between groups; Table S3. Primer sequences for qRT-PCR of 19 selected differentially expressed genes.