Transcriptome analysis of molecular mechanism to pyrethroids resistance in Spodoptera litura

Spodoptera litura is a destructive agricultural pest and had evolved resistance to multiple insecticides, especially pyrethroids. At present, the resistance mechanism to pyrethroids remains unclear. Four field-collected populations, namely CZ, LF, NJ and JD, were identified to have high resistance to pyrethroids, and the resistant ratio ranged from 40.8 to 1764.0-fold. To explore pyrethroids resistance mechanism, the transcriptome between pyrethroid-resistant (LF and NJ) and pyrethroid-susceptible (GX) populations were compared by RNA-sequence. Results showed that multiple differential expressed genes were enriched in metabolism-related GO terms and KEGG pathways. 35 up-regulated metabolism-related genes were screened to verify by qRT-PCR. Consistent up-regulation of 13 unigenes, including 3 P450s, 4 GSTs, 1 UGTs, 4 COEs and 1ABC, were verified in the additional pyrethroids resistant populations CZ and JD. The expression level of CYP3 and GST3, which were homologous to CYP321A8 and GST1, respectively, showed good correlation with their pyrethroids resistance level among CZ, LF, NJ and JD populations. While the expression level of CYP12, CYP14, COE4 and ABC5 showed good correlation with their pyrethroids resistance level in at least three populations. UGT5 had the highest expression level among the tested UGT genes in the four pyrethroids resistant populations.


Background
The common cutworm, Spodoptera litura (Fabricius), is a polyphagous agricultural pest worldwide and caused enormous losses to many economical crops like cotton, soybean, tobacco and vegetables. Insecticides resistance developed rapidly in the last decades due to its extensively application. To date, S. litura has evolved high resistance to multiple conventional and new insecticides, which led to both failure control of pest and the destruction of natural enemies [1][2][3][4][5][6].
Insecticides resistance is the consequence of pest evolutionary adaptation to insecticides selection, a genetic phenomenon, which is conferred by target-site-based resistance (TSR) and non-target-based resistance (NTSR) mechanisms [7]. TSR is endowed by gene mutations in target enzymes, which can be easily detected by sequencing. Multiple mutations in sodium channel conferring resistance to pyrethroid insecticides had been identified in insect species, such as L1014F, L1014S, M918V, L925I, E435K, C785R, P1999L, F1538I, A1410V, M1524I, D1549V, T929I and M918T [8][9]. NTSR is achieved by mechanisms reducing insecticide concentration reaching the target-site, including metabolism-based resistance, which enhanced detoxification rate of insecticide by cytochrome P450 monooxygenase (P450), glutathione S-transferase (GST), carboxylesterases (COE), UDP-glycosyltransferase (UGT) and ATP-binding cassette transporters (ABC) [10]. Synergism is the primary step to identify the involvement of metabolism-related enzyme by piperonyl butoxide (PBO), diethyl maleate (DEM) and S,S,Stributylphosphorothioate (DEF), which are inhibitors of P450, GST and COE, respectively [7]. The determination of putative gene expression and detoxification enzyme activity is also employed to explore the mechanism of metabolism-based resistance. However, these methods only establishing the resistance phenotype is metabolic involved, fails to comprehensively identify the individual genes responsible for insecticide resistance.
The transcriptome is the complete set of all transcripts within a cell or tissue, that reflect the mRNA expression levels of different genes in a given sample [11]. RNA-Sequence (RNA-Seq) provided opportunity for the comprehensive analysis of pest transcriptome without reference genome in greater depth than ever before [12]. Although multiple researches have investigated the resistance of S. liturato various insecticides, molecular mechanism of metabolic resistance in S. litura is still not well defined due to lack of detailed genetic and sequence information [13][14]. In addition, pyrethroid insecticides were widely applied to control S. litura and other agricultural pests, for its high efficiency and low toxicity to mammalian animals. Yet excessive usage has imposed pyrethroids resistance the most severe situation [15]. To minimize the evolution of resistance and prolong the life time of insecticide, the molecular mechanisms for pyrethroids resistance in this pest should be declared. In this study, RNA-seq was performed by Illumina HiSeq™ 2000 sequencing platform and the transcriptome of pyrethroid-resistant and pyrethroidsusceptible populations were compared to provide new insights into mechanism of metabolism-based resistance to pyrethroids in S. litura.

Toxicity of insecticides to S. litura
Eleven insecticides were chosen to test the toxicity to the larvae of GX, LF and NJ populations. High level resistance to fenvalerate, beta-cypermethrin and cyhalothrin were found, which were 809.5-, 141.3-, 1764.0-fold in LF population and 322.0-, 40.8-, 951.0fold in NJ population, respectively, compared with that in GX population (Table 1-3). LF and NJ populations showed no resistance to phoxim, profenofos, chlorpyrifos and emamectin benzoate, with resistance ratio (RR) lower than 5.0-fold. However, low resistance had been developed to chlorantraniliprole, bromine cyantraniliprole, imidacloprid and methomyl, with RR ranged from 4.0 to 19.8-fold. The results above indicated that field-collect populations of S. litura had high level resistance to pyrethroid, which was urgent to be managed. (Table 1-3 Insert here) To determine the involvement of detoxification enzyme for pyrethroids resistance in S. litura, the synergists of P450, COE and GST, namely PBO, DEF and DEM, were used before pyrethroids treatment. Results showed that PBO increased pyrethroids toxicity distinctly in both LF and NJ populations, while DEF and DEM increased fenvalerate and cyhalothrin toxicity slightly and had no influence on beta-cypermethrin toxicity ( Illumina sequencing and reads assembly RNA-Seq was conducted with three cDNA libraries of GX, LF and NJ populations, respectively. Clean reads ranging from 51429712 to 60901328 per cDNA library were obtained by removing reads containing adapter, ploy-N and low quality reads from raw data ( Table 4). According to the sequence analysis, almost 97% raw reads had Phred-like quality score at the Q20 level (an error probability of 1%), with GC content around 47%.
Clean reads were used to assemble 140003 transcripts and 82014 unigenes ranging from 201 to 29587 bp (Table 5). Pearson correlation coefficients of different RNA samples in the same population ranged from 0.868 to 0.916, which was higher than 0.525 to 0.648 between different populations (Additional file 1).
( ( Table 6 Insert here) GO assignments were used to predict gene function. Total 82014 unigenes were classified into 54 subgroups, including molecular function (MF), biological process (BP) and cellular component (CC). In BP, cellular process was the largest subgroup, followed by metabolic process and single-organism process. In CC, the largest two subgroups were cell and cell part, respectively. For MF, binding and catalytic activity formed the largest two subgroups  The expression of candidate DEGs in additional two field populations CZ and JD are another two field populations with high resistance level to pyrethroid insecticides. Compared with GX population, CZ showed 84.5-fold resistance to fenvalerate, 11.5-fold resistance to beta cypermethrin and 678.5-fold resistance to cyhalothrin (Table   8), while JD showed 9123.5-fold resistance to fenvalerate, 642.0-fold resistance to beta cypermethrin and 6986.0-fold resistance to cyhalothrin (Table 9).  (Table 10). (

Disscusion
Insecticides continue to be key weapon for insect control in agriculture. However, the evolution of insecticides resistance due to strong selection imposed by widespread insecticides usage threatens the weapons' availability. Studying the molecular mechanism of insecticides resistance, especially for metabolism-based resistance, provides valuable insights into how genetic, physiological and biochemical changes in insects and opens windows for improving future pest control. In the present study, the transcriptome of GX, NJ and LF populations were sequenced by Illumina platform with high criteria of data quality (Table 4), enabling assembling of 140003 transcripts and 82014 unigenes (Table 5-6). GO and KEGG enrichment analysis of DEGs indicated that LF and NJ populations have relative high ability in oxidation-reduction activity, hydrolase activity and P450 involved metabolic process than GX population (Figure 4, 6). Recently, Cheng et al. [16] conducted the genome sequencing of S. litura, and 138 P450, 47 GST, 110 COE and 54 ABC genes were annotated, which was consistent with our transcriptome results (Table 11). In addition, the validation of genes expression level showed that RNA-Seq results were well coincident with that of qRT-PCR (Table 10). Therefore, the RNA-Seq results were reliable and this information will facilitate further researches on pyrethroids resistance and other physiological processed in S. litura. (Table 11 Insert here) Bioassay result in LF and NJ populations of S. litura showed high level resistance to fenvalerate, beta cypermethrin and cyhalothrin (Table 1-3). However, resistance to beta cypermethrin was much lower than that to fenvalerate and cyhalothrin in both LF and NJ populations, which might result from different sets of detoxification enzymes to each insecticide. No cross-resistance between organophosphates and pyrethroids was found in LF or NJ populations, which is similar with the result reported by Huang et al [13]. While Tong et al.found high level resistance to both organophosphates and pyrethroids in five field S. litura populations during 2010 to 2012 [14]. This opposite phenomenon could be explained by metabolic-based activation and passivation of insecticides. Multiple studies had revealed a direct link between elevated P450 activities and organophosphate and pyrethroids resistances in insects [17]. However, some P450s can also be involved in the activation of insecticides, chemically modifying insecticides to biologically active forms, as is the case for some organophosphorus insecticides [7]. For example, chlorpyrifos can be converted to chlorpyrifos-oxon by P450 and chlorpyrifos-oxon is 100-fold toxicity than chlorpyrifos [18][19]. Thus different sets of up-regulated P450s in insects could cause cross resistance or negative cross resistance between organophosphates and pyrethroids.
Cytochrome P450 monooxygenases can mediate resistance to most insecticides for their genetic diversity, broad substrate specificity and catalytic versatility. Go enrichment of up-regulated DEGs in LFvsGX and NJvsGX groups showed that oxidation-reduction process and metabolic process were significantly enriched ( Figure 4). Furthermore, P450s involved metabolism of xenobiotics and drug metabolism pathways were also enriched in KEGG enrichment analysis ( Figure 6). Besides, P450s inhibitor (PBO) significantly increased the toxicity to pyrethroids in LF and NJ populations ( Figure 1). Therefore, up-regulated P450s should be the major reason for pyrethroids resistance in LF and NJ populations. The qRT-PCR results indicated that most selected P450 genes were significantly over expressed in the four resistant populations, especially for CYP2, CYP3, CYP10, CYP12 and CYP14 (Table   10). Furthermore, the expression level of CYP3, CYP12 and CYP14, which were annotated as CYP321A8, CYP6AE43 and CYP6AB31, respectively, showed good correlation with pyrethroids resistance level among GX, CZ, NJ, LF and JD populations. CYP321A8 and CYP6AB31 could be induced by deltamethrin, lambda-cyhalothrin and lambda-cyhalothrin, respectively, and had high expression level in midgut and fat body [20][21][22]. CYP6AE43 of S. frugiperdacould only be detected in the midgut and malpighian tubules [23].However, there is building evidence that the key tissues for the metabolism of most compounds are midgut, malpighian tubules and fat body [23][24]. Therefore, P450 genes that had higher expression levels in these tissues may be involved in xenobiotic detoxification. Hence, most up-regulated P450s might be the set of pyrethroids detoxified genes in S. litura, especially for CYP3, CYP12 and CYP14 in consideration of the higher expression level and detoxify function.
COE were classified as Phase I detoxification enzymes and perform hydrolysis reaction at the first step of detoxification. It has been demonstrated that COE were involved in insecticide resistance, especially for organophosphates, pyrethroids and carbamates [25].
For instance, the over-expression of carboxylesterase E4 in Myaus persicae enhanced degradation and sequestration of pyrethroids, organophosphates and carbamates [26].
According to Go enrichment analysis of up-regulated DEGs, hydrolase activity was enriched in LF and NJ populations (Figure 4). Validation in qRT-PCR revealed that the 4 selected COEs were higher expressed in resistance populations, especially for COE4 (Table  10). Although only slight synergism to fenvalerate and cyhalothrim by DEF was observed ( Figure 1), COE may also play an important role in pyrethroids resistance in consideration of the enrichment analysis and high expression level in the four resistance populations.
The products of Phase I reactions often become substrates for Phase II enzymes, such as GST and UGT, that add glutathione and glycosyl groups, respectively. Insect GSTs are a superfamily of detoxification enzymes and had been reported to be associated with resistance to pyrethroids, organophosphorus and organochlorine insecticides [27].In contrast to organochlorine and organophosphorus resistance, GSTs are not involved in the direct metabolism of pyrethroids but rather are involved in the sequestration of pyrethroids and/or the detoxification of lipid peroxidation products induced by pyrethroids [28][29]. Functional study of GSTs in Nilaparvata lugens showed that NlGSTD1 had peroxidase activity and confers pyrethroids resistance via gene amplification [29]. Two S.  [39]. Our study provides a comprehensive understanding that CYPs, COEs, GSTs, UGTs, and ABC transporters may associated with pyrethroids resistance in S. litura, which will greatly extend our understanding of metabolism-based insecticide resistance mechanism, and facilitate the further identification of pyrethroids resistance involved genes. Quality control, Transcriptome assembly and Gene functional annotation Clean data were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. Q20, Q30, GC-content and sequence duplication level of the clean data were calculated and all the downstream analyses were conducted based on clean data with high quality. Transcriptome assembly was accomplished using Trinity [40] with min_kmer_cov set to 2 by default and all other parameters set default. Gene function was annotated based on Nr, Nt, Pfam, Swiss-Prot database and KOG/COG, with a significant threshold of E-value ≤ 10 -5 . The GO terms for functional categorization were analyzed using Blast2go software with the E-value threshold ≤ 10 -6 . The pathway assignments were carried out by sequence searches against the KEGG database, with the E-value threshold ≤ 10 -10 .

Differential expression analysis of unigenes
Differential expression analysis between the three populations was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p-value<0.05 found by DESeq were assigned as differentially expressed.

GO enrichment analysis of the DEGs was implemented by the GOseq R packages based on
Wallenius non-central hyper-geometric distribution [41]. KOBAS software was used to test the statistical enrichment of DEGs in KEGG pathways [42].

Validation of the RNA-Seq DEGs
To validate the transcriptional results from RNA-Seq and computational analysis, 35 unigenes were selected for quantitative RT-PCR (qRT-PCR). EF1α and RPL10 were used as housekeeping genes to normalize the expression of each gene [43].First strand cDNA was  Resistance ratio (RR) = LD 50 value of field population/ LD 50 value of GX population