Transcriptome profile analysis reveals the emamectin benzoate-induced genes associated with olfaction and metabolic detoxification in Spodoptera exigua Hübner (Lepidoptera: noctuidae)

Abstract

Emamectin benzoate (EMB) is one of the high-efficacious pesticide commonly used to control the ferocious economic crop pest S. exigua in the field. The results of S. exigua life table treated with three lethal doses of EMB (LC₅, LC₂₀ and LC₅₀) indicated that three lethal doses of EMB could all inhibited the proliferation of S. exigua population effectively by long-lasting inhibit the survival of the S. exigua larvae and per female fecundity significantly (p < 0.05), and the inhibition effects of key life parameters were affected by EMB stress in a dose-dependent manner. 67,123 unigenes with an N50 value of 1941 bp were assembled by RNA-seq data analyses and 44,113 were expressed more than 12 samples. Furthermore, 3517, 3139 and 5135 differentially expressed genes were identified under different lethal doses of EMB treatment. Finally, 16 hub genes and three olfaction-related genes were identified by WGCNA analysis as biomarkers for the EMB response, and the cAMP-dependent proteins, P450 metabolic-related genes and ABC transporter genes might have been involved in response to EMB stress in S. exigua.

KEYWORDS:

• Spodoptera exigua
• emamectin benzoate
• RNAseq
• WGCNA
• pesticide stress
• olfaction-related genes

Introduction

The beet armyworm, Spodoptera exigua (Hübner) which originated from south-east Asia has now become an important worldwide agricultural pest (Hafeez et al. 2020). As a polyphagous pest, S. exigua has been identified to be detrimental to over 130 host plants species, as it causes considerable losses to economic crops, such as Chinese cabbage, soybean, peanut, onion, corn, scallion, ginger, pepper, sorghum, tobacco and cotton (Moulton et al. 2002; Zheng et al. 2011; Zhou et al. 2011; Hafeez et al. 2021).

The key stage for the effective control of S. exigua is before the 3rd instar stage, because the larvae at this stage are small, have poor drug resistance, and are easy to kill. After the 3rd instar, the resistance of larvae against pesticides increases. S. exigua also exhibits a phenomenon known as false death, which makes this pest difficult to control (Zhang et al. 2011). In field, various control measures have been applied to deal with S. exigua damage, such as biological, physical, chemical as well as mechanical control (Akter et al. 2019). Among these methods, chemical control, which relies on the spraying of different pesticides such as carbamate (methomyl), pyrethroid (fenpropathrin), organophosphate (chlorpyrifos), and biopesticide (Emamectin benzoate, abbreviated as EMB), is the major strategy (Liburd et al. 2000; Ishaaya et al. 2002; Huang et al. 2021).

Currently, chlorantraniliprole, deltamethrin, cyhalothrin, indoxacarb, pyrethroid, spinetoram, and EMB are widely used to control the S. exigua in China (Wang et al. 2019b; Lai et al. 2011). Due to the excessive and frequent application of pesticides, S. exigua has developed different resistance levels against most conventional pesticides around the world (Hartman and Clair 1999). For example, manual selection of S. exigua carried out in the laboratory with chlorantraniliprole resulted in an increase in resistance by 11.8-fold after continuous screening for 23 generations, and by 128.6-fold compared to the susceptible strains kept in the laboratory (Lai et al. 2011). Another study indicated that S. exigua developed a 309.5-fold resistance to lambda-cyhalothrin, and selection in the laboratory resulted in a 423-fold increase in resistance (Tian et al. 2018). Since S. exigua has developed high resistance to these conventional pesticides, there is an urgent need for highly effective pesticides in the field.

Previous studies have shown that among the different conventional insecticides, EMB is an excellent pesticide that can be used to control S. exigua in the field efficiently (López et al. 2010; ElSheikh 2015; Moustafa et al. 2016). As a widely used macrocyclic lactone insecticide, EMB ((4^″R)-4^″-deoxy-4^″-(methylamino) avermectin B1 benzoate) is a biopesticide extracted from the naturally-occurring compound, ivermectin, that has high efficiency in controlling many economic crop pests (Mushtaq et al. 1996; López et al. 2010; Mokbel and Huesien 2020). The action mode of EMB is as a glutamate-gated chloride channel allosteric modulator (IRAC 2019), which disturbs the nervous and muscular system of insects, ultimately resulting in their deaths (Ishaaya et al. 2002; Zaka et al. 2014). EMB possesses excellent insecticidal potency, which can be up to 80,400 more effective compared to conventional pesticides such as fipronil, indoxacarb, cypermethrin, and fenvalerate, used on economic crops (Argentine et al. 2002; Ishaaya et al. 2002; López et al. 2010).

Insects mainly rely on the olfactory system to identify a variety of environmental chemicals, including pesticides and other exogenous substances. A study reported that the binding affinity of two A. lepigone pheromone binding proteins (AlepPBP2 and AlepPBP3) to phoxim is higher than that of other insecticides (Zhang et al. 2020a). After the amino acid residues of the protein were changed by site-directed mutagenesis, the binding affinity of AlepPBP2 and AlepPBP3 to phoxim decreased, indicating that they should play a key role in the AlepPBPs/phoxim interaction. Therefore, the relationship between the effect of pesticides on insects and the sense of smell of insects is worthy of further study. ABC transporters exist in all areas of life and have the function of natural excretion. Ryan C. Gott et al. found that it is very reasonable for insects to use these transporters to excrete insecticides and produce drug resistance (2017). Studies have shown that several members of the insect ABC transporter superfamily play a very important role in pesticide resistance by inhibiting the accumulation of intracellular pesticides and their metabolites. Several members of the insect ABC transporter superfamily inhibit the accumulation of intracellular pesticides and their metabolites. It plays a very important role in pesticide resistance (Buss and Callaghan 2008).

The resistance to EMB due to its extensive application, has been documented in different pests, such as Amrasca devastans (Argentine et al. 2002), Musca domestica (Abbas et al. 2015; Khan et al. 2016), Spodoptera litura (Shad et al. 2010; Zaka et al. 2014; Wang et al. 2018), and also S. exigua (Ishtiaq et al. 2014; Abbas et al. 2015). A field study conducted from 2008 to 2010 in three regions in Shandong province (China) revealed that S. exigua had developed varying degrees of resistance to EMB. Yet, the molecular mechanism underlying the evolution of this resistance in S. exigua against EMB remains unclear (Zhou et al. 2011). Therefore, in order to systematically identify the target genes in S. exigua responsible for its resistance to EMB, this study used transcriptome sequencing, bioassay, and correlation analysis to elucidate the underlying molecular pathways that are involved in the responses of S. exigua treated with different lethal doses of EMB.

Materials and methods

Tested insects

EMB has both gastric toxicity and contact killing properties as its mode of action (MoA). The detailed information about its purity (95.5% AI,) was provided by Hebei Veyong Bio-Chemical Co. Ltd. (Source) (Sun et al. 2016, 2014). A regression equation, y = 5.8037 + 1.3285x, which described the responses of the 3rd instar larvae of S. exigua to EMB toxicity was developed in our previous study. Also determined was the final 50% lethal concentration (LC₅₀) of EMB, which was = 0.248 mg•L^–1 (SE = 0.04, the 95% confidence interval was 0.18–0.35 mg•L^–1, and the R-value was 0.9971). The formula given above was used to prepare different lethal doses of EMB (LC₅ = 0.0144, LC₂₀ = 0.0578, and LC₅₀= 0.248 mg•L^–1) which were incorporated into artificial diets to feed the 3rd instar larvae of S. exigua. Acetone and tween80 solution treatments were used as solvent contrast control (AC), while water was used as a blank control (CK).

Pre-experimental results indicated that concentrations lower than LC₅, had little effects and concentrations higher than LC₅₀, produced no results for most life parameters. Therefore, the lethal concentrations, LC₅, LC₂₀, and LC₅₀ were chosen for this study. Samples made up of the 3rd instar larvae of S. exigua of similar sizes were selected and fed on artificial diets mixed with different doses (LC₅, LC₂₀, and LC₅₀) of EMB, AC and CK. Each treatment was replicated six times, each repeat had 30 dishes and each dish had five larvae, resulting in every treatment containing at least 900 larvae of S. exigua.

The original colony of S. exigua was provided by the Plant Protection Institute of Jiangsu Academy of Agricultural Sciences (China). S. exigua larvae were reared on an artificial diet and did not have any contact with pesticides in the climate room. The room was set at 27 ± 1°C, RH 65 ± 5% humidity and 14: 10 d/night light regime. A 10% honey water was supplied to the adults (Dai 2015; Zhao et al. 2018). After feeding on the diet mixed with EMB for 72 h, surviving insects were switched to a non-toxic fresh artificial diet. Data on larval mortality, oviposition and hatching ability were recorded following the methods described by Zhao et al. (2018) and Sun et al. (2016) with slight modifications. For example, within 72 h from the beginning of the experiment, individual insects which did not move when their bodies were touched repeatedly with a brush dipped in water were considered dead. Finally, we also considered insects that did not successfully pupate as dead. We defined ultimate mortality as the mortality of larvae from the beginning of treatment to the end of the 5th instar larvae, that is, the total mortality of the larvae under different concentrations of EMB. In addition, the following formula was used to calculate ‘growth inhibition rate’: [(WAT_CK–WBT_CK)–(WAT_LC–WBT_LC)]/(WAT_CK–WBT_CK). (WAT_CK is weight after treatment under blank control; WBT_CK is weight before treatment under blank control; WAT_LC is weight after treatment under lethal concentration treatment; WBT_LC is weight before treatment under lethal concentration treatment). Data analyses were carried out with the SAS v.9.0 software (SAS Institute, Cary, NJ, USA). Data were mean ± standard deviation and the differences between means were analyzed using a one-way analysis of variance (ANOVA), followed by Tukey’s honest significant difference (HSD) test. Samples of S. exigua fed with the diet mixed with EMB for 72 h were collected for later transcriptome sequencing and quantitative RT–PCR analysis. The sample from each treatment contained at least 100 larvae of S. exigua.

RNA isolation and sequencing

Total RNA was extracted from S. exigua larvae treated with different concentrations of EMB (CK, AC, LC₅, LC₂₀, LC₅₀), with the SV Total RNA Isolation kit (Promega, Madison, WI, USA) following the manufacturer’s instructions (Tian et al. 2018). Magnetic beads with Oligo (dT) were used to enrich eukaryotic mRNAs to construct the sequencing library, after which a fragmentation buffer was added to break the mRNAs into short pieces. Six random base primers were used to synthesize the first strand of cDNA using mRNA as templates; after which the double-strand cDNA was synthesized by adding buffer, dNTPs, DNA polymerase I and RNase H. The AMPure XP beads were used to select the fragment size. Finally, PCR amplification was performed, and the AMPure XP beads were used to purify the PCR product to obtain the final library.

Qubit 2.0 was used to initially quantify the constructed library and was diluted to 1.5 ng•μl^–1. An Agilent 2100 bioanalyser was used to detect the insert size of the library. The Q-PCR method was used to quantify the effective concentration of the library (> 2nM) after achieving the expected insert size (Wang et al. 2019). Then, the different libraries were pooled for sequencing on the Illumina HiSeq platform, based on their effective concentration and the required raw data.

Quality control and functional annotation of unigenes

Adapters and low-quality reads were first removed using fastp (Chen et al. 2018). Using the clean reads, a de novo transcriptome assembly was performed using Trinity, and the longest transcripts of genes were selected as unigenes (Grabherr et al. 2011). The CDS sequences were generated using TransDecoder and then were aligned to protein databases, including the NCBI non-redundant protein database, Swiss-Prot database, KOG and KEGG database using Blastx with E-values < 1E-5. GO annotation of the unigenes was performed using Blast2GO and InterProScan5 software (Conesa et al. 2005; Jones et al. 2014).

Analysis of differentially expressed genes (DEGs)

For each sample, high-quality reads were mapped to unigenes using Tophat and Bowtie (Langmead et al. 2009; Trapnell et al. 2009; Kim et al. 2013). Cufflinks calculated the expression level (FPKM, fragment per kilobase per million reads). Differentially-expressed genes (DEGs) between groups were identified using cuffdiff with FDR (false discovery rate) <0.01 and a 2-fold change (Roberts et al. 2011).

GO and KEGG pathway analysis

Based on the GO and KEGG annotation results, a hypergeometric test was used to analyze the functional enrichment of DEGs. The dhyper function in the R language was used to calculate the p-value in the enrichment analysis. Also, the multiple-test corrected p-value (q-value) was calculated using the Benjamini-Hochberg function implemented in the R library. Q value <0.05 was established as the significance threshold.

Co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was used to construct a scale-free network from gene expression data (Langfelder and Horvath 2008). Also, a hybrid dynamic tree-cutting algorithm was used to allocate modules, which included at least 15 genes. After studying the correlation between the gene expression modules and phenotype characteristics, the correlation coefficient, and the p-value of correlation analysis were used to determine whether the module was related to the phenotype. Also, an MM (module membership) value was used to identify the hub genes in the modules. The Cytoscape software was used to visualize these networks (Killcoyne et al. 2009).

qRT–PCR data analysis

Total RNA of S. exigua samples treated with different concentrations of EMB (CK, AC, LC₅, LC₂₀, LC₅₀) was extracted using the SV total RNA isolation kit (Promega, Madison, WI, USA) following the manufacturer's instructions. An MMLV Reverse Transcriptase (Promega, Madison, WI, USA) treated with ribonuclease H (TCytoscapeakara, Tokyo, Japan) was used to synthesize the cDNA library, and spectrophotometry was used to quantify the library. The primers used for real-time quantitative reverse transcription (RT)-PCR (SYBR Green I) were designed based on eight randomly-selected genes from the RNA-seqdata (Table 1). Also, the S. exigua housekeeping gene, β-actin, was used as an endogenous reference for data normalization (GenBank accession no. JN616391). The primers for β-actin were 5′-CCAGCCTTCCTTCTTGGGTAT-3′ and 5′- AGGTCCTTACGGATGTCAACG-3′ (Zhao et al. 2018). The other 15 random sequences that were used to verify differential transcriptome sequence results in this study are listed in Table 1.

Table 1. Fifteen RT-Q-PCR primer sequences used to verify differential transcriptome sequencing results in this study.

Gene Code	Upstream primer (5’-3’)	Downstream primer (5’-3’)
1 c10037_g1	5'-GATAAGCCCATGAGGAATG-3'	5'-TAAAGACGCCGAGAACAC-3'
2 c10097_g1	5'-ATCGCCATGATACGCGTG-3'	5'-GAACCTTAACATTTCCAAA-3'
3 c10172_g1	5'-CAGTGAGGCAGTGATGTTTAT-3'	5'-ACCATGTCCAATCCGTGT-3'
4 c10441_g1	5'-AAAAATGCGAGCGAATCAAC-3'	5'-TGTGCCGACGAAGACGAG-3'
5 c10014_g1	5'-GAGATGGGAGTTGAACAGTC-3'	5'-CCAAAGTATGCCAAGAAAT-3'
6 c10095_g1	5'-GGCAGCCTCTTCAT65ACATC-3'	5'-GGTTTCCGTTATTATCCTCC-3'
7 c10035_g1	5'-TTGAACCTGGAACACGA-3'	5'-GCACTACGAGCAACCTCT-3'
8 c10291_g1	5'-CATACGTCGGTTCACGC-3'	5'-TTTCTGGCATAGGCTTCA-3
9 c10010_g1	5'-TTTTATACCGTGCCTTTAC-3'	5'-ACACTGGGTCACTTTGC-3'
10 c10070_g1	5'-ATGTCCACAGGAACACGCAATG-3'	5'-AGGGTCACAGGCCCAAACG-3'
11 c10160_g1	5'-CTAGACTATACTGCTGTAGATTT-3'	5'-TTTCTTACTCCCAACGA-3'
12 c1007x_g1	5'-TCAAGGCAGCCGATGAT-3'	5'-CAGGGATTGGTTGTGGTTC-3'
13 c10310_g1	5'-CATAGAGCAACCCAAGTAG-3'	5'-GACATCGGTCATCATCAA-3'
14 c10510_g1	5'-CACCTGCTCGTAGAAGTCG-3'	5'-GCCTCATTGTCCCATAGTT-3'
15 c10188_g1	5'-AACGGACTCTGAGATAGG-3'	5'-CATAGGCACTTTGCTGTA-3'

Primer 5.0 software was used to design all EST-specific primers. The reactions were carried out using the SYBR Premix Ex Taq Kit (Takara, Tokyo, Japan) and were run on the Bio-rad iCycler real-time quantitative RT–PCR detection system. Melting curve analysis and gel electrophoresis were also used to confirm these genes. The reactions of each treatment were replicated four times; non-template control reactions were performed in triplicate for each primer pair. The expression levels of each gene in S. exigua under the different treatments (expressed as the relative quantification (RQ) values), were calculated using the 2^−ΔΔCt method (Livak and Schmittgen 2001). Significant differences in gene expressions were detected with Duncan’s new multiple range method.

Results

Physiological indices of S. exigua under different treatments

After S. exigua larvae were treated with different doses of EMB, the S. exigua all showed a series of long-lasting effect phenotypes, such as non-feeding, growth retardation, slower pupal development and reduction in oviposition. Specifically, the growth inhibition rate, the final mortality rate and the mortality rate after 72 h of larvae treated with different lethal doses of EMB were all significantly higher compared to the CK and AC controls (p < 0.05). By contrast, larval weight was significantly lower than that in the control groups (p < 0.05), and no significant difference was found between the two controls (p > 0.05) (Table 2). The mortality after 72 h treatment of EMB was close to the 5%, 20%, and 50% mortality rates predicted by the LC₅, LC₂₀, and LC₅₀ respectively. However, the LC₅, LC₂₀, and LC₅₀ values were higher than 5%, 20%, and 50% of final mortality rates recorded respectively. This indicated that the EMB treatment had a sustained lethal effect on the S. exigua larvae; specifically, an inhibitory effect on the development and survival of the S. exigua larvae.

Table 2. Effects of different lethal doses of EMB on key life-history parameters of S. exigua larvae.

	Different treatment
Key life history parameters of S. exigua larvae	CK	AC	LC5	LC20	LC50
Weight before treatment (mg)	7.494 ± 0.203a	7.436 ± 0.247a	7.410 ± 0.247a	7.500 ± 0.226a	7.535 ± 0.223a
Weight after treatment (mg)	34.736 ± 1.266a	34.536 ± 0.945a	28.238 ± 0.731b	19.308 ± 1.026c	13.836 ± 1.355d
Growth inhibition rate (%)	0.000 ± 0.000a	0.460 ± 2.787a	23.510 ± 1.553b	56.681 ± 1.953c	76.942 ± 3.528d
Mortality after 72 h treatment (%)	2.833 ± 0.950a	3.000 ± 1.395a	8.333 ± 1.667b	26.667 ± 4.125c	51.500 ± 6.303d
Ultimate mortality (%)	3.500 ± 1.087a	3.833 ± 1.728a	26.170 ± 6.734b	45.500 ± 5.157c	78.000 ± 6.169d
Larval period (d, start at 3rd)	8.372 ± 0.919b	8.428 ± 0.952b	9.440 ± 1.412b	11.538 ± 2.772a	12.116 ± 3.081a

Note: Data are mean ± standard error; Different lowercase letters indicate a significant difference obtained by ANOVA, followed by a Tukey’s HSD test (p < 0.05).

The life span of S. exigua significantly increased (p < 0.05) following an increase in lethal doses (LC₅ to LC₅₀) compared to the CK and AC controls, while the eclosion rate of S. exigua significantly decreased (p < 0.05) (Table 3). Similar downward trends were recorded for egg hatching rates, intrinsic rate of increase, per female fecundity following increases in lethal doses (LC₅ to LC₅₀). The intrinsic rate of increase and per female fecundity showed significant decreases after treatment with different lethal doses of EMB (p < 0.05) compared to the controls (Table 3). There were no significant differences in egg hatching rate, larvae and total periods between the S. exigua treated with EMB at LC₅ dose and two controls (p > 0.05). However, these were significantly lower (egg hatching rate) or higher (larvae and total periods) in S. exigua treated with EMB at LC₂₀ and LC₅₀ doses than that in two controls or S. exigua treated with EMB at LC₅ dose (p < 0.05). In addition, there was no significant difference in all of the key life-span parameters between the blank (CK) and solvent (AC) controls (p > 0.05) (Tables 2 and 3). These results indicated that key life parameters, such as larval final mortality rate, per female fecundity and intrinsic growth rate in S. exigua were affected by EMB stress in a dose-dependent manner; and different lethal doses of EMB could long-lasting inhibited the development and survival of the S. exigua larvae, and significantly reduced per female fecundity and intrinsic growth rate of S. exigua (p < 0.05), therefore effectively inhibited the proliferation of S. exigua population.

Table 3. Effects of different lethal doses of EMB on key life-history parameters of S. exigua.

	Different treatment
Key life history parameters of S. exigua	CK	AC	LC5	LC20	LC50
Life-span (d, start at 3rd)	25.306 ± 2.011b	25.182 ± 1.920b	26.836 ± 3.107b	29.936 ± 4.911a	31.000 ± 5.567a
Adult emergence rate (%)	78.89 ± 5.70a	79.25 ± 4.13a	71.39 ± 3.57b	65.92 ± 5.82c	60.30 ± 7.44d
Per female fecundity	613.45 ± 32.10a	601.85 ± 29.90a	467.46 ± 22.36b	313.48 ± 39.82c	246.82 ± 27.79d
Hatching rate of eggs (%)	70.99 ± 2.30a	70.22 ± 5.56a	72.29 ± 4.66a	65.69 ± 2.79b	56.75 ± 2.46c
Intrinsic rate of increase (R)	0.250 ± 0.033a	0.247 ± 0.049a	0.171 ± 0.028b	0.146 ± 0.029c	0.093 ± 0.018d

Note: Data are mean ± standard error; Different lowercase letters indicate a significant difference obtained by ANOVA, followed by a Tukey’s HSD test (p < 0.05).

De novo assembly of transcriptome and functional annotation of unigenes

The Illumina HiSeq^TM 4000 sequencing platform was used to construct transcriptome sequence containing 15 cDNA libraries. Consequently, a total of 90 Gb of 150-bp paired raw reads were generated in this study. De novo transcriptome assembly was performed using trinity which generated 67,123 unigenes with an N50 value of 1941 bp.

A total of 26,358 unigenes were annotated in at least one public database, such as NR, NT, COG, GO, KEGG, Swissport, and Interpro protein database (Additional file 1). Out of this total, 16,654 unigenes were annotated to gene ontology terms. These unigenes were found in 49 functional groups: 23 groups in the ‘biological process’, 16 groups in the ‘cellular component’, and ten groups in the ‘molecular function’. In addition, the top five species from which the most BLAST hits were obtained included Bombyx mori, Danaus plexippus, Microtus ochrogaster, Helicoverpa armigera, and Tribolium castaneum (Figure 1(A)).

Figure 1. (A) Species distribution of unigenes annotation. (B) KEGG annotation of unigenes.

In the KEGG database, a total of 8591 unigenes were assigned to 32 functional categories; the largest category was signal transduction (971 unigenes), followed by carbohydrate metabolism (685 unigenes), translation (654 unigenes), and then other pathways (Figure 1(B)).

Differential expression and GO enrichment analysis

To determine the expression of S. exigua unigenes at different EMB treatments, the clean RNA-seq reads were mapped to the 67,123 unigenes. Gene expression levels were calculated as FPKM values. Finally, we found that 44,113 of 67,123 expressed in at least 12 samples. In this study, differential gene expression analysis was conducted between the different EMB treatments and the control (AC and CK) with fold change levels ≥2 with a FDR ≤0.05. Removal of the DEGs between the water (CK) and the solvent controls (AC) resulted in the identification of 3517 DEGs between AC and LC_5. Of these 1,622 unigenes were down-regulated, while 1895 genes were up-regulated under LC₅ treatment (Additional file 7). In addition, 3139 DEGs were identified between the AC and LC₂₀. Of these, 1644 genes were down-regulated, while 1495 were up-regulated following the LC₂₀ treatment (Additional file 8). Also, 5135 unique DEGs were found between the AC and LC₅₀ treatment. Out of these 2351 were down-regulated, while 2784 were up-regulated following the LC₅₀ treatment (Figure 2(A), Additional file 9). More importantly, the number of DEGs was the highest in samples treated with the LC₅₀ of EMB.

Figure 2. (A) the up-regulated genes compared to down-regulated genes in LC after removal of the differentially-expressed genes between AC and CK (B) The enrichment GO terms of DEGs.

To understand the functions of the DEGs, the GO terms of DEGs were performed by GO enrichment analysis. As shown in Figure 2(C), ‘growth’ (GO: 0040007) was enriched between LC₅₀ and AC. This indicated that the growth of S. exigua was inhibited by the EMB treatment. We also found that the EMB treatment can cause concentration-specific gene expression patterns. GO enrichment analysis results indicated that the DEGs were also related to ‘carbohydrate metabolic process’ (GO: 0005975), ‘cell proliferation’ (GO: 0008283), and ‘cellular metabolic process’ (GO: 0044237), but these GO terms were enriched in different EMB treatment concentrations.

Devlopment and immunity-related pathway analysis by KEGG

KEGG pathway enrichment analyses was carried out to understand the functions of DEGs. A total of 21 KEGG pathways were enriched in the DEGs between LC₅ and AC. One of these was the peroxisome pathway which is related to insect development and immunity. It is a remarkbale fact that KEGG pathways of ‘apoptosis’, ‘ABC transporters’, ‘RNA degradation’, ‘olfactory transduction’, ‘folate biosynthesis’ and ‘lysosome pathways’ were significantly enriched in DEGs between LC₂₀ and AC. These pathways are responsible for insect development, reproduction and immunity. As the Figure 3 shows, the Lysosomal-associated membrane protein 1 (LAMP-1) gene was up-regulated under EMB treatement and three cathepsins genes were up-regulated under EMB treatement. Also, 17 pathways which including the ribosome, folate biosynthesis, and phagosome pathways were enriched in the DEGs between LC₅₀ and AC. These are related to insect development and immunity. These results indicated that the regulation of the expression of these genes involved in insect development and immunity could be negatively affected by lethal concentrations of EMB.

Figure 3. Lysosome pathway and the expression of the related DEGs. In the heatmaps, there are five samples in each row, which represent the expression level of the gene in the five samples of AC, CK, LC₅, LC₂₀, LC₅₀. Red represents a high level of expression, while green represents a low level of expression.

Effects of gene expression levels on ABC transporters and insect olfaction

KEGG pathway analysis of the transcriptome data showed that the ABC transporter was an important gene family involved in the response of the S. exigua to different lethal stress of EMB (Figure 4). The ABC transporter-related genes c15483_g1(ABCB4, ATP-binding cassette subfamily B member 4), c15483_g2(ABCB1, ATP-binding cassette subfamily B member 1) and c17259_g2(ABCA3, ATP-binding cassette subfamily A member 3), were highly expressed under LC₂₀ and LC₅₀ treatment (p < 0.05), and only c1698_g1 and c22475_g1 were highly expressed under LC₅₀ treatment (p < 0.05).

Figure 4. ABC transporter pathway and the expression of the related DEGs. In the heatmaps, there are five samples in each row, which represent the expression level of the gene in the five samples of AC, CK, LC₅, LC₂₀, LC₅₀. Red represents a high level of expression, while green represents a low level of expression.

Olfaction has an important role in many aspects of insect behavior, such as foraging. The analyses revealed several DEGs involved in the olfactory transduction pathway under EMB treatment (Additional file 7-9). One of such was the cAMP-dependent protein kinase catalytic subunit-like gene that was up-regulated in the S. exigua larvae treated with EMB. While the cGMP-dependent protein kinase G gene was down-regulated under EMB treatment. Also, the calmodulin gene, c21765_g1, in the olfactory transduction pathway and cAMP-dependent protein kinase catalytic subunit gene, c20528_g1, were all highly expressed under LC₅₀ treatment (p < 0.05), and only the calmodulin gene, c4625_g2 was highly expressed under LC₅ treatment (p < 0.05).

Weighted correlation network analysis and identification of hub genes

Furthermore, we performed weighted gene co-expression network analysis (WGCNA) to explore the dynamic patterns of S. exigua gene expression and identify genes associated with complex phenotypes. WGCNA divided the DEGs into 42 different modules that depicted by a branch of different colors (Figure 5(F)). Association of the detected modules and the physiological indices showed that five modules (orange, blue, royal blue, skyblue3 and light yellow), containing 668 unigenes, exhibited a strong correlation between gene expression and some key life-history parameters (|correlation coefficient| > 0.6, p < 0.01), such as growth inhibition rate, larval period, mortality after 72 h treatment, ultimate mortality, life-span, adult emergence rate, per female fecundity, hatching rate of eggs and intrinsic rate of increase. These suggested that these DEGs played an important role in S. exigua when in response to EMB stress.

Figure 5. (A–E) Network of genes in the blue, yellowlight, skyblue3, orange and royalblue modules. (F) Association of modules with key life-history parameters. Forty-two modules were identified by WGCNA analysis. For each module, a module eigengenes (ME) was selected as representative of the module gene expression. The correlation coefficients and p values were marked on the diagram.

Hub genes, defined as genes which are mostly connected to other genes in the modules, are more important than other genes in the cluster. According to this principle, four, three, three, three, and three hub genes in orange, blue, royal blue, skyblue3, and lightyellow modules were identified, respectively (Figure 5(A–E), Additional file 2-6). For further elucidation of the functional properties of the genes in important modules, KEGG pathway enrichment analysis was also conducted (Figure 6). The results indicated that genes in the orange module were enriched in protein digestion and absorption and carbohydrate digestion and absorption pathways, which indicated that the food digestion and absorption of S. exigua were affected by EMB treatment. The yellow light module was positively related to the mortality after 72 h treatment and negatively correlated with the hatching rate of eggs. The genes in lightyellow module were enriched in the mRNA surveillance pathway, fat digestion, absorption pathway and drug metabolism – cytochrome P450. This module included two P450 metabolic-related unigenes, which may played an important role in the regulation of detoxification function. Unigene c13436_g1 was significantly down-regulated under LC₂₀ and LC₅₀ treatment (p < 0.05), and unigene c21930_g4 was significantly up-regulated under LC₅ treatment (p < 0.05). In addition, the blue module, which was significantly correlated with seven key life-history parameters, was mainly enriched in pathways apoptosis, toll-like receptor signaling pathway, and DNA replication.

Figure 6. KEGG pathway enrichment analysis of gene in (A) orange, (B) skyblue3, (C) blue, (D) royal blue, and (E) light yellow.

qRT-Data analysis

A total of eight genes were selected according to their significance of differentially expression to confirm the gene expression data generated from the quantitative RT–PCR results. They showed that the c10037 was responsible for one kind of anti-peronosporomycetal agent synthesis; the c10172 and c10095 genes were responsible for 4-hydroxyproline betaine 2-epimerase (hpdD) and norathyriol-6-O-beta-D-glucoside (phytochemical compounds, br08003, https://www.genome.jp/kegg/kegg2.html) synthesis; the c10291 and c10070 genes had critical roles in alizarin 2-methyl ether (br08003) and isogentisin (br08003) synthesis; the c10160 and c10310 genes had important roles in myosmine (br08003) and chimaphilin (br08003) synthesis; the c10441 gene was responsible for 4-p-coumaroylquinic acid synthesis. The transcriptional changes in these eight genes between RNA-seq and RT-qPCR are described in Figure 7(a–h). A significant correlation (R² = 0.79) was also observed between the results of RT-qPCR and the RNA-seq analysis (Figure 7(i)), indicating the reliability of the RNA-seq data results.

Figure 7. In total, eight genes were randomly selected to confirm the gene expression data generated from the quantitative RT-PCR results.

Discussion

EMB is a highly efficient pesticide that is commonly used for pest control in the field (López et al. 2010; Moustafa et al. 2016). However, S. exigua has begun to develop resistance to EMB. A previous study revealed no significant difference in the key life parameters (such as spawning rate or larvae mortality) between a conventional strain with EMB resistance and a normal strain (p > 0.05) (Che et al. 2015). Simultaneously, the D value calculated from LC₅₀ of the two-parent strains and their F1 progeny was 0.6, which confirmed incomplete dominance of resistance to EMB in the resistance strain. It also indicated the return of the physiological indices to normal levels, even when the resistance of S. exigua against EMB appeared (Che et al. 2015). However, in this study, treatment with different lethal doses of EMB resulted in a significantly higher larval mortality and a prolonged larval stage (ultimate mortality of larvae stage) (p < 0.05). They also had a significant negative effect on adult parameters, such as female fecundity, adult emergence rate and egg hatching rates (p < 0.05). These showed that EMB had a great long-lasting control effect on conventional stains of S. exigua. However, the resistant strains of S. exigua may have evolved physiological or molecular mechanisms to overcome the toxicity of EMB (Che et al. 2015).

Similar to our results, other reports also indicated that EMB had a significant impact on the growth and spawning of pests. For example, Muhammad et al. (2018) reported a significant negative effect on the development, feeding potential and body weight in Paederus fuscipes treated with LC₃₀ of EMB. Also, treatment of Spodoptera littoralis with EMB resulted in a significant reduction in body weight and decreased survival rate (ElSheikh 2015). Further, treatment of Mamestra brassicae with different doses of EMB resulted in prolonged development, sustained mortality, and decreased egg production (Moustafa et al. 2016). The above results indicated that pests treated with EMB showed similar inhibitory phenotypes in growth, development, and population reproduction. In addition to EMB, another study showed that the effect of inhibitory phenotypes in growth, development, or population reproduction of S. exigua was also shown after its larvae treated with Heliothis virescens ascovirus 3 h (HvAV-3 h) (Yu et al. 2019, 2018) or cadmium stress (Su et al. 2020). However, expressed genes of S. exigua larvae treated with EMB, HvAV-3 h or cadmium stress presented a huge difference at the transcriptome level, indicating that when coping with different stresses, even if the similar phenotypes may be shown in S. exigua, different corresponding mechanism at the molecular level may be existed.

In order to further and more accurately discover the primary cues involved in resistance to EMB, we used entire life history parameters to elucidate the physiological phenotype of S. exigua under different lethal stresses of EMB. A previous study reported by Che et al. 2015 identified these cues by only using resistance strains (Che et al. 2015). Compared with our study, the inhibitory phenotypes in growth, development, and population reproduction were not shown in the results using EMB resistant population of S. exigua by Che et al. 2015, indicating that some molecular mechanism must be regulated to help S. exigua defend against the initial stress of EMB, so as to restore its normal phenotype. Therefore, in this study, gene expression networks of S. exigua larvae treated with three different lethal concentrations of EMB was evaluated. Our results revealed that a series of important metabolic pathways (such as ABC transporters, olfactory recognition, and metabolic detoxification) may be involved in the response of the S. exigua against EMB stress and may further participate in the later evolution of resistance.

Previous studies indicated that ABC transporters play an important role in pests resistance to pesticides and other toxic compounds (Epis et al. 2014; Pignatelli et al. 2018; Zuo et al. 2018). Specifically, some studies have reported that ABC transporter genes were up-regulated in mosquito larvae after insecticide exposure (Figueira-Mansur et al. 2013; Mastrantonio et al. 2017). Also, the expression of ABCC, which belongs to ABC subfamilies, was increased in Nilaparvata lugens nymphs treated with sub-lethal concentrations of triazophos (Bao et al. 2010). In addition, another study also reported that the ABC transporter gene (ABCG3) was involved in imidacloprid detoxification in Q-type Bemisia tabaci (He et al. 2019). Also, ABC transporter genes were the important genes for S. exigua defense against insecticidal protein, such as 11 ABC transporter unigenes were reported up-regulated in Cry1Ca-treated S. exigua (Ren et al. 2020). In our study, KEGG pathway analysis of transcriptome data showed that the ABC transporter was also an important gene family involved in the response of the S. exigua to the different lethal stress of EMB. ABC transporter genes such as c15483_g1, c15483_g2 and c17259_g2 were all highly expressed under LC₅, LC₂₀ and LC₅₀ treatment. By contrast, c13729_g1 and c17287_g1 were highly expressed under LC₂₀ and LC₅₀ treatment and only c1698_g1(ABCA5) and c22475_g1 (ABCB6) were highly expressed under LC₅₀ treatment. Thus, these seven ABC transporter genes are very likely to have an important role in the evolution of resistance to EMB in the S. exigua and should be further studied.

The insect olfactory system plays an important role in directing insects towards safe food, away from toxic compounds, and locating mating partners or oviposition sites (Sato and Touhara 2009; Leal 2013; Pelosi et al. 2018; Yang et al. 2020). Silencing of the odorant receptor co-receptor (BmOrco) gene was in Bombyx mori, resulted in anti-feeding behavior and reduction in mate selection ability (Liu et al. 2017). After the olfactory-related gene Orco was disrupted, olfactory responses in both adults and larvae of Anopheles coluzzii were severely disrupted, especially for olfactory-related behaviors, such as host-seeking or oviposition in adult females (Sun et al. 2020). Further, current studies have shown that the key proteins and genes of olfactory systems have also become important molecular response targets to insecticide stress (Sato and Touhara 2009; Leal 2013; Pelosi et al. 2018; Yang et al. 2020). For example, Athetis lepigone responded to pesticides with general-odorant binding proteins (GOBPs). Previous ligand binding assays demonstrated that AlepGOBP2 in A. lepigone showed high bindings with sex pheromones, maize plant volatiles and insecticides frequently used in the field (Zhang et al. 2020). Furthermore, a study of GOBPs in Plutella xylostella indicated that they might be involved in guiding larvae to their pheromones for foraging better food (Zhu et al. 2016). In addition, two A. lepigone pheromone-binding proteins AlepPBP2 and AlepPBP3 had higher binding affinities to phoxim than to other insecticides, which indicated that they play crucial roles in the AlepPBPs/phoxim interactions (Zhang et al. 2020).

Many genes in S. exigua, such as odorant receptors (ORs), odorant binding proteins (OBPs) involved in olfactory recognition had also been identified in adults (Du et al. 2018). However, the information about olfactory recognition related gene in S. exigua larvae was still insufficient, especially the olfactory recognition related genes in response to the pesticide stress. In view of larval stage is the main chemical control period of S. exigua, identifying genes related to the special olfactory recognition of the larvae in response to pesticides was the key way to clear the resistance evolution of S. exigua against pesticides (Zhang et al. 2012; Jin et al. 2015; Xiong et al. 2019). In our study, the S. exigua treated with EMB, showed series of physiological phenotypes such as reduced egg-laying rate, refusal to feed, and growth retardation, which could be influenced by changes in the expression of olfactory genes. The calmodulin gene c21765_g1 in the olfactory transduction pathway and cAMP-dependent protein kinase catalytic subunit gene c20528_g1, were highly expressed under LC₅₀ treatment. By contrast, the calmodulin gene c4625_g2 was highly expressed under LC₅ treatment. These three olfactory genes may have important roles in the normal behavior of S. exigua, such as mating, oviposition, feeding and escaping from toxic compounds, and consequently may be special molecular targets in S. exigua against EMB stress.

Detoxification enzymes are involved in the evolution of pesticide resistance (Ketterman et al. 2011; Tang et al. 2014) . Insects usually use P450 detoxification enzymes to metabolize pesticides in their bodies (Yang et al. 2013; Scanlan et al. 2020; Wu et al. 2020). In the present study, we explored two genes involved in the P450 metabolic pathway by transcriptome sequencing. The results showed that, the gene, c13436_g1, was significantly down-regulated under LC₂₀ and LC₅₀ treatment, while gene c21930_g4 was significantly up-regulated under LC₅ treatment compared to AC samples. The significant decrease in expression of gene c13436_g1 reduced the expression of related genes in the P450 metabolic pathway. This may have resulted in the malfunctioning of this pathway which may have resulted in the death of S. exigua larvae. Moreover, the significant up-regulation of gene c21930_g4 under LC₅ treatment indicated that under low stress, a significant increase in the expression of some detoxification enzymes could overload the expression of the metabolic detoxification pathways to achieve detoxification. However, high selection pressures under LC₂₀ and LC₅₀ could cause this detoxification enzyme system to exceed its range, to cause a decreased expression of this gene under these two treatments.

In conclusion, the present study showed that key life parameters, such as larval final mortality rate, per female fecundity and intrinsic growth rate in S. exigua were affected by EMB stress in a dose-dependent manner. Different lethal doses of EMB (LC₅, LC₂₀, and LC₅₀) could all effectively inhibited the proliferation of S. exigua population by long-lasting inhibition the survival of the S. exigua larvae, and significantly reduction of per female fecundity and intrinsic growth rate of S. exigua (p < 0.05). In addition, this study also highlighted many genes responsible for EMB resistance evolution by differential transcriptome analysis. Subsequently, genes involved in the key metabolic pathways were explored, and 16 hub genes were identified as biomarkers (Additional file 10). These results help to elucidate the mechanisms underlying the resistance of S. exigua against EMB and are very significant for a comprehensive understanding of the evolution of this resistance.

Authors’ contributions

YS, HYD, LXB, YJZ, and XDK conceived and designed the experiments. YS, WZ, BJ, and SMT contributed to paper writing. YS, JZ, SMT, BSL, and HYD conducted the experiment. YS, SMT, and XDK contributed to the data analysis.

Availability of data and materials

The data are presented in the manuscript and the supporting materials. The raw reads data are submitted to the Short Read Archive (SRA) and BioProject accession number PRJNA523075 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA523075).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [grant number 31871964, 31401753, 31901889], the Special Funds for Supporting Innovation and Entrepreneurship for Returned Oversea-students in Anhui Province [grant number 2020LCX035], Anhui Provincial Natural Science Foundation [Youth project 2008085QC148], the Natural Science Foundation of Jiangsu Province [grant number BK20130717], Anhui Provincial Key Laboratory of Molecular Enzymology and Mechanism of Major Diseases, Key Laboratory of Biomedicine in Gene Diseases and Health of Anhui Higher Education Institutes.