Profiling of MicroRNAs in Midguts of Plutella xylostella Provides Novel Insights Into the Bacillus thuringiensis Resistance

The diamondback moth (DBM), Plutella xylostella, one of the most destructive lepidopteran pests worldwide, has developed field resistance to Bacillus thuringiensis (Bt) Cry toxins. Although miRNAs have been reported to be involved in insect resistance to multiple insecticides, our understanding of their roles in mediating Bt resistance is limited. In this study, we constructed small RNA libraries from midguts of the Cry1Ac-resistant (Cry1S1000) strain and the Cry1Ac-susceptible strain (G88) using a high-throughput sequencing analysis. A total of 437 (76 known and 361 novel miRNAs) were identified, among which 178 miRNAs were classified into 91 miRNA families. Transcripts per million analysis revealed 12 differentially expressed miRNAs between the Cry1S1000 and G88 strains. Specifically, nine miRNAs were down-regulated and three up-regulated in the Cry1S1000 strain compared to the G88 strain. Next, we predicted the potential target genes of these differentially expressed miRNAs and carried out GO and KEGG pathway analyses. We found that the cellular process, metabolism process, membrane and the catalytic activity were the most enriched GO terms and the Hippo, MAPK signaling pathway might be involved in Bt resistance of DBM. In addition, the expression patterns of these miRNAs and their target genes were determined by RT-qPCR, showing that partial miRNAs negatively while others positively correlate with their corresponding target genes. Subsequently, novel-miR-240, one of the differentially expressed miRNAs with inverse correlation with its target genes, was confirmed to interact with Px017590 and Px007885 using dual luciferase reporter assays. Our study highlights the characteristics of differentially expressed miRNAs in midguts of the Cry1S1000 and G88 strains, paving the way for further investigation of miRNA roles in mediating Bt resistance.


INTRODUCTION
Bacillus thuringiensis (Bt), a class of spore-forming gram-positive bacterium that can produce different insecticidal crystal proteins (Cry and Cyt toxins), has been widely used as an entomopathogen for pest control (Wu et al., 2008;Raymond et al., 2010;Palma et al., 2014). The application of Bt toxins cannot only increase crop yields and bring substantial economic benefits, but also reduce environmental pollution caused by the abuse of chemical insecticides (Bravo et al., 2011). Bt Cry toxin has specific control efficiency on lepidopteran pests, and the mode of action involves toxin solubilization, proteolytic activation, interaction with midgut proteins of insects, formation of a prepore oligomeric structure, facilitation of the insertion into cell membrane, and creation of an ionic pore that kills midgut cells (Bravo et al., 2011). There are several resistance strategies taken by insects to counter Bt toxins, including mutations in receptor [e.g., cadherin, aminopeptidase N (APN) or ATPbinding cassette (ABC) transporters], alteration in Cry toxin activation and binding ability to gut membrane, sequestration of toxins by glycolipid moieties or esterases, and the elevation of insects immune responses (Pardo-Lopez et al., 2013). So far, more than 22 cases of field-evolved Bt resistance in several insect species have been documented, potentially reducing the control efficacy of Bt toxins toward agricultural pests (Janmaat and Myers, 2003;Carrière, 2017, 2019;Calles-Torrez et al., 2019). Therefore, to delay the resistance evolution, it is urgent to understand the mechanisms causing Bt resistance and provide novel targets for pest control.
The diamondback moth (DBM), Plutella xylostella, a notorious pest of cruciferous crops, has caused severe economic losses globally. One of the main challenges of DBM control is its rapid evolution of resistance against a wide range of insecticides, including Bt-based products (Furlong et al., 2013). Considering its severe insecticide-resistance, short life cycle and host range expansion , microRNA (miRNA)-mediated efficient and environmentally friendly approaches have been proposed to combat DBM (Vaschetto and Beccacece, 2019).
MiRNAs are a class of endogenous non-coding RNAs ranged from 19 to 24 nt in length and play crucial roles in the post-transcriptional regulation (Bartel, 2004;Cullen, 2004;Asgari, 2013). MiRNAs regulate gene expression by partially or completely binding their seed sequence region with the 3 untranslated region (3 UTR) of corresponding target genes (Bartel, 2004). A large number of miRNAs have been identified from multiple insect species, participating in insect development (Liu et al., 2018(Liu et al., , 2020bXu et al., 2019;Lim et al., 2020), reproduction (Lucas et al., 2015;He et al., 2016;Ling et al., 2017), behavior (Cristino et al., 2014;Chen and Rosbash, 2016;Niu et al., 2019), and host-pathogen interaction (Singh et al., 2014;Dubey et al., 2019;Wu et al., 2019). Recently, a growing body of evidence has illustrated the pivotal role of miRNAs in responses to environmental chemical and pathogen exposures, which attracts great interests in the field of insect toxicology. Insecticide resistance is predominantly caused by the rising metabolic detoxification of insecticides and declining in the target sites sensitivity of the insecticides (Liu, 2015). MiRNAs enhance insecticide resistance by negatively regulating the expression of detoxification-related genes or altering insect susceptibility to various kinds of insecticides. In Culex pipiens pallens, the pyrethroid-resistance is modulated by miR-92a via suppressing CpCPR4 expression (Ma et al., 2017), miR-278-3p and miR-932 are also involved in the pyrethroid resistance by targeting CYP6AG11 and CpCPR5, respectively (Lei et al., 2014;Liu et al., 2016). In addition, it has been reported that miR-13664 could interact with CpCYP314A1 and subsequently regulate deltamethrin resistance . Another study confirms the significant role of miR-2∼13-71 cluster in deltamethrin resistance by regulating the expression of CYP9J35 and CYP325BG3 . In DBM, miR-7a and miR-8519 regulate chlorantraniliprole resistance via overexpressing ryanodine receptor (RyR) (Li X. et al., 2015). MiR-2b-3p and miR-14-5p are proposed to inhibit the detoxification pathways by suppressing the transcript levels of CYP9F2 and CYP307a1, respectively (Etebari et al., 2018). MiR-189942 directly targets the ecdysone receptor (EcR) isoform B and increased the tolerance to fufenozide (Li X. et al., 2020). Besides, miR-998-3p is implicated in Bt resistance through targeting ABCC2, a putative receptor for Bt Cry1Ac toxin in three lepidopterans P. xylostella, Helicoverpa armigera, and Spodoptera exigua . While these studies provide insight into the critical roles of miRNAs in regulating metabolic resistance to insecticides, most of them are focused on the chemical insecticide. The mechanism of miRNAmediated Bt resistance still remains limited.
A previous study has profiled the dynamic expression patterns of miRNAs in response to Bt at different time courses using the whole body of DBM, indicating the potential role of miR-2b-3p in regulating insect immunity . In the present study, we aimed to expand the numbers of identified miRNAs in specific tissues, i.e., midguts of the Bt resistant and susceptible DBM stains and explore the potential mechanisms of miRNAmediated regulation in Bt resistance. The differentially expressed miRNAs in midguts of both strains were putatively linked to Bt resistance. One of the differentially expressed miRNAs, novel-miR-240, which showed an inverse correlation with its target genes, was then selected for dual luciferase assays to validate the interaction with its target genes. Our work is important for further understanding the potential roles of miRNAs in the evolution of Bt resistance in DBM field populations, providing novel insights for their putative applications in integrated pest management.

Insects Rearing
The insecticide-susceptible DBM strain (Geneva 88, G88) was originally collected from the New York State Agricultural Experiment Station in 1988 and maintained on an artificial diet without exposure to insecticide. A colony of DBM, which was established from a crucifer field in Florida in 1992, was used to investigate its resistance to B. thuringiensis subsp. kurstaki spray formula. Specifically, it was challenged with the transgenic broccoli expressing Cry1Ac toxin, and the survived moths were maintained as a Cry1Ac-R strain (Zhao et al., 2002). The Cry1Ac-R resistance strain was further selected with 1,000 µg/ml Cry1Ac protoxin in 2016 and then reared on an artificial diet without additional insecticide exposure for more than 50 generations. This strain was thereafter named as Cry1S1000 (Liu et al., 2020a). Both the G88 susceptible and Cry1S1000 resistance strains were reared with the same approach reported before (Xu et al., 2020). In brief, artificial diet and 10% honey solution were, respectively, used to feed larvae and adults at 25 ± 1 • C, with a 16:8 h light: dark cycle and 65 ± 5% relative humidity.

Sample Collection and Small RNA Library Construction
Fresh midguts dissected from 60 4th instar larvae of the Cry1S1000 and G88 strains were collected in separate tubes and the midgut contents were gently scraped off. The midgut samples were immediately immersed in liquid nitrogen and subsequently frozen in -80 • C until further analysis. To generate a small RNA (sRNA) library, the miRNeasy Mini Kit (Qiagen, Germany) was used for total RNA extraction from the midgut samples. RNA concentration and purity were measured using the NanoDrop 2000 (Thermo Fisher Scientific, United States) and Agilent Bioanalyzer 2100 system (Agilent Technologies, United States), respectively. To generate sRNA libraries, 3 µg RNA per sample was used for adaptor ligation in both of 5 and 3 ends using the NEBNext Multiplex Small RNA Library Prep Kit for Illumina (NEB, United States). Then, the M-MuLV Reverse Transcriptase (RNase H − ) was used for the first strand cDNA synthesis. The PCR amplification was performed and the fragments of 140∼160 bp were recovered and sequenced using the Illumina HiSeq 2500 sequencing platform (Illumina Inc., United States).

Bioinformatic Analysis
To obtain the clean reads, low-quality reads, reads containing more than 10% unknown base N, without 3 adaptor sequence, shorter than 18 nt, or longer than 30 nt and 3 adapter reads were removed from the raw sequences, and the identical reads were collapsed in order to remove redundancy. After that, the trimmed sequences were aligned to Silva, GtRNAdb, Rfam, and Repbase databases to filter the ncRNAs including ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA), and repetitive RNAs. The remaining reads mapped to the reference genome 1 with Bowie software were used for further miRNA identification. Known miRNAs were identified by aligning with mature sequences in miRBase (release22) 2 with one mismatched allowed. For novel miRNAs identification, the miRDeep2 package was utilized to obtain potential precursor sequences. Prediction of novel miRNAs was completed based on the precusor stucture energy information and the distribution infromation of reads on precursor sequence (Friedlander et al., 2012).

Expression Patterns of miRNAs Based on Transcripts per Million
To identify the miRNAs differentially expressed in the Cry1S1000 and G88 strains, transcripts per million (TPM) was used to estimate miRNA expression levels. The normalization formula was: normalized expression = (number of mapped reads for each miRNA/total number of mapped reads) × 10 6 . The DESeq2 algorithm (Love et al., 2014) was performed to estimate the differences in reads frequencies comparing the Cry1S1000 to G88 strains, and p-value was adjusted using the Benjamini-Hochberg false discovery rate (FDR) procedure (Puoliväli et al., 2020). Only those miRNAs with a log2 fold-change > 1.58 or < -1.58 in expression and a FDR ≤ 0.01 were determined to be significantly differentially expressed between the Cry1S1000 and G88 strains.

Target Gene Prediction and Enrichment Analysis
To predict and analyze the potential target genes of differentially expressed miRNAs, miRanda (Rehmsmeier et al., 2004), and TargetScan (Lewis et al., 2005) algorithms with default parameters were employed. The DBM genome sequence (DBM-DB) 3 was used as a reference and the coding sequence (CDS) was searched for predicting miRNA targets. To screen more reliable results, only those genes predicted by both algorithms were retained. The predicted miRNA targets were annotated based on NR, Swiss-Prot (Apweiler et al., 2004), Gene Ontology (GO) (Ashburner et al., 2000), COG (Tatusov et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004), KOG (Koonin et al., 2004) and Pfam databases using BLAST program with a cut-off E-value = 10 −5 . To further investigate the functions of the putative target genes, the software Database for Annotation Visualization and Integrated Discovery (DAVID) 4 was used to conduct the GO enrichment consisted of biological processes, cellular components and molecular functions. Pathway analysis was based on the KEGG database that determine the pathways associated with the target genes of the differentially expressed miRNA. The DBM genome database was used as a background to determine the most enriched GO terms with a corrected p-value ≤ 0.05 as the threshold. The KEGG analysis was performed and a corrected p-value was set as 0.05 for identifying significantly enriched pathways within the predicted differentially expressed miRNA target dataset.

Expression Patterns of miRNAs and Target Genes Based on RT-qPCR
Total RNA was extracted from midguts of 20 fourth instar larvae using a PF miRNA Isolation Kit (Omega, United States), and the concentration was measured by NanoDrop 2000 (Thermo Fisher Scientific, United States). The miRNA reverse transcription reactions were conducted using the miScript II RT kit (Qiagen, Germany), and the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, United States) was used to synthesize the first-strand cDNA of target genes.
Reverse transcriptase-quantitative PCR (RT-qPCR) reactions of miRNAs were performed using the miScript SYBR Green PCR kit (Qiagen, Germany) with miRNA-specific forward primers (Supplementary Table 1), and the RT-qPCR reactions of target genes were performed using the Eastep qPCR Master Mix (Promega, United States. All amplification reactions were performed on the Bio-Rad Realt-time PCR system (Bio-Rad, United States) with three technical repeats and three independent biological replicates. The normalized expression levels were calculated with 2 − Ct method using the internal references U6 snRNA for miRNAs normalization and RPL32 for target genes normalization. The significance was analyzed using Student's t-test with SPSS 25.0 program. All the reagents and kits used in the current study were applied according to the manufacturer's instruction unless stated specifically.

Dual Luciferase Reporter (DLR) Assay
Novel-miR-240 agomir (5 -UCCUCAAUAUCAUAUUCC UCGC-3 ) and negative control (NC) agomir (5 -UUCUCCGAACGUGUCACGUTT-3 ) were synthesized by Sangon Biotech (Shanghai, China). An agomir is a sequnce of double-strand RNA with special chemical modification, and it can mimic the function of miRNA and up-regulate the expression of endogenous miRNA. The luciferase reporter plasmids were constructed by cloning the wild type (WT) and the mutated (MUT) target sequences of Px017590 and Px007885 into the pmirGLO vector (Promega, United States). The fragments of Px017590 and Px007885 containing WT and MUT sequences were then confirmed by sequencing. To obtain high transfection effciency and low background expression, HEK293T cell line was used for the DLR assay. HEK293T cells were cultured in a 24-well plate and co-transfected with the reporter plasmid and the miRNA agomir or NC agomir using the Attractene Tranfection Reagent (Qiagen, Germany). Each well containing 1 µg reporter plasmid and 50 nM miRNA agomir or NC agomir. The activeties of the Firefly and Renilla luciferases were preformed by using a Dual-Glo Luciferase Assay System (Promega, United States) at 48 h post-transfection. Firefly luciferase activity was normalize to Renilla luciferase activity.

Overview of the sRNA Libraries
To analyze the role of miRNAs in the Bt resistance of DBM, we performed the deep sequencing for identification and characterization of miRNAs in the larval midguts of the Cry1S1000 and G88 strains. Using the high-throughput sequencing, a total of 125,558,480 raw reads (≥17.6 million per library) were obtained. After removing low-quality reads and 3' adaptor, only the reads ranging from 18 to 30 nt were kept. Clean reads from three Cry1S1000 libraries (15,964,505, 18,860,033, and 16,118,802) and three G88 libraries (12,107,465, 13,815,968, and 18,445,775) were obtained, respectively. The clean reads were aligned with Silva, GtRNAdb, Rfam, and Repbase databases and then mapped to the reference genome using the Bowtie package. A total of 9,200,300 mapped reads were retained and used to predict mature miRNAs in subsequent analyses ( Table 1).

Identification of Known and Novel miRNAs
Bioinformatic analysis was carried out to identify miRNAs expressed in the larval midguts of the Cry1S1000 and G88 strains. By searching against miRBase, 76 miRNAs in our libraries were identical to mature miRNA sequences and annotated as known miRNAs. Additionally, the miRDeep2 software package (Wen et al., 2012) was used to predict novel miRNAs by exploring the potential precursor sequences and estimating their randfold value. In total, 361 novel miRNAs were identified ( Table 1). The length distribution of both known and novel miRNAs showed a peak at 22 nt, which is the standard size of animal miRNAs (Figures 1A,B). The first base preference of identified miRNAs was also investigated, and the result showed a dominant bias toward uracil (U). This bias was particularly evident in both known miRNAs (20 and 22 nt in length) and novel miRNAs (19 nt and 22 nt in length). Calculation of the miRNA nucleotide bias at each position demonstrated that U and A (adenine) were more common than C (cytosine) and G (guanine) in almost all positions, especially at base 14 ( Figures 1C,D).
Next, we classified 178 known and novel miRNAs into 91 families based on the sequence conservation across different species, while the remaining 259 miRNAs could not be classified into any families. Six main families, which hold relatively more members compared with other sets, were discovered. To be specific, 12 identified miRNAs belonged to the miR-4864 family and ten miRNAs belonged to the miRNA-8517 family. It was also estimated that the miR-279 and miR-9 families, respectively, contained six miRNAs, while the other two miRNA families (miR-2 and miR-185) were each represented by five miRNAs. In comparison, the remaining families were predicted to have less than five miRNA members each (Figure 2). All of the 178 identified miRNAs and their miRNA families were listed in Supplementary Table 2.
In addition, the miRNA expression levels were assessed by TPM formula. The most abundant known miRNA was pxy-miR-279b-3p, followed by pxy-miR-750, pxy-miR-8494-5p, pxy-miR-306, and pxy-miR-8532-3p, while the five most highly expressed novel miRNAs were novel-miR-1, novel-miR-246, novel-miR-176, novel-miR-28, and novel-miR-132. Interestingly, the average expression value of the listed novel miRNAs was significantly higher than the known miRNAs. The most abundant known and novel miRNAs in the libraries were represented in Figure 3 and the details of all the identified miRNAs were listed in Supplementary Table 3.

Differentially Expressed miRNAs in the Cry1S1000 and G88 Strains
The expression levels of known and novel miRNAs were compared after normalization by TPM formula described above.  Twelve differentially expressed miRNAs including nine downregulated and three up-regulated miRNAs were found in the Cry1S1000 strain, compared to the G88 strain. In these miRNAs, novel-miR-210 and novel-miR-48 were predicted to show the strongest up-regulation and down-regulation trends, respectively (with the highest Log2-fold change value among listed miRNAs) ( Table 2). Additionally, novel-miR-97 and novel-miR-237 were classified into the miR-973 and the miR-587 families, respectively, while others did not belong to any known miRNA family. Furthermore, the cluster analysis was performed to investigate the expression patterns of 12 differentially expressed miRNAs. Among all tested miRNAs, novel-miR-240, novel-miR-237, novel-miR-48, and novel-miR-97 showed significantly lower expression levels in the Cry1S1000 strain than in the G88 strain, while novel-miR-210 displayed a relatively higher expression level in the Cry1S1000 strain (Figure 4). It was noted that novel-miR-118 and novel-miR-24 shared the same sequence with novel-miR-97, novel-miR-157 shared the same sequence with novel-miR-210, and the sequence of novel-miR-266 was the same with novel-miR-225. Thus, these miRNAs sharing the same sequences were combined and, respectively, named as novel-miR-97, novel-miR-210, and novel-miR-225 hereafter.

Target Gene Prediction of the Differentially Expressed miRNAs
A total of 7,647 target genes of the miRNAs identified in our study were predicted by miRanda and TargetScan, of which 7,485 target genes were annotated with NR, Swiss-Prot, GO, COG, KEGG, KOG, and Pfam databases (Finn et al., 2014).
To further explore the roles of the differentially expressed miRNAs, 44 genes potentially targeted by these 12 differentially expressed miRNAs were selected and represented in Table 3. Among them, the number of predicted target genes for each miRNA ranged from zero to 27, which implying a complex regulatory network between miRNAs and their target genes. Most miRNAs had multiple target genes. Conversely, some of the genes were also modulated by multiple miRNAs. Additionally, the annotation using the DBM-DB database showed that these putative transcripts were likely involved in multiple biological processes. Some of the differentially expressed miRNAs were predicted to regulate genes encoding chitinase A1, cuticle protein 6, transporting P-type ATPase, alyl-tRNA synthetase and cell wall protein, which might be closely related to the defense of external toxicant, energy transportation and metabolism and likely be involved in the insect immunity and insecticide resistance.

GO Enrichment and KEGG Pathway Analysis of Target Genes of Differentially Expressed miRNAs
GO annotation enrichment was performed to evaluate the presumptive functions of target genes of differentially expressed miRNAs (Figure 5). GO term analysis for these target genes revealed seven cellular components, five molecular processes, and seven different biological processes, among  Potential miRNA, sequences of differentially expressed miRNAs were similar to those in other species but differ in some nucleotide positions. * Means miRNA sequences were not found in other species in miRbase (V22) with E-value cutoff < 3.
which cellular process (GO:0009987), metabolic process (GO:0008152), single-organism process (GO:0044699), and biological regulation (GO:0065007) were dominant. Additionally, most of the target genes were classified in membrane (GO:0016020), binding (GO:0005488), and catalytic activity (GO:0003824), implying the putative roles of these differentially expressed miRNAs in transmembrane transport and translation. A KEGG pathway analysis was also conducted to elucidate the biological interpretation of genes targeted by differentially expressed miRNAs (Figure 6). The result exhibited that Px008940 and Px016226 were enriched in the Hippo signaling pathway, Px007475 and Px008286 were enriched in aminoacyl-tRNA biosynthesis and Px006958 was enriched in MAPK signaling pathway, which were associated with energy metabolism in insects. In addition, the remaining target genes were enriched in the metabolic pathways such as ubiquitin mediated proteolysis, ribosome biogenesis in eukaryotes, protein processing in endoplasmic reticulum and fatty acid metabolism (Supplementary Table 4).
The potential target genes of differentially expressed miRNAs were investigated with two miRNA targets prediction algorithms, miRanda and TargetScan. A total of 34 target genes were profiled using RT-qPCR. Twelve target genes (Px013363, Px014239, Px008286, Px011169, Px002074, Px009200, Px011234, Px006256, Px017590, Px007885, Px015953, and Px010315) were significantly induced in the Cry1S1000 strain, among which Px006256 and Px015953 were siginificantly up-regulated for more than 4 times, while others were increased for 1∼3 times. Moreover, only four target genes (Px001355, Px002066, Px008303, and Px000596) showed the significantly lower expression levels in the Cry1S1000 strain compared to those in the G88 strain. The expression levels of Px017590 and Px007885 exhibited opposite trends with their corresponding novel-miR-240. But there was no predicted target genes for the differentially expressed novel-miR-270 and novel-miR-116 (Figure 8).

DLR Validation of the Interaction Between Novel-miR-240 and Its Target Genes
DLR assay was used to assessed whether the novel-miR-240 could regulated the expression of Px017590 and Px007885 by acting on their predicted binding sites. The target sequences (∼40 bp DNA fragments) containing the binding sites in the CDS of Px017590 and Px007885 were sythesized and cloned into the pmirGLO vector (pmirGLO-Px017590-WT and pmirGLO-Px007885-WT), respectively (Figures 9A,B). When novel-miR-240 agomir were co-transfected with pmirGLO-Px017590-WT or pmirGLO-Px007885-WT in HEK293T cells, the normalized firefly luciferase activity were siginificantly declined by 30 or 50% compared to the NC agomir control. Furthermore, we also designed the mutant fragments by altering the bases in the seed binding regions of Px017590 and Px007885 (pmirGLO-Px017590-MUT and pmirGLO-Px007885-MUT). As expected, the luciferase reporter activity was not affected by novel-miR-240 agomir, when it was co-transfected with pmirGLO-Px017590-MUT or pmirGLO-Px007885-MUT in HEK293T cells (Figures 9C,D). These results revealed that novel-miR-240 could regulate the expression of Px017590 and Px007885 by binding to their seed regions and therfore inhibiting mRNA translation.

DISCUSSION
MiRNAs, a class of endogenous non-coding RNAs, are considered to be modulators of gene expression at the posttranscriptional level (Lee et al., 2004). In recent years, with the development of the in silico miRNA identification platform such as mirMachine (Cagirici et al., 2021), miRanda (Rehmsmeier et al., 2004) and RNAhybrid (Mohebbi et al., 2021), many miRNAs have been characterized from multiple organisms including plants, vertebrates and invertebrates following the integration of various tools in bioinformatics, genetics, molecular biology and biochemistry. In insects, miRNAs played indispensable roles in development (Shen et al., 2020), behavior (Niu et al., 2019), immunity (Li R. et al., 2019), host-virus interaction (Singh, 2020), and resistance to multiple insecticides FIGURE 4 | Cluster analysis diagram showing the differentially expressed miRNAs in the resistant (Cry1S1000) and susceptible (G88) strains. Clustering was performed with log 10 (TPM + 1E-6) values. Samples of the G88 strain, S01-S03; samples of the Cry1S1000 strain, S04-S06. Columns represent different samples, while rows indicate different miRNAs. Red blocks represent the higher expressed miRNAs, green blocks represent the lower expressed miRNAs. Li X. et al., 2020). In addition, the high-quality genome offers the opportunity to explore the putative biological function of identified miRNAs (Robertson et al., 2018). DBM is a notorious agriculture pest, understanding the roles played by miRNAs in insecticides resistance can provide a significant theoretical basis for finding novel targets for pest control.
Here, we successfully constructed sRNA libraries from Bt resistant and susceptible strains, and identified a total of 437 miRNAs including 76 known and 361 novel miRNAs belonging to 91 miRNA families. Importantly, 12 differentially expressed miRNAs (three up and nine down) were revealed in the Cry1S1000 strain compared to the G88 strain. We also predicted the putative target genes for each of these 12 differentially expressed miRNAs, and the target genes were enriched in 43 GO Px000280 Elongation of very long chain fatty acids protein 7 Px016697 Putative L-ribulose-5-phosphate 3-epimerase sgbU terms based on GO and KEGG pathway analyses. Furthermore, the expression patterns of these differentially expressed miRNAs and their potential target genes were investigated with RT-qPCR. Finally, the dual luciferase assay was carried out to confirm the interaction between novel-miR-240 and its target genes Px017590 and Px007885.
From the identified miRNAs, a member of the miR-279 family, miR-279b-3p, was listed as the most abundant miRNA in our sRNA libraries. Previously, it has been found that miR-279, another member of the miR-279 family, was also listed as an abundantly expressed miRNA in DBM larvae after destruxin A injection and predicted to regulate the immunity-related genes (Shakeel et al., 2018). MiR-306, a common abundantly expressed miRNA identified in our study, has been shown to associated with Bt Cry1Ab resistance in Ostrinia furnacallis (Xu et al., 2015).
Additionally, a total of 78 known and novel miRNAs were classified into 91 miRNA families. MiRNAs sharing the same seed region were classified as a miRNA family, which likely indicated that they targeted the same genes and performed similar biological functions. In the current study, novel-miR-133 and novel-miR-250 were classified into the let-7 family, which is functionally conserved from insects to humans and participated in multiple biological processes. In insects, modulating the abundance of let-7 altered the tolerance of Myzus persicae nicotianae to nicotine . Let-7 was also required in many developmental processes, such as sleep homeostasis (Goodwin et al., 2018), age-dependent behavioral changes (Behura and Whitfield, 2010), developmental transitions of egg hatching, molting, pupation, adult eclosion  and the regulation of ecdysone levels (Liu et al., 2007).
Novel-miR-93 and novel-miR-125 were classified into the miR-14 family, which is conserved in insects and is closely related to olfactory and chemoreception (Shan et al., 2020). Another study showed that miR-14-3p, another member of the miR-14 family, could tightly control the ecdysone signaling pathway by regulating EcR and E75 (Luo et al., 2020). A total of five miRNAs (novel-miR-142, novel-miR-214, novel-miR-219, novel-miR-222, and novel-miR-253) identified in our study belong to the miR-2 family. The conserved miR-2 family was an invertebrate-specific family of miRNAs and involved in insect metamorphosis (Lozano et al., 2015), oogenesis (Song et al., 2019) neural development (Marco et al., 2012), as well as deltamethrin resistance in Cx. pipiens pallens .
Pxy-miR-252 belonged to the miR-252 family. It has been found that miR-252 could directly repress the mbt expression to control the developmental growth of Drosophila (Lim et al., 2019) and was also involved in dengue virus replication in Aedes albopictus (Yan et al., 2014). Pxy-miR-274 was classified into the miR-274 family. Inhibition of miR-274-3p could facilitate B. mori cytoplasmic polyhedrosis virus (BmCPV) replication (Wu et al., 2017). Novel-miR-78 was classified into the miR-276 family. The miR-276 family might be involved in the spirotetramat resistance of Aphis gossypii Glover  and the dengue virus (DENV) replication in C6/36 cells (Su et al., 2019). Pxy-miR-277 belonged to the miR-277 family, which was conserved in insects and participated in multiple physiological processes such as restoring immune homeostasis of the IMD pathway (Li R. et al., 2020), controlling metamorphosis (Shen et al., 2020), lipid metabolism and reproduction (Ling et al., 2017). MiR-277-3p was also a part of the molecular toolkit regulating reproductive diapause in Culex pipiens (Meuti et al., 2018) and Helicoverpa zea (Reynolds et al., 2019). Hence, we speculate that the pxy-miR-277 might be invloved in diverse biological processes in DBM and indirectly affect detoxification of Bt Cry toxin. Pxy-miR-279a, pxy-279b-5p, and pxy-279-3p were classified to the miR-279 family. It has been reported that miR-279 was essential to maintain circadian rhythm (Sun et al., 2015) and reproductive responses by females to male sex peptide in D. melanogaster (Fricke et al., 2014). So, pxy-miR-279a, pxy-279b-5p, and pxy-279-3p might play the similar roles with miR-279 in insects. Pxy-miR-306 and pxy-miR-308 belonged to miR-306 and miR-308 families, respectively. MiR-279, miR-306, and miR-308 were also predicted to be involved in the Destruxin A-responsive immunity process in DBM (Shakeel et al., 2018), indicating that these miRNA families might be closely associated with Bt resistance. Novel-miR-21 was classified into the miR-317 family, and miR-317 was found to be involved in the process of pupation in B. dorsalis . In Drosophila, miR-87 was an important regulator of dendrite regeneration (Kitatani et al., 2020), so we speculated that novel-miR-199 as a member of the miR-87 family, might paly the similar roles in DBM. In summary, these miRNA families identified in the current study were widely involved in diverse biological processes in insects, some of the miRNA families also took part in insect immunity. Investigating the functions of these conserved miRNA families and their representative members will provide a better understanding of miRNA-mediate post-transcriptional regulation in insects.
It has been previously shown that the differentially expressed miRNAs among the susceptible and resistant strains might be related to insecticides resistance in DBM. For instance, the expression of miR-276-5p and miR-8530-5p were 2-3 times higher in the chlorantraniliprole resistant DBM strain compared to the susceptible strain, indicating that these two miRNAs likely participated in the response of DBM to chlorantraniliprole (Zhu et al., 2017). The expression level of miR-278-3p was 4.2 fold lower in the deltamethrin-resistant strain (DR-strain), which induced the pyrethroid resistance in Cx. pipiens pallens (Lei et al., 2014). The other two conserved miRNAs, miR-932 and miR-285 were estimated to be 1.8-and 2.8-fold higher, respectively, in the DR-strain and were involved in pyrethroid resistance Tian et al., 2016). In this study, the expression levels of novel-miR-240, novel-miR-270 and novel-miR-116 in the resistant Cry1S1000 strain were significantly lower than that in the susceptible G88 strain, other miRNAs such as novel-miR-210, novel-miR-274, novel-miR-288, novel-miR-25, and pxy-miR-8522 represented significantly higher expression levels in the resistant Cry1S1000 strain, and novel-miR-48, novel-miR-97, novel-miR-237, and novel-miR-225 exhibited no significant difference between these two strains. Moreover, the differentially expressed novel-miR-240, novel-miR-270, and novel-miR-116 did not belong to any known miRNA family, indicating that there was no miRNA shared the same seed regions or derived from the same arm of the stem-loop with these two miRNAs (Meyers et al., 2008). While novel-miR-270 had a similar sequence with ame-miR-9864-5p, suggesting that they may be the same miRNA in different species and play similar roles in insects. Our results suggested that the expression of miRNAs varies among different FIGURE 6 | The most enriched KEGG pathways based on the genes targeted by differentially expressed miRNAs in the resistant (Cry1S1000) and susceptible (G88) strains. The X-axis shows the number of genes annotated to the pathway and the ratio of the number of annotated genes to the total genes annotated, and the Y-axis shows the pathway names. The different color of column indicates different type of KEGG pathway.
FIGURE 7 | Relative expression levels of differentially expressed miRNAs in midguts of the resistant (Cry1S1000) and susceptible (G88) strains by RT-qPCR (A), and up-/down-regulated expression of the miRNAs based on Illumina sequencing (B). The expression levels of miRNAs were normalized by U6. Statistical significance was analyzed using one-way ANOVA. The asterisks represent significance, where two asterisks indicate p < 0.01, and one asterisk indicates p < 0.05. (B) The orange column above the X-axis shows up-regulated miRNAs, while the purple column below the X-axis shows down-regulated miRNAs based on Illumina sequencing. strains, and these differentially expressed miRNAs may regulate detoxification to Bt Cry1Ac toxin in a complicated manner by targeting detoxification-related genes.
It should be noted that novel-miR-240 represents downregulated expression in the Cry1S1000 strain relative to the G88 strain, whereas its target genes Px017590 and Px007885 were upregulated in the Cry1S1000 strain. The DLR assay results implied the involvement of novel-miR-240 in the regulation of Px017590 and Px007885, indicating that this miRNA might modulate Bt resistance of DBM via negatively regulating its target genes.
Furthermore, GO annotation exhibited that most of the predicted target genes of differentially expressed miRNAs were enriched in biological process (cellular process, metabolism process, single-organism process, and biological regulation), cellular component (membrane, cell, cell part, and membrane part) and molecular function (binding and catalytic activity), which might be related to the insect immunity and the metabolism to the Bt Cry1Ac toxin. The target genes of differentially expressed miRNAs were mostly enriched in metabolic-related pathways such as Hippo signaling pathway, minoacyl-tRNA biosynthesis and MAPK signaling pathway. And in a previous study, MAPK signaling pathway was found to lead to Bt Cry1Ac resistance by mediating differential expression of APN and other midgut genes in DBM . The DBM-DB database was used for the annotation of target genes. Px015953 and Px000596, the target genes of novel-miR-25, were predicted to be associated with chitinase and cuticle protein, which may be involved in the process of DBM combating the exogenous Bt Cry1Ac toxin.
Chitinase, an insect molting enzyme, is a potential biopesticide that degrades chitin to low molecular weight, soluble and insoluble oligosaccharides. Previous research has shown that FIGURE 8 | Relative expression levels of target genes in midguts of the resistant (Cry1S1000) and susceptible (G88) strains based on RT-qPCR. The expression levels of target genes were normalized by RPL32. Statistical significance was analyzed using one-way ANOVA. The asterisks represent significance, where three asterisks indicate p < 0.001, two asterisks indicate p < 0.01, and one asterisk indicates p < 0.05. chitinases could degrade the vital structures of insects such as their peritrophic membranes and cuticles, thus facilitating the entry of the pathogens into the tissues of susceptible insects (Brandt et al., 1978). Chitinases also facilitated the penetration of host cuticle by entomopathogenic fungi. Analysis of the top 20 differentially expressed genes showed that the chitinase ChiA2 was up-regulated after Beauveria bassiana infected nonnatural hosts such as H. armigera and Clanis bilineata (Liu et al., 2021). These chitinase-related genes indicated that when entomopathogenic fungi transfer to non-original hosts, the fungi changed the metabolic response of hosts and used the novel infection strategies to break the barrier of different cuticle chitin components to better adaption to new hosts. Therefore, the chitin based bioformulation may act as an effective pesticide against Lepidopterous pests. A prior work revealed that the expression of cuticle proteins were siginificantly induced at 6 h posttreatment of Cry1Ac toxin (Qin et al., 2021). Genes encoding cuticle proteins were induced during B. mori bidensovirus (BmBDV) infection . The down-regulation of some Cardinium-responsive miRNAs enchanced the cuticle FIGURE 9 | Dual luciferase reporter assay confirmed the interaction between novel-miR-240 and its target genes in vitro. (A) Predicted target sites of novel-miR-240 in the CDS of Px017590. (B) Predicted target sites of novel-miR-240 in the CDS of Px007885. WT, wild-type sequence; MUT, mutant sequence. (C) The interaction between novel-miR-240 and Px017590. (D) The interaction between novel-miR-240 and Px007885. Statistical significance was analyzed using one-way ANOVA followed by Tukey's multiple comparisons. The three asterisks indicate p < 0.001 and "ns" represents no significance. proteins in whitefly Bemisia tabaci (Gennadius) biotype Q and modulated insect reproduction and growth . In addition, RNA-seq of Cryptolestes ferrugineus individuals from phosphine susceptible and resistant populations revealed the significantly increased expression levels of nine cuticular protein genes in the resistance population, demonstrating that the cuticular protein genes might be related to phosphine resistance (Chen et al., 2021). In the current study, the target genes of identified differentially expressed miRNAs, especially the chitinase and cuticle protein, might involved in the insects immune response and the detoxification of xenobiotics including the Bt Cry1Ac toxin. Although the RT-qPCR results showed that some of the identified differentially expressed miRNAs exhibited positively correlate with their potential targets, the target genes of these differentially expressed miRNAs were still worthy for further investigation, where functional characterizations of the target genes will provide novel insights of insects immunity and detoxification mechanisms.

CONCLUSION
In the present study, we identified 76 known and 361 novel miRNAs from midguts of the Bt resistant and susceptible DBM strains. Among the identified miRNAs, 178 were classified into 91 miRNA families, and the potential functions of these miRNA families have also been illustrated. Expression of 12 differentially expressed miRNAs and their potetial target genes were profiled, revealing that these miRNAs may play vital roles in detoxicating Bt Cry1Ac toxins by regulating their corresponding target genes. Moreover, the interaction between novel-miR-240 and its target genes Px017590 and Px007885 have been confirmed by the DLR assay. However, further experiments are still required to demonstrate the functional pathways which the differentially expressed miRNAs and their target genes are involved in. Our study contributes to the field for further investigation of the miRNA-mediated gene regulation at the post-transcriptional level in insect resistance and the RNAi-based pest management.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI SRA; PRJNA754494.