Toxicities of insecticides to S. litura
Eleven commonly used insecticides were applied to the larvae of GX, LF and NJ populations to test their toxicities. When compared with the insecticide-susceptible population, named as GX, higher resistant level to fenvalerate, beta-cypermethrin and cyhalothrin were observed, which were 809.5-, 141.3-, 1764.0-fold in LF population and 322.0-, 40.8-, 951.0-fold in NJ population, respectively (Table 1). 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, cyantraniliprole, imidacloprid and methomyl, with RR ranging from 4.0 to 19.8-fold. The results above indicated that field-collect populations of S. litura had high levels of pyrethroid resistance, which was urgent to be managed.
To determine the involvement of detoxification enzymes for pyrethroid resistance in S. litura, PBO, DEF and DEM, which were the synergists of P450, COE and GST, respectively, were used before pyrethroids treatment. Results showed that pyrethroid toxicities were distinctly increased in both LF and NJ populations by PBO. Fenvalerate and cyhalothrin toxicities were slightly increased by DEF and DEM, while beta-cypermethrin toxicity was not influenced by DEF and DEM (Figure 1).
Illumina sequencing and reads assembly
RNA-Seq was conducted with the cDNA libraries of GX, LF and NJ populations. 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 (Additional file 1). More than 96.5% raw reads had Phred-like quality score at Q20 level (an error probability of 1%). Clean reads were used to assemble 140003 transcripts and 82014 unigenes ranging from 201 to 29587 bp (Additional file 2). Pearson correlation coefficients of RNA samples in the same population ranged from 0.868 to 0.916, which was higher than that in different populations (Additional file 3).
Functional annotation of unigenes
Seven databases, NCBI non-redundant protein sequences database (Nr), NCBI non-redundant nucleotide sequences data-base (Nt), Protein family (Pfam), Clusters of Orthologous Groups of proteins (KOG/COG), Swiss-port, KyotoEncyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) database were used for unigenes functional annotation. There are 47185 (57.53%) unigenes annotated in at least one database, and 8974 (10.94%) unigenes annotated in all databases. Among the annotated unigenes, 50.24% hit in Nr database, 39.15% in GO database and 38.87% in Pfam database (Additional file 4). More than 90% annotated unigenes had nucleotides similarity above 60% (Additional file 5). The results of Nr database annotation indicated that S. litura has 41.3% sequence similarity to Bombyx mori, 16.2% to Danaus plexippus, 14.1% to P. xylostella, 1.9 % to Helicoverpa armigera and 1.6% to Papilio xuthus (Additional file 6).
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. For MF, binding and catalytic activity formed the largest two subgroups (Additional file 7).
KEGG database was used for understanding high-level functions and utilities of the biological system from molecular-level information. 20164 unigenes were divided into five groups, including cellular processes (A), environmental information processing (B), genetic information processing (C), metabolism (D), and organism systems (E). 3002 unigenes (14.89%) were involved in A; 3023 unigenes (14.99%) were involved in B; 2954 unigenes (14.65%) were involved in C; 5563 unigenes (27.59%) were involved in D; 5622 unigenes (27.88%) were involved in E (Additional file 8).
Differential gene expression among LFvsGX and NJvsGX groups
In order to explore the molecular mechanism of pyrethroid resistance, differentially expressed genes (DEGs) were selected according to q-value<0.05 and |log2 (Fold change)|>1. Compared with the transcriptome of GX population, there are 18043 DEGs (10063 up-regulated and 7980 down-regulated) in LF population and 20370 DEGs (11147 up-regulated and 9223 down-regulated) in NJ population, respectively. And 11162 DEGs were found in both LFvsGX and NJvsGX groups, including 5720 up-regulated and 5442 down-regulated genes (Additional file 9).
GO and KEGG enrichment of the DEGs
To further analyze the function of DEGs, GO and KEGG enrichment analyses were conducted. Based on GO enrichment, the up-regulated unigenes in LFvsGX group were enriched into multiple GO terms belonging to 20 BP, 5 CC and 20 MF; while 20 BP and 20 MF in NJvsGX group. In BP, the top three largest GO terms were metabolic process, single-organism process and oxidation-reduction process in both LFvsGX and NJvsGX groups. In MF, the top three largest GO terms were catalytic activity, hydrolase activity and oxidoreductase activity in LFvsGX group, while catalytic activity, oxidoreductase activity and peptidase activity in NJvsGX group (Figure 2A, B).
Much less unigenes were significantly enriched in the down-regulated GO terms in both LFvsGX and NJvsGX groups, including chitin metabolic process and chitin binding, which indicated that different epidermal penetration rate may be involved in pyrethroid resistant S. litura populations (Figure 2C, D).
Among the top twenty enriched KEGG pathway terms in LFvsGX and NJvsGX groups, fourteen terms were shared in both groups, including metabolism of xenobiotics by cytochrome P450 and drug metabolism-cytochrome P450 pathways (Figure 3). The results of GO and KEGG enrichment suggested that metabolism was the main enriched terms of DEGs and should be a main reason for pyrethroid resistance in S. litura.
Validation of the candidate DEGs
Based on high reading count and significant high expression levels in RNA-Seq, 14 P450, 5 COE, 5 GST, 5 UGT and 6 ABC up-regulated unigenes were screened out and validated by quantitative RT-PCR (qRT-PCR). Among the 35 selected DEGs, 30 unigenes (including 14 P450, 5 GST, 5 UGT, 4 COE, and 1 ABC) were both significantly up-regulated by qRT-PCR validation in LFvsGX and NJvsGX groups (Table 2). CYP12 (annotated as CYP6AE43) had the highest expression level among all the P450 unigenes, with 358.27-fold in LF and 237.95-fold in NJ population. UGT5 (annotated as UGT33T2) had the highest expression level among all the UGT unigenes, with 309.49-fold in LF and 3236.54-fold in NJ population.
The expression of candidate DEGs in additional two field populations
CZ and JD are another two field collected 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, while JD showed 9123.5-fold resistance to fenvalerate, 642.0-fold resistance to beta cypermethrin and 6986.0-fold resistance to cyhalothrin (Table 3).
Based on the pyrethroid RR, resistance level of JD>LF>NJ>CZ>GX could be concluded. To further validate the role of candidate DEGs in pyrethroid resistance, the expression level of the up-regulated DEGs in LF and NJ populations were also determined by qRT-PCR in CZ and JD populations. Results indicated that 13 unigenes (CYP3, CYP12, CYP14, GST1, GST2, GST3, GST4, UGT5, COE1, COE2, COE4, COE5 and ABC5) were also significantly higher expressed in the two additional populations. UGT5 had the highest expression level among all the tested UGT genes in the four pyrethroid resistant populations. The expression level of CYP3 and GST3 among the four field-collected populations conformed well with their pyrethroid resistance level. While the expression level of CYP12, CYP14, COE4 and ABC5, showed good correlation with their pyrethroid resistant level in at least three field-collected populations (Table 4).