Global identification of microRNAs associated with chlorantraniliprole resistance in diamondback moth Plutella xylostella (L.)

The diamondback moth (DBM), Plutella xylostella (L.), is one of the most serious cruciferous pests and has developed high resistance to most insecticides, including chlorantraniliprole. Previous studies have reported several protein-coding genes that involved in chlorantraniliprole resistance, but research on resistance mechanisms at the post-transcription level is still limited. In this study, a global screen of microRNAs (miRNAs) associated with chlorantraniliprole resistance in P. xylostella was performed. The small RNA libraries for a susceptible (CHS) and two chlorantraniliprole resistant strains (CHR, ZZ) were constructed and sequenced, and a total of 199 known and 30 novel miRNAs were identified. Among them, 23 miRNAs were differentially expressed between CHR and CHS, and 90 miRNAs were differentially expressed between ZZ and CHS, of which 11 differentially expressed miRNAs were identified in both CHR and ZZ. Using miRanda and RNAhybrid, a total of 1,411 target mRNAs from 102 differentially expressed miRNAs were predicted, including mRNAs in several groups of detoxification enzymes. The expression of several differentially expressed miRNAs and their potential targets was validated by qRT-PCR. The results may provide important clues for further study of the mechanisms of miRNA-mediated chlorantraniliprole resistance in DBM and other target insects.

Sequencing of miRNAs from P. xylostella. Three small RNA (sRNA) libraries for CHS, CHR and ZZ were constructed, and 45,297,627, 47,566,639 and 40,906,561 raw reads were generated from the sRNA library for CHS, CHR and ZZ, respectively. The low-quality sequences, reads without a 3′ adaptor and reads that were less than 18 nt were eliminated; subsequently, 31,952,995, 33,561,726, and 29,46,332 clean reads were respectively obtained for the three strains and used for further analysis ( Table 2). The length distributions are displayed in Identification of known miRNAs. Before this research, 235 miRNAs and 462 miRNAs in P. xylostella had already been identified by Etebari et al. 21 and Liang et al. 22 , respectively, in 2013. However, a number of same or similar miRNAs in these two publications were named differently. Therefore, we collated all the miRNAs reported in these two references first, and then all clean sequences generated from this study were mapped against the resulting miRNAs. As a result, 199 known mature miRNA were obtained. Then the pre-miRNA sequences of these mature miRNAs were aligned to those reported by Etebari et al. 21 and Liang et al. 22 and mapped to the latest version of P. xylostella genome, and 121 confident pre-miRNA sequences were conformed, which produced 172 of 199 identified mature miRNAs (Table S1). However, the pre-miRNA sequences of the rest 27 conserved miRNAs were not detected in the current P. xylostella genome. Considering the incomplete assembly of this DBM genome version, these 27 conserved miRNAs were also used for further analysis (Table S2).
Identification of novel miRNAs. After removal of the identified known miRNA sequences mentioned above, the rest of the clean sequences were processed to remove any non-coding RNAs, protein-coding RNA fragments and repeated sequences. Finally, 20,801,228 (CHS), 22,526, 217 (CHR) and 19,823,507 (ZZ) unannotated clean sequences were obtained and used to predict novel miRNAs ( Table 2).
According to the unannotated sequences, miRNA precursor sequences and structures were predicted and identified using miRDeep2 (a probabilistic algorithm based on the miRNA biogenesis model) 23 , and a total of 51 potential novel miRNAs were initially predicted from the three libraries (Table S3). RNAfold was used to confirm the structures of the predicted miRNAs 24 . After keeping only novel miRNAs with a rand fold P-value ≤ 0.05 and a miRDeep2 score ≥ 3, 30 potential novel microRNAs were retained and used for expression analysis, and 24 were determined to have complementary star miRNA sequences (most of the star miRNAs have a low copy number). The length of the novel miRNAs ranged from 18 to 25 nt. Novel miRNAs were named based on their positions in the P. xylostella genome (Table S3). Secondary structures for some potential miRNA precursors with high miRD-eep2 scores are shown in Fig. 2. The most abundant miRNAs in P. xylostella. MiR-Bantam, miR-10 and miR-281 were the three most abundant miRNAs in each of the three libraries. These three miRNAs were also highly expressed in P. xylostella second instar larvae reported by Etebari et al. 21 . Ten of the most highly expressed miRNAs are listed in Table 3, and a total of 20 miRNAs had high expression with a mean count number > 10,000, including 2 novel miRNAs, pxy-novel-117_9740 and pxy-novel-95_8740 (Tables S1, S2, S3).

Analysis of differentially expressed miRNAs.
To systematically identify chlorantraniliprole resistance associated miRNAs, a differential expression analysis was performed among the three strains using the sequencing results. In total, 20 known and 3 novel miRNAs were identified as differentially expressed between CHR and CHS; 13 were significantly down-regulated, and 10 were significantly up-regulated (Table S4, Fig. 3A). In addition, 80 known and 10 novel miRNAs were differentially expressed between ZZ and CHS, 89 were significantly down-regulated, and only 1 were up-regulated (Table S5, Fig. 3B). Compared to CHS, 9 known and 2 novel miRNAs were found to be differentially expressed in both CHR and ZZ, 10 were down-regulated in both resistant strains, except pxy-miR-8491-5p, which was up-regulated ( Fig. 3C, Table 4). Overall, most of the differentially expressed miRNAs were down-regulated in the resistant strains.
Target prediction and annotation of differentially expressed miRNAs. Usually, miRNA functions by binding to its target mRNAs, therefore annotating the potential targets of differentially expressed miRNAs, and is very important in defining their roles in chlorantraniliprole resistance. For 23 significantly differentially expressed miRNAs between CHS and CHR, a total of 5,384 targets were predicted using miRanda, and 5,241 targets were identified using RNAHybird (Table S6, Fig. 4A). Similarly, for 90 significantly differentially expressed miRNAs between CHS and ZZ, 23,129 and 20,111 targets were predicted with miRanda and RNAHybird, respectively (Table S7, Fig. 4B).
To make the prediction results more reliable, the miRNA targets predicted by both miRanda and RNAhybrid were considered the final target genes. Finally, 242 miRNA-mRNA pairs between CHS and CHR (Table S8, Fig. 4A), and 1,276 miRNA-mRNA pairs between CHS and ZZ were supported by both algorithms (Table S9, Fig. 4B). Furthermore, all of the binding sites predicted by miRanda and RNAhybrid between miRNAs and their mRNA targets were counted. For different expressed miRNAs and their potential target mRNAs between CHS and CHR, 24% of the binding sites predicted by the two algorithms were almost identical (0 nt difference), and only 7% showed shifts more than 5 nt (Fig. 4C). Similarly, for different expressed miRNAs and their potential target mRNAs between CHS and ZZ, 23% of the binding sites predicted by the two algorithms were almost identical, and 9% showed shifts more than 5 nt (Fig. 4D).
A set of miRNAs were found to target several families of important genes that are often involved in insecticide resistance, such as cytochrome P450 [25][26][27][28] , esterase 29 , GSTs 30 , ABC transporter family protein 31,32 , cuticle protein 33 , glutamate-gated chloride channel 34,35 and superoxide dismutase (SOD) 36 . The target genes of some selected miRNAs are listed in Table 5. The final predicted targets for the differentially expressed miRNAs were used for qRT-PCR analysis, and their GO annotations and KEGG pathway mapping results are listed in Table S8 and Table S9, respectively.

Quantitative RT-PCR validation of differentially expressed miRNAs and their potential targets.
To validate the expression profiles of the differentially expressed miRNAs identified from the small RNA sequencing, a number of miRNAs were randomly selected for quantitative RT-PCR (qRT-PCR) assays.
To analyze the correlation between the expression levels of miRNAs and their potential targets, the relative expression of three miRNAs (pxy-miR-8533-3p, pxy-miR-8534-5p and pxy-miR-375-5p) down-regulated in both CHR and ZZ and their corresponding targets, including larval cuticle protein LCP-30, cytochrome P450 6B6 and cytochrome P450 4G15, were verified through qRT-PCR. The expression of the 3 selected miRNAs were all down-regulated, which shared a similar trend with the sequencing results, while the expression of their corresponding targets were all up-regulated (Fig. 8). All 3 selected miRNAs showed a significant negative correlation with their targets.

Discussion
In recent years, an increasing number of papers on insect miRNA have been published 37 . A total of 3,824 mature miRNAs belonging to 26 species of insects have already been deposited in miRBase. MiRNA plays very important roles in the growth and development of insects, such as the reproductive process 38 germ cell development 39 neurogenesis 40,41 , wing development [42][43][44] , phenotypic plasticity 45 and muscle growth [46][47][48] . More recently, miRNAs were also found to be involved in insecticide resistance [49][50][51] . To systematically identify miRNAs associated with chlorantraniliprole resistance in P. xylostella, small RNAs from a susceptible strain (CHS) and two resistant strains (CHR, ZZ) were sequenced using Illumina sequencing technology in this study. The differentially expressed miRNAs among CHS, CHR and ZZ were analyzed, and their target genes were also predicted. A total of 229 miRNAs were identified in the three libraries, of which 199 miRNAs had already been reported in P. xylostella before, and 30 novel miRNAs were predicted using miRdeep2 for the first time.
Some highly conserved miRNAs, such as miR-let-7, miR-8, miR-9, miR-184, miR-278 and miR-bantam, which play essential roles in many types of insects 37 , were also discovered with high expression in the three P. xylostella strains, implying important regulatory roles in P. xylostella.
Twenty-three miRNAs were identified to be differentially expressed between CHR and CHS. Because CHR was established from CHS by successive selection with chlorantraniliprole (i.e., they have same genetic background) and have been reared under same laboratory conditions with CHS, the 23 differentially expressed miRNAs are likely to be associated with chlorantraniliprole resistance. Between ZZ and CHS, 90 differentially expressed miRNAs were identified. The ZZ strain was a field strain, and it had developed high levels of resistance to several other commonly used insecticides, such as beta-cypermethrin, abamectin, spinosad and indoxacarb (unpublished data from a local plant protection station), in addition to chlorantraniliprole. Each of these insecticides kills P. xylostella with distinctive modes of action. Therefore, the 90 differentially expressed miRNAs likely result from of the comprehensive effects of these different insecticides as well as other environmental factors.
When the differentially expressed 23 and 90 miRNAs were put together, we found that 11 of them overlapped. The overlapped miRNAs are likely to be involved in chlorantraniliprole resistance because they were differentially expressed in both laboratory-selected and field-collected resistant strains. However, due to the complex insecticide resistance mechanisms in the ZZ strain, the 11 miRNAs may also reveal common resistant mechanisms to other insecticides. In fact, some of them have been reported to be associated with several insecticides Figure 2. Predicted secondary structure of three selected novel miRNAs. The entire sequence represents pre-miRNAs, the red represents mature miRNA, and the purple represents miRNA*.   49 . In this study, we found that an analog of miR-375 was differentially expressed between chlorantraniliprole susceptible and the two resistant P. xylostella strains, though with low abundance. These results, together with one P450 mRNA, imply that miR-375 has a high possibility of involvement in the regulation of insecticide resistance. The unique differentially expressed miRNAs in the CHR strain may reveal a unique mechanism of resistance to chlorantraniliprole; therefore, we should pay more attention to these 10 miRNAs in our follow-up work. The unique differentially expressed miRNAs in the ZZ strain, especially those that have different expression profiles compared with CHR, are more likely to be involved in other insecticide resistance and not chlorantraniliprole resistance. Some of them have been reported to be associated with deltamethrin, Cry1Ab or fenpropathrin resistance, such as miR-210, miR-965, miR-981, miR-1, miR-306 and miR-281 [49][50][51] . Because experimental validation of miRNA targets is still a major challenge, in silico prediction is widely used for the identification of potential miRNA targets 52 . In the current study, the target genes of the differentially expressed miRNAs among CHS, CHR and ZZ were predicted using two different algorithms, miRanda and RNAhybrid. Only the target genes supported by both algorithms were retained, and thus, the predicted results were relatively more credible. Although we did not find a ryanodine receptor in the predicted results, many other  insecticide resistance associated genes were discovered, such as P450, esterase, GSTs, ABC transporter family protein, cuticle protein and SOD genes, of which several cytochrome P450, multidrug resistance-associated protein 4 (ABCC4) and SOD genes have already been identified to be differentially expressed among the susceptible strain and 3 different chlorantraniliprole-resistant P. xylostella strains 36 . Esterase FE4 29 , GSTT1 30 and larval cuticle protein LCP-30 genes 33 were also confirmed to be involved in resistance to other insecticides. In addition, some other metabolic enzymes were predicted to be associated with the differently expressed miRNAs, but there was no evidence for their involvement in insecticide resistance yet. Interestingly, several of them were discovered to be involved in detoxification process in mammal, such as UDP-glucuronosyltransferase 2A1 53 .

MiRNA
In this study, a number of miRNAs were speculated to be involved in chlorantraniliprole resistance in P. xylostella based only on their differential expression profiles and their predicted target genes. To further reveal the roles of these miRNAs in chlorantraniliprole resistance, gain-and loss-of-function experiments should be carried out in vivo and could include, for example, the up-or down-regulation of the expression of a miRNA by injecting its synthetic mimics or inhibitors and suppressing the expression of its target mRNAs by RNA interference (RNAi). This should be followed by the evaluation of resistance levels in treated larvae using bioassays.

Conclusion
This paper presents the first study to systematically screen miRNAs associated with chlorantraniliprole resistance in P. xylostella. In this study, we identified 199 known miRNAs and 30 novel miRNAs in three DBM strains. A set of differentially expressed miRNAs among CHS, CHR and ZZ were obtained and considered highly likely to be associated with chlorantraniliprole resistance. The potential targets of the differentially expressed miRNAs were also predicted, and many of the target genes were related to detoxification processes. These results may guide us in further investigating the mechanisms of miRNA-regulated chlorantraniliprole resistance and may provide novel insights into resistance management in P. xylostella.

Materials and Methods
Insects. The laboratory susceptible DBM strain (CHS) was collected in the vegetable fields of Beijing and reared in our laboratory without exposure to any insecticide for more than 10 years. The chlorantraniliprole-resistant strain (CHR) was derived from the CHS strain by successive selection with chlorantraniliprole for more than 60 generations, and the field-resistant strain (ZZ) was collected in the vegetable fields of Zhangzhou, Fujian province, in southeastern China in 2015. All stages of P. xylostella were maintained at 27 ± 1 °C, an RH of 50-60% and a photoperiod of 16 h light/8 h dark on radish seedlings (Raphanus sativus L.). P. xylostella adults were provided with 10% (W/V) honey solution and allowed to mate and oviposit on the radish seedlings.
Bioassay. The Leaf-dip method 54 was used in this study. Cabbage leaves were dipped in the required chlorantraniliprole concentrations for 10-15 s and then allowed to air dry in the shade. A 0.1% (v/v) Triton X-100 water solution was used as a control. Approximately 20-25 third instar larvae were transferred onto each leaf, and three replications were used for each concentration. The mortality was assessed after four days of treatment. LC 50 values were calculated using POLO-Plus 2.0 software (LeOra Software Inc., Berkeley, CA).

MiRNA
Target id Annotation Related insecticide References  Table 5. Potential target genes of miRNAs that have already been identified to be involved in insecticide resistance.

Figure 5. qRT-PCR validation of significantly differentially expressed miRNAs between CHS and CHR.
Different lowercase letters (a and b) represent significant differences by t-test (P < 0.05). The same applies below.
Small RNA library construction and sequencing. RNA samples were prepared from the three DBM strains (each sample containing 30-50 third-instar P. xylostella larvae). Trizol Reagent was used to isolate the total RNA from each sample according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA). RNA   Analysis of sequencing data. Data files from each of the three libraries were used for analysis. The raw data of this study were deposited in the NCBI Short Read Archive (SRX1968414, SRX1968415 and SRX1952874). Clean reads were screened from the raw data after processing out low-quality reads, adapters, and sequences of fewer than 18 nucleotides. Clean reads were mapped to the P. xylostella genome 55 (version 2, http://iae.fafu.edu. cn/DBM/index.php) to analyze the distribution using bowtie software 56 . All P. xylostella miRNAs reported by Etebari et al. 21 and Liang et al. 22 were collated and re-annotated first, then clean sequences generated in this research were used to search against the resulting miRNAs. All miRNAs identified in this step were considered known miRNAs.
All remaining clean sequences were subjected to the Rfam database to remove known noncoding RNA families, including rRNA, scRNA, snoRNA, snRNA and tRNA (http://www.sanger.ac.uk/Software/Rfam/ftp.shtml), and were then searched against known genes of P. xylostella to discard degraded fragments; the search was concluded in RepBase to remove repeated sequences (http://www.girinst.org). Unannotated clean sequences that did not match any of the above databases were further used to analyze and predict novel miRNAs using miRDeep2 software. Both a false positive rate (FPR) and true positive rate (TPR) were used to assess the predicted results, and RNAfold was used to confirm the structures of the predicted miRNAs. Differential Expression Analysis of miRNAs. The abundance of miRNAs identified in the three libraries was first normalized using the tags per million reads (TPM) method: TPM = (number of mapped reads for each miRNA/total number of mapped reads) × 10 6 . The log 2 (TPM ratios) among the three libraries was calculated, and the P-value was calculated using the Audic Claverie statistic. The miRNAs with |log 2 (TPM ratios)| ≥ 1 and P-value < 0.05 were regarded as differentially expressed among the three P. xylostella strains.
Target prediction and annotation of differentially expressed miRNAs. MiRNA usually regulates gene expression through binding to the 3′ untranslated region (3′ UTR) of target mRNAs, so 3′ UTR annotation information was first extracted from the genome database of P. xylostella for target prediction. The potential target genes of differentially expressed miRNAs were predicted and analyzed using two different types of software, miRanda 57 and RNAhybrid 58,59 . For each prediction method, high efficacy targets were selected by the following criteria: (1) miRanda: total score ≥ 140, total energy ≤ − 25 kcal/mol; (2) RNAhybrid: P-value < 0.05, mfe ≤ − 25 kcal/mol. Verification of differentially expressed miRNAs and their potential targets by quantitative real-time PCR. Quantitative real-time PCR was performed to experimentally validate the relative expression levels of the identified miRNAs and their potential targets. Total RNA was extracted from the same samples used for deep sequencing. The first-strand cDNA of mature miRNA and mRNA were synthesized using a miScript II RT kit (Qiagen, Germany) and a PrimeScript ™ RT reagent Kit with gDNA Eraser (Perfect Real Time) (Takara Biotechnology, Dalian, China), respectively, following the manufacturer's instructions. qRT-PCR analysis was carried out using SYBR Premix Ex Taq (Takara Biotechnology, Dalian, China). Each reaction was performed in an ABI 7500 Real Time PCR system (Applied Biosystems) with three biological replicates. The expression levels for miRNA and mRNA were normalized to U6 snRNA and ribosomal protein L32 mRNA, respectively. The relative expression levels of the miRNAs and targets were calculated using the 2 -ΔΔCt method 60 . All primers used in this study are listed in Table S10.