Identification of transcriptome and fluralaner responsive genes in the common cutworm Spodoptera litura Fabricius, based on RNA-seq

Background Fluralaner is a novel isoxazoline insecticide with a unique action site on the γ-aminobutyric acid receptor (GABAR), shows excellent activity on agricultural pests including the common cutworm Spodoptera litura, and significantly influences the development and fecundity of S. litura at either lethal or sublethal doses. Herein, Illumina HiSeq Xten (IHX) platform was used to explore the transcriptome of S. litura and to identify genes responding to fluralaner exposure. Results A total of 16,572 genes, including 451 newly identified genes, were observed in the S. litura transcriptome and annotated according to the COG, GO, KEGG and NR databases. These genes included 156 detoxification enzyme genes [107 cytochrome P450 enzymes (P450s), 30 glutathione S-transferases (GSTs) and 19 carboxylesterases (CarEs)] and 24 insecticide-targeted genes [5 ionotropic GABARs, 1 glutamate-gated chloride channel (GluCl), 2 voltage-gated sodium channels (VGSCs), 13 nicotinic acetylcholine receptors (nAChRs), 2 acetylcholinesterases (AChEs) and 1 ryanodine receptor (RyR)]. There were 3275 and 2491 differentially expressed genes (DEGs) in S. litura treated with LC30 or LC50 concentrations of fluralaner, respectively. Among the DEGs, 20 related to detoxification [16 P450s, 1 GST and 3 CarEs] and 5 were growth-related genes (1 chitin and 4 juvenile hormone synthesis genes). For 26 randomly selected DEGs, real-time quantitative PCR (RT-qPCR) results showed that the relative expression levels of genes encoding several P450s, GSTs, heat shock protein (HSP) 68, vacuolar protein sorting-associated protein 13 (VPSAP13), sodium-coupled monocarboxylate transporter 1 (SCMT1), pupal cuticle protein (PCP), protein takeout (PT) and low density lipoprotein receptor adapter protein 1-B (LDLRAP1-B) were significantly up-regulated. Conversely, genes encoding esterase, sulfotransferase 1C4, proton-coupled folate transporter, chitinase 10, gelsolin-related protein of 125 kDa (GRP), fibroin heavy chain (FHC), fatty acid synthase and some P450s were significantly down-regulated in response to fluralaner. Conclusions The transcriptome in this study provides more effective resources for the further study of S. litura whilst the DEGs identified sheds further light on the molecular response to fluralaner.


Background
The common cutworm, Spodoptera litura Fabricius (Lepidoptera: Noctuidae), is a destructive polyphagous pest with a broad host plant range of more than 100 species of crops and vegetables [1]. It is widely distributed around the world and has evolved high resistance to most conventional insecticides, including carbamates, pyrethroids and organophosphates [2,3]. As a result, farmers are using more insecticides leading to serious environmental pollution [4]. Therefore, exploring alternative and environmentallyfriendly insecticides is one important approach to reduce the economic loss from S. litura.
The transcriptome has been widely used as a powerful tool for exploring the molecular mechanisms of genes that respond to insecticides and toxins. For instance, Cui et al. (2017) explored the responsive genes of the Asian corn borer (Ostrinia furnacalis Guenée) to flubendiamide [10]. Song et al. (2017) found that immunity-, metabolism-and Bt-related genes in S. litura midgut were responsive to Vip3Aa toxin, highlighting that trypsin was possibly involved in Vip3Aa activation [11]. Xu et al. (2014) found that in the carmine spider mite, T. cinnabarinus (Boisduval), 10 gene categories were putatively involved in insecticide resistance [12]. Marco et al. (2017) found that the defensome family genes of the mosquito, Anopheles stephensi Liston, responded to toxicants during insecticide exposure [13].
In the present study, the transcriptomes of control and fluralaner-treated S. litura were sequenced with the Illumina HiSeq Xten (IHX) platform (Illumina, Inc., San Diego, CA) and assembled according to the reference genome [14]. Hundreds of new genes were annotated by the NR, GO, COG and KEGG databases. The differentially expressed genes (DEGs) between control and fluralaner-treated S. litura were identified. In particular, detoxification enzyme genes and insecticide-targeted receptor genes were investigated, and fluralaner responsive genes relating to detoxification and development were also analyzed.

Results
RNA-seq, sequence assembly, transcript analysis and functional classification The S. litura transcriptomic data were generated from 9 different libraries on the IHX platform. After removing the reads with adaptor and low quality, 70.97 Gb clean reads were produced with Q 30 > 94.15%. After alignment, clean reads (85.71-88.22%) were successfully mapped to the S. litura genome, of which 85.83% mapped to exons and the others mapped to intergenic regions or introns. Finally, 16,926 genes were assembled in the S. litura transcriptome. Subsequently, 16,572 genes, including 451 new genes, were successfully annotated according to COG (6142), GO (7887), KEGG (6812), and NR (16,554) databases. In the transcriptome, 12,707 annotated genes were > 1000 bp long whilst 3850 annotated genes were between 300 and 999 bp.
Through comparing the S. litura transcripts against NR databases, most annotated genes (≥ 99.26%) were highly homologous to Lepidopteran genes. To be specific, 15,884 (95.95%) genes annotated by NR databases were homologous with those of S. litura, followed by 238 (1.44%) and 147 (0.89%) genes that were homologous to those of H. armigera Hübner and Heliothis virescens (Fabricius), respectively. Moreover, homologous genes identified by significant Blast hits from nine insect species are shown in Fig. 1a. According to the COG database, 6142 genes were annotated and classified into 25 specific categories and each gene was classified into at least one category ( Fig. 1b and Additional file 1). According to the GO database, 7887 annotated genes were primarily divided into three categories of "cellular component", "molecular function" and "biological process", and were further classified into 51 functional sub-categories, including 21, 16 and 14 sub-categories in the "biological process", "cellular component" and "molecular function" categories, respectively ( Fig. 2 and Additional file 1). According to the KEGG database, 6812 annotated genes were involved in 154 pathways with "purine metabolism" (174 genes, 2.55%), "carbon metabolism" (154 genes, 2.26%) and "peroxisome" (143 genes, 2.10%) being the largest three pathways (Additional file 2).

Expression profile analysis of fluralaner responsive genes
In general, the transcriptome can represent the response of gene expression upon exposure to chemicals [10,15]. To determine the influence of fluralaner on the expression profile of genes, the change of genes was analyzed using the Fragments Per Kilobase of transcript per Million fragments mapped (FPKM) method [16]. A total of 3725 DEGs were found in LC 30 and LC 50 fluralaner-treated groups. The Venn diagram shows the number of shared and exclusive DEGs at each group, and a total of 3275 (1695 up-regulated and 1580 down-regulated) or 2491 (1500 up-regulated and 991 down-regulated) DEGs were found after exposure to LC 30 and LC 50 fluralaner compared with the control group, respectively. 100 DEGs were shared between the LC 30 and LC 50 treated groups (Fig. 3). 1209 and 406 DEGs were exclusive to LC 30 and LC 50 exposure to fluralaner, respectively, whilst 2041 genes DEGs were shared between all treated groups (Fig. 3).   Subsequently, the DEGs in LC 30 and LC 50 fluralanertreated groups were analyzed by GO, COG and KEGG databases. 1842 and 1297 DEGs annotated by GO databases were divided into 50 and 48 sub-categories, respectively. The major sub-categories were "membrane" with 592 (32.14%) and 455 (35.08%) DEGs, respectively, "catalytic activity" with 823 (44.68%) and 580 (44.72%) DEGs, and "metabolic process" with 814 (44.19%) and 531 (40.94%) DEGs (Fig. 2). 1238 and 880 DEGs were annotated by the COG database and divided into 24 categories. The "general function prediction only" was the largest category with 187 (15.11%) and 131 (14.89%) DEGs, followed by 151 (12.20%) and 128 (14.55%) DEGs in the "posttranslational modification, protein turnover, chaperones" category (Additional file 12). 1728 and 1207 DEGs were annotated by the KEGG database and classified into 128 and 124 pathways with 62 and 41 DEGs being in the pathway "purine metabolism" in LC 30 and LC 50 fluralanertreated groups, respectively (Additional file 13).

Discussion
In the present study, the transcriptome of S. litura untreated and treated with fluralaner was assembled based on the S. litura genomic database resulting in 16,572 genes being annotated according to NR, COG, GO, and KEGG databases. The number of annotated genes were almost consistent to those (16,161 genes) in a de novo transcriptome of S .litura [17], whereas more than those (11,692 genes) based on sequence homology of S. litura [18] and the reference genome (15,317 genes) [14].
In insects, insecticides often induce the increase of metabolic activity of P450s, GSTs and CarEs [19,20]. P450s are members of an important metabolic enzyme superfamily, which are involved in the metabolism of xenobiotic compounds such as insecticides and plant secondary metabolites [21]. Generally, there are four large clades of insect P450 genes that are CYP2, CYP3, CYP4 and CYPM [22]. In the present study, less genes of CYP3 and CYP4 clades and almost equal number of genes of CYP2 and CYPM clades were identified when compared to those from the reference genome (61, 58, 11 and 8 genes of CYP3, CYP4, CYPM and CYP2 clades, respectively) [14]. The CYP3 and CYP4 are the largest clades and are mainly responsible for xenobiotic metabolism and insecticides resistance [22]. In this study, 90 genes accounting for 84.11% of the total P450 genes identified were clustered into the CYP3 and CYP4 clades, which were consistent with those containing 86.23% P450 genes in the reference genome [14]. In this study, 16 P450 genes belonging to the CYP4, CYP6, CYP9 and CYP12 sub-families were overexpressed in fluralaner-treated groups (Additional file 14). Similarly, butene-fipronil, phoxim, cycloxaprid and chlorantraniliprole could significantly enhance the expression of CYP304A, CYP49A1, CYP4C1 and CYP9E2 in Drosophila melanogaster Meigen, B. mori, Periplaneta americana (Linnaeus) and Plutella xylostella Linnaeus, respectively [23][24][25][26]. Also, CYP4V2, CYP6A13-like and CYP12A2 were overexpressed in the insecticide resistant populations of Bemisia tabaci (Gennadius), Aphis glycines Matsumura and M. domestica, respectively [27][28][29]. Conversely, reduced CYP6B7 was found to increase larval susceptibility of H. armigera to fenvalerate [30]. Methanol was found to up-regulate the expression of CYP4D2 in D. melanogaster larvae [31] and CYP4G15 is probably important for the metabolism of endogenous compounds in the central nervous system of D. melanogaster [32]. We therefore speculate that up-regulating expression of P450 genes belonging to the CYP4, CYP6, CYP9 and CYP12 sub-families are likely involved in fluralaner metabolism in S. litura.
GSTs are a class of multifunctional detoxification enzymes that play an important role in the metabolism of insecticides [33]. The total number of GSTs genes identified in this study was less than those in the reference genome (47 GSTs) [14] but was more than those in B. mori (23 GSTs) (Additional file 6). In both fluralanertreated groups, six GST DEGs were detected and only "gene7753" belonging to the σ clade was up-regulated, to 5.10-and 4.61-fold in LC 30 and LC 50 treatment, respectively. The σ clade of GSTs play protective roles in processing xenobiotic compounds in insects [34,35]. Similar to the A. aegypti GST-2 gene, which exhibited more than 8-fold mRNA levels in resistant strains than those in susceptible strains [36], the S. litura GST-2-like gene (gene7753) was over-expressed after exposure of fluralaner (Additional file 14) indicating that GST-2 may participate in metabolizing fluralaner. CarEs can hydrolyze a diverse range of carboxylates. In T. cinnabarinus, the expression levels of CarE-6 increased up to 1.64-fold under fenpropathrin-exposure [37] and the protein level of CarE-1 was significantly increased in a chlorpyrifos-resistant strain of the small brown planthopper, Laodelphax striatellus Fallén [38]. Under LC 30 fluralaner-treatment in S. litura, venom CarE-6-like (gene8407), CarE 1C-like (gene5053) and liver CarE B-1-like (gene15885) were up-regulated by 2.64-, 4.06-and 2.57-fold, respectively (Additional file 14), which indicates that these CarEs may be involved in the metabolism of fluralaner.
Fluralaner mainly acts on ionotropic GABARs [39,40], therefore the identified ionotropic GABAR subunit genes would be useful for exploring the relationship between fluralaner and the S. litura GABAR. As shown in Additional file 10, three genes (gene14064, gene14068 and gene15324) were classified as resistance to dieldrin (RDL) subunits whilst the others were assigned to ligand-gated chloride channel homolog 3 (LCCH3) (gene11829), CG8916 (gene11828) and an unclassified clade (gene6139) (Additional file 10). However, it was notable that gene6139, annotated as a GABAR δ subunit, was possibly not a typical GABAR gene according to the phylogenetic tree analysis (Additional file 10). Similarly, "gene9312" was annotated as a GluCl subunit, "gene1980" and "gene2657" were annotated as nAChR subunits, and "gene555", "gene13723", "gene15777" and "gene4992" were annotated as AChE subunits in the NR database but not in the phylogenetic tree (Additional file 10 and Additional file 11). Therefore, we speculated that the phylogenetic analysis is very necessary while annotating the new genes. Notably, the 8916 subunit is considered as a GABAR-like gene in the phylogenetic tree analysis (Additional file 10), in agreement with a previous study in C. suppressalis [41], despite there still lacking functional evidence to date [42]. Previous studies have found that insecticide-treatment can affect expression of RDL in various insect species [43,44]. The over-expression of RDL1 (gene14064) and RDL2 (gene15324) after LC 30 and LC 50 exposure of fluralaner were similar to the relative mRNAs expression levels of GABAR from Leptinotarsa decemlineata (Say) treated with a sublethal concentration of fipronil [45]. GluCl mediates inhibitory synaptic transmission [46] and is identified as the secondary target of fluralaner [39]. Similar to many other species, such as D. melanogaster [47], Apis mellifera Linnaeus [46] and Tribolium castaneum (Herbst) [48], only one GluCl gene was identified in the S. litura transcriptome (Table 1 and Additional file 10).
Sublethal doses of fluralaner was found to reduce larval body weight, decrease pupation and emergence, and cause notched wings of adults in S. litura [9]. The DEGs of CHT 10, JHE, Tsl and Jhamt genes in the present study may mediate these effects in response to fluralaner exposure. The expression of CHT 10, JHE, Tsl and Jhamt were changed in the fluralaner-treated group (Additional file 14). As is known, CHT are important enzymes for chitin degradation and reconstruction, and play important roles in the shedding of the old cuticle and peritrophic matrix turnover in insects [49]. JHE is the primary JH-specific degradation enzyme and indirectly regulates the JH titers [50]. Tsl is not just a specialized cue for Torso signaling but also acts independently in the control of body size and the timing of developmental progression [51]. Jhamt can transfer a methyl group from S-adenosyl-L-methionine to the carboxyl group of JH acids to produce active JHs in the corpora allata in insects (Jhamt in D. melanogaster). In line with the present study, knockdown of CHT10 induced lethal phenotypes, developmental arrest and high mortality in Nilaparvata lugens (Stål) [49]. Also, depletion of B. mori JHE by CRISPR/Cas9 resulted in the extension of larval stages [50] and Tsl mutation could lead to developmental delay in D. melanogaster [51]. Furthermore, overexpression of Jhamt resulted in a pharate adult lethal phenotype in D. melangaster [52]. Therefore, the transcriptional change of these related genes may contribute to the abnormal growth and development of S. litura.

Conclusions
In conclusion, transcriptome analysis of S. litura based on the reference genome was provided in the present study. The identified DEGs may help uncover the possible molecular mechanisms underlying responses to fluralaner in S. litura. In particular, P450s may be involved in the detoxification of fluralaner in vivo. This study, therefore, may facilitate the identification of genes involved in fluralaner resistance and may offer useful information for exploring novel insecticides to control S. litura.

Insects and chemicals
A laboratory strain of S. litura was reared as described in our previous report [9]. Fluralaner (a.i. > 99%) was purified from Bravecto® (Merck & Co. Inc., Isando, South Africa) as described in Sheng et al. (2019) [53]. The 3rd instar larvae of S. litura were chosen for the experiment and three treatments (45 larvae per treatment) were used in the present study. Larvae in the control group were fed with artificial food containing acetone and in the fluralaner-treated groups with artificial food containing different sublethal concentrations [LC 30 (0.59 mg of fluralaner per kg of artificial food, mg·kg − 1 ) and LC 50 (0.78 mg·kg − 1 ), respectively] of fluralaner dissolved in acetone. After 24 h, 15 alive and normal S. litura from each treatment were randomly collected, equally divided into three 1.5 mL Eppendorf tubes, respectively, immediately frozen in liquid nitrogen and stored at − 80°C until use.
Total RNA extraction, cDNA library construction and IHX sequencing Total RNA was isolated from whole body of S. litura from each group (the control and fluralaner-treated groups) by TRIzol™ Reagent (Invitrogen, Carlsbad, CA). The concentration and integrity of total RNA were measured using NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA) and with the RNA Nano 6000 Assay Kit by Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA), respectively. RNA sequencing libraries were generated from each sample and the sequencing was carried out by Bio-Marker (Beijing, China) following the manual of the NEB Next® Ultra™ RNA Library Prep Kit (NEB, Ipswich, MA) in Illumina. Briefly, 1 μg of total RNA from each sample was enriched using magnetic beads with Oligo (dT) and broken into short fragments, which was carried out using divalent cations under elevated temperature by adding NEB Next First Strand Synthesis Reaction Buffer (5x) (NEB). These short fragments of mRNA served as templates for synthesis of the first-strand complementary DNA (cDNA) with random hexamer primers and Moloney Murine Leukemia Virus (M-MuLV) Reverse Transcriptase (NEB). Using buffer, dNTPs, RNase H and DNA Polymerase I, the second-strand cDNA was subsequently synthesized. Remaining overhangs were converted into blunt ends via exonuclease/ polymerase activities. After adenylation of 3′ ends of DNA fragments, NEB Next Adaptor with hairpin loop were ligated for hybridization. The purified cDNAs were subjected to end repair by adding an adenosine triphosphate (A) base to the 3′ end and connected with sequencing adaptor. Suitable fragments of cDNAs were extracted by the AMPure XP system (Beckman Coulter, Beverly, CA). Three microliter of USER Enzyme (NEB) was used for size selection, and adaptor-ligated cDNAs at 37°C for 15 min followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High Fidelity DNA polymerase (NEB), Universal PCR primers, and Index (X) Primer (NEB). PCR products were purified with the AMPure XP system (Beckman Coulter) and library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent).
The clustering of the index coded samples was carried out on a cBot Cluster Generation System (Illumina) using TruSeq PE Cluster Kit v4-cBot-HS (Illumina) according to the manufacturer's instruction. After cluster generation, the prepared libraries were sequenced on an Illumina platform and paired end reads were generated.

Mapping and annotation of the clean reads
Raw reads in the fastq format were first processed by inhouse perl scripts. Subsequently, clean reads were obtained from raw reads by removing reads with low quality, adaptor and poly-N. The intergenic regions and introns were identified by comparing the reads to the general feature format (GFF) file of the reference genome [14]. The clean reads with a perfect match or with only one mismatch were mapped to the reference genome of S. litura [54] by the Hisat2 software [55]. The mapped reads were assembled into the transcriptome by the StringTie method [56] and all the genes were functionally annotated and classified by Blast [57] according to the databases of Clusters of Orthologous Groups (COG) [58], Gene Ontology (GO) [59], Kyoto Encyclopedia of Genes and Genomes (KEGG) [60] and Non-redundant protein (NR) [61]. In addition, new genes from the unannotated sequences encoding more than 50 amino acids were identified and analyzed by Blast using the above mentioned databases.
Referring to the previous studies [10,21,62], all of the confirmed protein sequences were aligned with S. litura, B. mori, D. melanogaster and A. mellifera, then the phylogenetic trees were mapped with 1000 bootstrap replications using the neighbor-joining method to evaluate the branch strength of each tree using MEGA 7 [63]. The phylogenetic trees were further annotated via the Evol-View online tool [64]. Referring to previous publications [46,48,65], the phylogenetic tree of cys-loop ligand-gated ion channel proteins was constructed using Clustal X [66] and TreeView software [67].

Identification and validation of DEGs
The relative expression levels of genes among the control, LC 30 -and LC 50 -treated groups were analyzed using DEGseq method [68], which provided statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The fold change (FC > 2) and false discovery rate (FDR < 0.01) from corrected P-value were used as a judgment standard for the significant differences in gene expression between two samples [15]. The DEGs were annotated, classified and analyzed by databases of COG, GO, KEGG, and NR. The RT-qPCR assay was performed to validate the randomly selected 26 DEGs using TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) (Takara Biotechnology (Dalian) Co., Ltd., Liaoning, China) on a Quant Studio™ 6 and 7 Flex Real-Time PCR System (Life Technologies Corporation, Carlsbad, CA) according to the procedures from Liu et al. (2018) [9]. The housekeeping gene of EF-1α was selected as an internal standard to eliminate variations in mRNA and cDNA quantity and quality [69]. Three technical replications per biological replicate were conducted and the relative mRNA expression levels of genes were analyzed using the 2 -△△Ct method [70]. The primers used in this study are listed in Additional file 15. All values were expressed as means ± standard error (SE) and the statistical figures were constructed by Microsoft Excel (Microsoft, Redmond, WA).