A systemic study of indoxacarb resistance in Spodoptera litura revealed complex expression profiles and regulatory mechanism

The tobacco cutworm, Spodoptera litura, is an important pest of crop and vegetable plants worldwide, and its resistance to insecticides have quickly developed. However, the resistance mechanisms of this pest are still unclear. In this study, the change in mRNA and miRNA profiles in the susceptible, indoxacarb-resistant and field indoxacarb-resistant strains of S. litura were characterized. Nine hundred and ten co-up-regulated and 737 co-down-regulated genes were identified in the resistant strains. Further analysis showed that 126 co-differentially expressed genes (co-DEGs) (cytochrome P450, carboxy/cholinesterase, glutathione S-transferase, ATP-binding cassette transporter, UDP-glucuronosyl transferase, aminopeptidase N, sialin, serine protease and cuticle protein) may play important roles in indoxacarb resistance in S. litura. In addition, a total of 91 known and 52 novel miRNAs were identified, and 10 miRNAs were co-differentially expressed in the resistant strains of S. litura. Furthermore, 10 co-differentially expressed miRNAs (co-DEmiRNAs) had predicted co-DEGs according to the expected miRNA-mRNA negative regulation pattern and 37 indoxacarb resistance-related co-DEGs were predicted to be the target genes. These results not only broadened our understanding of molecular mechanisms of insecticide resistance by revealing complicated profiles, but also provide important clues for further study on the mechanisms of miRNAs involved in indoxacarb resistance in S. litura.

Analysis of mRNA sequencing data. In order to identify the mRNA expression profiles in different strains of S. litura, a total of three libraries, SS, InRS and FInRS, were constructed and sequenced by Illumina HiSeq 2500 platform. The samples from each strain were sequenced in triplicate. The SS, InRS and FInRS libraries were found to contain 163,943,684, 156,829,302 and 177,060,140 raw reads, respectively. After removing low-quality reads, adaptors and all possible contaminants, 161,116,162, 153,986,650 and 173,406,844 clean reads were obtained, respectively. Among them, 89.94%, 89.81% and 86.89% clean reads were uniquely mapped to the reference genome, respectively ( Table 2).
The DEGs were screened based on DEGSeq 2 analysis taking |log 2 (fold change)| ≥ 1 and P < 0.05 as the cut-off between resistant and susceptible strains. Comparing with the SS, 1576 up-and 1227 down-regulated, 2606 upand 1781 down-regulated genes were identified in the InRS and FInRS of S. litura, respectively (Fig. 1A). Venn diagram was generated using DEGs and it was depicted that 910 genes were co-up-regulated and 737 genes were co-down-regulated in the InRS and FInRS of S. litura (Fig. 1B). A cluster heat map was subsequently adopted to show co-DEGs in the InRS and FInRS of S. litura. co-DEGs were divided into two groups with four clusters, and the gene expression of the InRS is similar to that of the FInRS, but they were differentially expressed compared to the SS (Fig. 1C). Furthermore, GO analysis results showed that these co-DEGs were divided into three ontologies: biological process, cellular component, and molecular function, including 26 GO terms, and the most enriched functions were structural molecule activity, transmembrane transport, extracellular region, cofactor binding and lipid metabolic process ( Supplementary Fig. S1).  Table S4).   Table 3. Small RNA sequencing statistics.
Through the bioinformatics analysis, 91 known miRNAs and 52 novel miRNAs were identified in the SS, InRS and FInRS of S. litura (Supplementary Table S5). A total of 7 up-and 9 down-regulated, 14 up-and 21 down-regulated miRNAs were detected in the InRS and FInRS of S. litura, respectively ( Fig. 2A). Venn diagram showed that 5 miRNAs were co-up-regulated and 5 miRNAs were co-down-regulated in the InRS and FInRS of S. litura (Fig. 2B). The cluster heat map results of co-DEmiRNAs depicted that there are much more differences in miRNAs expressed in the SS, InRS and FInRS of S. litura (Fig. 2C).

Validation of expression profiles by qPCR.
To validate RNA sequencing results, we used qPCR to investigate the relative expression levels of randomly selecting 20 co-DEGs and all 10 co-DEmiRNAs. The results revealed that 19 of these co-DEGs (95.0%) (Fig. 3) and 9 of these co-DEmiRNAs (90.0%) (Fig. 4) were consistent  www.nature.com/scientificreports www.nature.com/scientificreports/ with RNA sequencing. The qPCR results validated the RNA sequencing results and increased the accuracy and reliability of the differentially expressed genes and miRNAs.
DEGs Encoding P450, CCE and GST. When considering the metabolic enzyme genes potentially involved in resistance, the strong response of the resistant strains against indoxacarb selection through transcription level modifications was confirmed, with several detoxification genes being over-expressed (Fig. 5). Studying P450 genes expression revealed that 24 P450 genes showed significant transcription level variations (19 genes co-up-regulated and 5 genes co-down-regulated) in the InRS and FInRS ( Table 4). The differently expressed P450 genes distributed in all P450 clans (clan 2, 3, 4, and M) in the InRS and FInRS. The fold change of these up-regulated P450 genes ranged from 2.11-fold (SlituP450-003) to 66.46-fold (SlituP450-018) and 2.78-fold (SlituP450-037) to 199.98-fold (SlituP450-085) in the InRS and FInRS, respectively (Table 4).
Another important metabolic enzyme related to insecticide resistance is CCE. The number of differently expressed CCE genes was less than the number of P450 genes. Comparing with the SS, there were 11 CCE genes having differential expression (7 genes co-up-regulated and 4 genes co-down-regulated) ( Table 4). Three of the up-regulated CCE genes belonged to lepidopteran esterase, one of these genes belonged to α-esterase and three of these genes belonged to integument esterase ( Table 4). The fold change ranged from 3.07-fold (SlituCOE009) to 37.20-fold (SlituCOE073) and 2.72-fold (SlituCOE111) to 10.01-fold (SlituCOE073) in the InRS and FInRS, respectively (Table 4).
GST was also an important detoxification enzyme playing crucial roles in insecticide metabolic resistance. In this study, there were 3 GST genes differentially expressed (2 genes co-up-regulated and 1 gene co-down-regulated) in the InRS and FInRS ( Table 4). The fold change of these two up-regulated GST genes (SlituGST20 and SlituGST38) were 4.24-, 72.65-fold and 29.47-, 17.30-fold in the InRS and FInRS, respectively, all belonging to class epsilon. (Table 4).
Other DEGs related to insecticide resistance. There are several other insecticide resistance-related genes that were up-regulated in the InRS and FInRS, such as 10 UDP-glucuronosyl transferase (UGT), 9 ATP-binding cassette (ABC) transporter and 1 aminopeptidase N (APN) genes. These genes were detected to be co-up-regulated in the InRS and FInRS (Supplementary Table S6).
There is also a sialin gene with 2.96-fold and 2.99-fold over-expressed in the InRS and FInRS, respectively (Supplementary Table S6). Sialin belongs to the anion/cation symporter (ACS) family, which is a large subfamily of the major facilitator superfamily (MFS) of transporters. Another family of up-regulated genes associated with insecticide resistance is serine protease (SP). Six SP genes were identified as over-expressed genes, with fold change ranged from 2.27-to 5.75-fold and 3.03-to 13.33-fold in the InRS and FInRS, respectively (Supplementary Table S6). In addition, a number of genes encoded cuticle protein (CP) were over-transcribed in the InRS and FInRS with fold change ranged from 2.97-to 674.37-fold and 3.20-to 872.19-fold among the 52 up-regulated CP genes (Supplementary Table S6). miRNAs target genes prediction and correlation analysis of miRNAs-mRNAs. According to the sequence information of the known and novel miRNAs, miRNAs target genes prediction were conducted. As a result, a total of 16,235 target genes were predicted in our study (Supplementary Table S7). Further analysis showed that 746 co-DEGs were predicted as all 10 co-DEmiRNAs target genes (Supplementary Table S8). To explore the potential function of target co-DEGs, GO annotation enrichment analysis was conducted. The result of GO annotation enrichment showed that most of these predicted target co-DEGs were mainly focused on structural molecule activity, structural constituent of cuticle and oxidoreductase activity ( Supplementary Fig. S2).

Discussion
Indoxacarb is a novel oxadiazine insecticide which has good field activity against a number of lepidopteran pests, as well as certain homopteran and coleopteran pests 18 . Indoxacarb can be metabolized by insect esterases or amidases to a N-decarbomethoxylated metabolite (DCJW), which is a more active sodium channel blocker than indoxacarb, resulting in paralysis and death of target pest species 19,20 . Because of its safety to mammals and non-target organisms, favorable environmental and residue properties, broad spectrum and rapid inhibition of insect feeding making indoxacarb a powerful new insecticide for crop protection 18 . However, due to intensive use of indoxacarb, many studies have shown that several insects have developed resistance to indoxacarb in recent years, including Choristoneura rosaceana 21 , Musca domestica 22 , Plutella xylostella 23,24 , Spodoptera exigua 25 and Helicoverpa armigera 26    www.nature.com/scientificreports www.nature.com/scientificreports/ The mechanisms mediating resistance to indoxacarb have been studied in several insects. Some mutations in the sodium channel gene have been shown to confer target site insensitivity to the neurotoxic effects of indoxacarb. Wang et al. 29 identified two point mutations (F1845Y and V1848I) in P. xylostella, and Gao et al. 25 identified one point mutation (L1014F) in S. exigua, which have been proven to confer high levels of resistance to indoxacarb. Even though enhanced target site insensitivity is important for the resistance of insects, detoxification enzymes are also important factors for the metabolism of insecticides. Shono et al. 22 showed that P450 was involved in the resistance to indoxacarb in M. domestica. Sayyed and Wright 30 , and Nehare et al. 31 found that esterase and GST were related to resistance to indoxacarb in P. xylostella. Gao et al. 25 also suggested that carboxylesterase and GST were major factors leading to indoxacarb resistance in S. exigua. What is more, it was demonstrated that the increased activities of carboxylesterase and P450 were important in conferring indoxacarb resistance in S. litura 28 . However, the indoxacarb resistance mechanism of S. litura remains unclear at molecular level.
In addition to the well-known detoxification gene families involved in insecticide resistance, we also revealed that other insecticide-related genes had significant higher expression levels in the resistant strains. UGT as biotransformation enzymes, widely distributed within living organisms and viruses, were presumed to originally participate in the detoxification process 37 . Overexpressed P. xylostella UGT2B17 38 and T. cinnabarinus UGT201D3 39 have been shown to be involved in insecticide resistance. In addition to enzymes related to metabolism and conjugation, a number of transporter families, of which ABC transporters are the best studied, also play an important role in xenobiotic tolerance 40 . Sun et al. 41 suggested that ABC transporters might be involved in resistance to multiple insecticides in Laodelphax striatellus. What is more, APN has been shown to function as Cry protein receptor in insects 42 and might be involved in the response to different classes of xenobiotics in S. litura 10 . In this case, the over-expression of 10 UGT, 9 ABC transporter and 1 APN genes may also associate with the detoxification process of indoxacarb in S. litura.
By revealing that several other genes with a broad range of biological functions are similarly affected by insecticides, our results suggest that the ability of S. litura to better tolerate insecticides might also be the consequence of the induction of other proteins involved in a wide range of functions. The major facilitator superfamily (MFS), belonging to secondary active membrane transporters, can transport a wide range of small solutes (including  www.nature.com/scientificreports www.nature.com/scientificreports/ inorganic ions, sugars, amino acids, and xenobiotics) in response to chemiosmotic ion gradients in humans. The roles of MFS on adapting stress from host plant shift and acaricide exposure were reported in Tetranychus urticae 43 . Sialin is a member of the MFS of transporters 44 . In the present study, up-regulation of sialin gene in S. litura may result in a higher efflux of insecticides out of S. litura cells and developing of insecticide resistance. Serine protease, secreted by pancreas in mammals, mainly function in two aspects: protein digest and hydrolysis as well as activation of all proenzymes secreted by pancreas 45 , which is important in the cellular and the humoral arm of invertebrate immune response. It had been reported that up-regulated SP genes were related to deltamethrin resistance in C. pipiens pallens 46 . Intracellular proteases might play a role in protein biosynthesis or modified conformation of enzymes as part of this induction process 47 . This mechanism may involve increased supply of precursor amino acids from proteolytic degradation products to the intracellular pool, prior to synthesis of detoxifying enzymes in S. litura following insecticide exposure. This hypothesis could be supported by another fact that 6 up-regulated SP genes were found in S. litura responded to indoxacarb in our study. Many studies on cuticle protein suggest that insecticide resistance could manifest as a slower rate of penetration due to higher protein and lipid content in the cuticle and/or greater sclerotization and was caused by increased cuticle genes expression 48 . The overexpression of so many CP genes in this study indicated that this gene family may be involved in indoxacarb resistance in S. litura.
A total of 737 genes were down-regulated (44.7% of total co-DEGs) in the resistant strains. This phenomenon is very understandable since all living organisms have limited energy inputs. Reduced expressions in these genes would save some energy for resistant S. litura since the other 55.3% of total co-DEGs were up-regulated. It is generally recognized that overexpressed genes may play more important roles than down-regulated genes in the insecticide resistance. The decreased expression of some metabolic detoxification genes might result from the responses to various endogenous and exogenous compounds, or path of physiological signals 49 . Many studies have deduced the balance mechanism of up-and down-regulated genes, including an adaptive homeostatic response that protects the cell from the deleterious effects of oxidizing species, such as nitric oxide and arachidonic acid metabolites from catalytic and/or metabolic enzymes 50 . In return, this also could be a pathological response to inflammatory processes and a need for the tissue to utilize its transcriptional machinery and energy for the synthesis of other components involved in insecticide resistance 51 .
Analysis of mRNA through high throughput expression profiling using transcriptome analysis methods has provided considerable advances in understanding the molecular base of resistance in insects 52 . In M. domestica, 1316 genes were identified as being up-regulated in the multiple insecticide resistant strain in comparison to the susceptible strains by transcriptome analysis, and the majority of these up-regulated genes fell within the structural classification of proteins (SCOP) categories of metabolism, regulation and intracellular processes 53 . A similar whole transcriptome study has been carried out in Culex quinquefasciatus. The results showed that 367 genes were found to be up-regulated in the highly permethrin-resistant strain, and all P450 genes were up-regulated by at least twice 54 . Bai et al. 55 also revealed that P450 genes (CYP9) were highly expressed in pesticide-exposed Cimex lectularius populations by transcriptome and qPCR analysis. In this study, 910 co-up-regulated and 737 co-down-regulated genes were identified in the indoxacarb-resistant strains through mRNA sequencing, and the roles of 126 detoxification-associated genes (107 co-up-regulated and 19 co-down-regulated) in indoxacarb resistance were systematically analyzed. These results provide clues to the identification of potential detoxification genes involved in indoxacarb resistance in S. litura. Although the two resistant strains have different genetic background, it can be more reasonably screened to obtain insecticide resistance-related genes that reflect the real situation in the field. A similar approach has been used to study the molecular mechanisms of spirodiclofen resistance in T. urticae 56 .
To date, more than 30,000 miRNAs have been found in over 100 organisms 57 . In C. pipiens, 100 known miR-NAs and 42 novel miRNAs were identified, and 28 miRNAs were differentially expressed in the susceptible and deltamethrin-resistant strains 14 . Seventy-five known miRNAs and 64 novel miRNAs were also identified in the susceptible and fenpropathrin-resistant strains of T. cinnabarinus, including 12 differentially expressed miRNAs 17 . In this study, we identified 91 known miRNAs and 52 novel miRNAs in S. litura, 10 of which were co-differentially expressed in the InRS and FInRS. This data has detected significantly more miRNAs than the early report by Rao et al. 13 , in which only 58 miRNAs were identified among different developmental stages in S. litura. The reasons that we can detect increased numbers of miRNAs in our study compared to Rao et al. 13 were the availability of S. litura whole genome, the improved methods and the upgraded miRBase database.
In general, miRNAs play important gene-regulatory roles by targeting the mRNAs of protein coding genes and repressing their post-transcriptional properties. In this way, down-regulation of a miRNA indicates increased activity of its target gene. It is an important step to identify the target genes of miRNAs for understanding their roles in gene regulatory networks. Therefore, we analyzed the relationship between indoxacarb resistance-related co-DEGs and co-DEmiRNAs in S. litura. Among 126 indoxacarb resistance-related co-DEGs, 33 up-regulated co-DEGs were predicted target genes of 5 down-regulated co-DEmiRNAs, 4 down-regulated co-DEGs were predicted target genes of 5 up-regulated co-DEmiRNAs. With the deepening of the study of miRNA function, some miRNAs have been shown to be associated with insecticide resistance by regulating resistance related genes in insects and mites. Hong et al. 14 and Lei et al. 15 indicated that miR-71 and miR-278-3p were involved in pyrethroid resistance by targeting CYP325BG3 and CYP6AG11 in C. pipiens and C. pipiens pallens, respectively. Zhang et al. 17 demonstrated that Tci-miR-1-3p could regulate the cyflumetofen resistance through TCGSTM4 in T. cinnabarinus. Ma et al. 58 also indicated that miR-92a regulates pyrethroid resistance through its interaction with a cuticular protein gene, CpCPR4 in C. pipiens pallens. Thus, our results indicated that co-DEmiRNAs of S. litura may function in the formation of the indoxacarb resistance through regulating the insecticide resistance-related co-DEGs.
In www.nature.com/scientificreports www.nature.com/scientificreports/ of S. litura. Our data revealed that indoxacarb selection strongly affected the transcription levels of 126 detoxification-associated genes (P450, GST, CCE, ABC transporter, UGT, APN, sialin, SP and CP) involved in indoxacarb resistance. Further analysis showed that 37 indoxacarb resistance-related co-DEGs were predicted to be the target genes of 10 co-DEmiRNAs, and indicated that these miRNAs may regulate the indoxacarb resistance through these indoxacarb resistance-related genes in S. litura. Overall, the present study help us better understand the indoxacarb resistance mechanisms and the regulatory mechanisms of miRNAs in S. litura at the molecular level. Functional analysis of these resistance related genes and verification of the regulatory relationships between these miRNAs and their target genes are needed to further elucidate the indoxacarb resistance mechanisms in S. litura.

Methods
Ethics statement. The laboratory population of S. litura was obtained from the Institute of Zoology, Beijing, China. The field population of S. litura was collected from the field crops of Changsha City, Hunan Province, China. There was no specific permission required for these collection activities because this insect is a kind of agriculture-harmful pest and distributes worldwide. We confirm that the field collection did not involve endangered or protected species. Bioassays were conducted with fourth-instar larvae of S. litura using the artificial diet dipping method 59 . Briefly, indoxacarb (15%, E.I. DuPont de Nemours and Co., Inc., Wilmington, DE, USA) was dissolved in sterilized water to at least 5 concentrations and the mortality was kept at 20-80%. The artificial diet was cut into an area of 2 cm 2 and a thickness of 5 mm, and dipped into the indoxacarb solution for 10 s, including sterilized water as control. These artificial diets were air dried at room temperature for 5-10 min. Then, fifteen larvae were placed on each treated artificial diet. Each dose was performed in three replicates. The bioassays were kept at 25 ± 2 °C, 65 ± 5% RH and 14:10(L:D) photoperiod. Mortality was recorded 24 h after exposure for indoxacarb. Larvae were considered dead if they failed to make a coordinated movement when prodded with a brush. LC 50 values were calculated via probit analysis using PoloPlus software (LeOra Software Inc., Berkeley, CA, USA).

RNA isolation.
Total RNAs were extracted from fourth-instar larvae of SS, InRS and FInRS using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), respectively. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity, concentration and integrity were measured using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA), Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively. mRNA sequencing and data analysis. Qualified total RNAs from three strains (three biological replicates for each strain) were used for preparing sequencing libraries using NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's recommendations. The library preparations were sequenced on an Illumina Hiseq 2500 platform and 125 bp paired-end reads were generated at the Novogene Bioinformatics Institute (Beijing, China) after the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. In order to obtain clean reads, the low-quality reads, adaptor sequence, ploy A or T or G or C and duplication sequence were removed from raw data. Then, clean reads were aligned to the reference genome of S. litura 10 using Hisat2 v2.0.5 60 .
Gene expression levels were estimated by FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced) 61 . DEGs analysis was performed using the DESeq2 R package (1.16.1) 62 . Genes with an adjusted P < 0.05 and |log 2 (fold change)| ≥ 1 were assigned as differentially expressed. Gene Ontology (GO) database 63 was used to identify functional modules of DEGs (P < 0.05). sRNA sequencing and data analysis. Qualified total RNAs from three strains (three biological replicates for each strain) were used as input materials for the sRNA libraries. Sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, USA.) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, total RNAs were ligated to 5′ and 3′ adaptors, then first strand cDNA was synthesized by reverse transcriptase. After PCR amplification of the cDNA, the amplified PCR products within 140-160 bp were separated and purified by a 8% polyacrylamide gel. Libraries quality were assessed on the Agilent Bioanalyzer 2100 system using DNA High Sensitivity Chips. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and 50 bp single-end reads were generated at the Novogene Bioinformatics Institute (Beijing, China). Clean reads were obtained by removing reads containing 'N' (an unrecognized base) at 10% or higher, with 5′ adapter contaminants, without 3′ adapter or the insert tag, containing ploy A or T or G or C and low-quality reads. After removing reads with a sequence shorter than 18 nt or more than 35 nt, clean reads were mapped to the reference genome of S. litura using Bowtie software 64 . Then, the matched sRNAs were compared with the mature miRNAs in miRBase 20.0 (http://www.mirbase.org/) to looking for known miRNA of S. litura, only perfectly matches were accepted and counted. Next, the remained sRNAs were compared with Repeatmasker and Rfam database (ftp://selab.janelia.org/pub/Rfam) to remove protein-coding genes, repeat sequences, ribosomal RNA (rRNA), small cytoplasmic RNA (scRNA), transfer RNA (tRNA), small nuclear RNA (snRNA) and small nucleolar RNA