Introduction
Insecticide resistance is now widespread and is rapidly increasing in pests across the world. This can obstruct their control in the field, therefore, threatens agriculture. Resistance is defined as the decreased susceptibility of a pest population to a pesticide that was previously effective at controlling it (Georghiou, Reference Georghiou1972; Perry et al., Reference Perry, Batterham and Daborn2011). Resistance to chemical insecticides involves two mechanisms: metabolic resistance and target-site resistance (Liu, Reference Liu2015). Metabolic resistance is manifested as the overexpression of detoxification enzyme genes. These enzymes are involved in three major metabolic detoxification gene families: cytochrome P450 (P450), glutathione S-transferases (GSTs), and carboxylesterase (CCEs). Target-site resistance is derived from the mutation or modification of genes that encoded target proteins. Traditional chemical pesticides such as dichlorodiphenyltrichloroethane (DDT) and pyrethroids are involved in sodium channels (Narahashi, Reference Narahashi1996). Moreover, there are many pathways involved in target-site insensitivity such as γ-aminobutyric acid (GABA) receptor, ryanodine receptor, and acetylcholinesterase (AChE).
The mulberry moth, Glyphodes pyloalis Walker, one of the most important mulberry pests, is widely distributed in major mulberry growing regions. G. pyloalis causes yield losses of mulberry and transmits some diseases to silkworm (Bombyx mori) (Watanabe et al., Reference Watanabe, Kurihara, Wang and Shimizu1988; Liu et al., Reference Liu, Su, Li, Li, Xu, Dai, Zhou and Wang2017). In recent years, the climate and improper use of pesticides have caused mulberry moth outbreaks, which brings serious losses to the local mulberry farmers. Thus, silk production is frequently attacked by G. pyloalis. Insecticides are the most common pest management strategy for controlling mulberry moth population in mulberry fields in China. Control of G. pyloalis has been based on the use of conventional insecticide. The abuse of insecticides in the main mulberry- growing regions of China has negative impacts on the environment. In addition, the strong selective pressure caused by conventional insecticide has led to an increase of G. pyloalis resistance to insecticides from different chemical groups, such as organophosphates, pyrethroids, and carbamates (Yazdani et al., Reference Yazdani, Sendi, Aliakbar and Senthil-Nathan2013; Piri Aliabadi et al., Reference Piri Aliabadi, Sahragard and Ghadamyari2016). However, to date, studies on the resistance of mulberry moth predominantly focused on the comparison of insecticide resistance of mulberry moth in different areas. Little is known about the relative expression of genes involved in insecticide resistance, and resistance mechanisms of mulberry moth.
The insect midgut is an important organ that plays multiple physiological function, such as food digestion and nutrient uptake (Zhang et al., Reference Zhang, Zhang, Li, Liu and Liu2016). In addition, the midgut plays a critical role in oxidative stress responses (Yanyuan et al., Reference Yanyuan, Xun, Wen, Qingjun, Shaoli, Zhaojiang, Baoyun, Xianchun, Xuguo and Youjun2014). It is also the center of insecticide detoxification and resistance development (Sparks et al., Reference Sparks, Blackburn, Kuhar and Gundersen-Rindal2013). In recent years, researchers have used new techniques and advances in genomic research to identify the mechanisms that govern insecticide resistance in pests (ME, Reference ME2008; Huang et al., Reference Huang, Shen, Chen, He, Zhou, Tan, Yuan and Yang2016). The present study adds insights to genome-wide transcriptional changes between the insecticide-susceptible and resistant G. pyloalis strains, a topic that is still rarely investigated. Genes involved in insecticide metabolic processes, including P450, GST, and CCE, were identified in the larval midgut of G. pyloalis. In addition, quantitative polymerase chain reaction (qPCR) analysis showed that CYP324A19, CYP304F14, CYP6AB10, CYP6AW1, GSTs5, CCE001c, and AChE-like were significantly induced in G. pyloalis larvae exposed to propoxur or phoxim. Thus, our findings provide an available genetic resource for further mechanistic studies of insecticide resistance of G. pyloalis.
Materials and methods
Insects
Two laboratory-selected strains of G. pyloalis were used in the experiment, an insecticide-susceptible strain and an insecticide-resistant strain. The susceptible G. pyloalis strain was reared on mulberry without insecticides treatment over the past 20 years. The resistant G. pyloalis was obtained from mulberry fields in Zhejiang (China) and was constantly exposed to propoxur and phoxim at least three times per year for about 20 years. All strains were reared at 25°C under a 12 h light-12 h dark photoperiod. Six males and six females on the third day of fifth instar larvae were randomly selected for subsequent analysis. Three biological replications were performed with each treatment.
RNA isolation and cDNA synthesis
Total RNA was isolated from the midguts of third day of fifth instar larvae using the RNAiso Plus Reagent (Takara, Kyoto, Japan). The quality and quantity of the RNA were assessed using a Nanodrop spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA). Genomic DNA was removed from total RNA by RNase-free DNase I (Takara, Kyoto, Japan). The cDNA was prepared by using the PrimeScript™ RT Reagent Kit (Takara, Kyoto, Japan) according to the manufacturer's instructions.
cDNA library construction and Illumina sequencing
The synthesized cDNA products were treated with buffer EB for end reparation and single nucleotide A (adenine) addition. Then, short fragments were connected with adapters. The suitable fragments were selected for the PCR amplification. The cDNA libraries were generated using PCR amplification. The quality of the sample library was assessed on the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and ABI StepOnePlus Real-Time PCR System. The cDNA library was sequenced using the Illumina HiSeq 4000 (New England Biolabs, Ipswich, MA, USA) platform and paired-end reads were generated.
De Novo transcriptome assembly
To obtain clean and high-quality reads for assembly, reads with only adaptors, high content (>5%) of unknown base (N), or low-quality reads were filtered. The clean reads were assembled into unigenes using the Trinity package with default parameters (v2.0.6). Trinity first configured the reads of a certain length (k-mers)into longer fragments, which are called contigs (k = 25). Next, to improve the efficiency, PCR duplication was removed by Trinity. The resulting sequences of Trinity were called transcripts. Then, gene family clustering was performed using TGICL to obtain final unigenes. Afterward, SSR detection, unigene expression analysis, and unigene functional annotation were performed. We mapped all clean reads to unigenes using HISAT (Kim et al., Reference Kim, Langmead and Salzberg2015), and then called single nucleotide polymorphisms (SNPs) with GATK (McKenna et al., Reference McKenna, Hanna, Banks, Sivachenko, Cibulskis, Kernytsky, Garimella, Altshuler, Gabriel and Daly2010). After the unreliable sites were filtered out, the final SNPs were obtained.
Functional annotation
The obtained unigenes were annotated to the following databases, NCBI non-redundant protein sequences (NR); Clusters of Orthologous Groups of proteins (KOG/COG); InterPro member database; Swiss-Prot (a manually annotated and reviewed protein sequence database); Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO). Finally, the corresponding annotation results were selected.
Differential gene expression in resistant and susceptible G. pyloalis strains
The differentially expressed genes (DEGs) between the insecticide-susceptible strain and insecticide-resistant strain of G. pyloalis were measured by calculating the number of reads per kilobase of transcript sequence per million mapped reads (RPKM). The P-value corresponded to the differential gene expression test, and the threshold P-value in multiple tests was determined by using the false discovery rate (FDR) method. We used FDR ≦0.001 and|fold change|≧2.0 as the parameters to judge the significant differences of gene expression.
Phylogenetic analysis
Multiple sequence alignments of the P450, GST and CCE amino acids were performed using MUSCLE in the MEGA6 program. A phylogenetic tree was constructed using the neighbor-joining method of the MEGA6 program. A 1000 bootstrap replication analysis was performed to evaluate the reliability of the cladogram (Tamura et al., Reference Tamura, Stecher, Peterson, Filipski and Kumar2013). The results were visualized by iTOL (http://itol.embl.de/).
Gene Ontology analysis of DEGs
We first mapped all DEGs to terms in the GO database, and then calculated gene numbers for every term. Moreover, we used a hypergeometric test to find significantly enriched GO terms in DEGs and compared them with the G. pyloalis transcriptome background. The formula for the calculation was as follows:
$$P = 1-\sum\limits_{i = 0}^{m-1} {\displaystyle{{\left( {\matrix{ M \cr i \cr}} \right)\left( {\matrix{ {N-M} \cr {n-i} \cr}} \right)} \over {\left( {\matrix{ N \cr n \cr}} \right)}}} $$Here, M indicates the number of all genes annotated to certain GO terms, and m is the number of DEGs in M. N indicates the number of all genes with GO annotation, and n is the number of DEGs in N. The calculated P-value takes the corrected P-value ≦0.05 as threshold.
Metabolic pathway analysis of DEGs
To further understand the biological function of genes, we performed metabolic pathway analysis of DEGs using Blast. KEGG is a major collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances. The significantly enriched metabolic pathways in DEGs were defined through metabolic pathway analysis. The formula was the same as that for GO analysis mentioned previously. In this equation, N is the number of all genes with KEGG annotation, and n is the number of DEGs in N. M indicates the number of all genes annotated to certain KEGG terms, and m is the number of DEGs in M.
Validation of data by quantitative real-time PCR (qRT-PCR)
To confirm the differential expression between insecticide-susceptible and insecticide-resistant strains of G. pyloalis, qRT-PCR was performed. Six DEGs, including four CYP genes and two GST genes, were randomly selected. The specific primers were designed using the Primer Premier 5 software (Premier, Canada). 1 µg RNA treated with RNase-Free DNase was used in a reverse transcription reaction. qRT-PCR was performed in a 7300 Real Time PCR System (Applied Biosystems®, Life Technologies). Ribosomal protein L32 (Rpl32) was used for internal control. The dates were analyzed with the 2−ΔΔCt method (Livak & Schmittgen, Reference Livak and Schmittgen2012). The experiment was performed three times for biological replicates.
Propoxur and phoxim treatment
To identify changes of the expression levels of detoxification enzyme genes after insecticide treatment, we selected two of the most frequently used insecticides in mulberry fields, propoxur (carbamate) and phoxim (organophosphorus). Ten milligrams of phoxim (Sigma-Aldrich, Burlington, MA, USA) and 10 mg of propoxur (Sigma-Aldrich, Burlington, MA, USA) were dissolved in 1 ml of acetone as stock solution, respectively. Then, the solutions were diluted with distilled water to 15 µg ml−1 as the working solution. The mulberry leaves were dipped in 300 ml working solution for 30 s and then air-dried. The mulberry leaves that were treated with distilled water served as control. The third day of fifth instar larvae of insecticide-susceptible G. pyloalis were fed with phoxim-treated leaves or propoxur-treated leaves for 24 h. qRT-PCR was performed as described above. Each treatment was biologically repeated for three times. The datasets were analyzed statistically using Student's t-test. The level of significance was set at P < 0.05.
Results
Transcriptome assembly and functional annotation
To obtain an overall view of the different gene expression profiles between insecticide-susceptible and insecticide-resistant strains of G. pyloalis, we performed transcriptomic analysis using larval midgut. The percentages of Q20 sequences in the two libraries were 96.28 and 96%, respectively (table S1). The data were assembled into 32,229 (table S2). For functional annotation, all assembled unigenes were matched against seven different databases (table S3). The results showed that the NR database had most of the matched unigenes.
Differential gene expression between susceptible and resistant G. pyloalis strains
Then, we compared the DEGs between resistant and susceptible G. pyloalis strains. We used FDR ≦0.001 and |fold change |≧2.0 as a threshold to filter out unigenes. Finally, we identified 1740 DEGs for insecticide-resistant strain vs. insecticide-susceptible strain (fig. 1), of which 643 were up-regulated and 1097 were down regulated in the insecticide-resistant strain.
Volcano plot of DEGs. The Y axis represents the statistical significance of the change of gene expression, and the X axis represents the fold change of gene expression in the different strains. The GP-S means the susceptible G. pyloalis and the control means the resistant G. pyloalis. The data was based on the discriminative significance values (P < 0.01) and relative expression (≥4) adopted. The red and blue points represent the significant up- and down-regulated genes, respectively.
GO analysis
To identify the putative functions of the annotated genes of the insecticide-resistant strains and the insecticide-susceptible strains, we performed GO classification and functional enrichment for DEGs. As a result, a total of 3878 unigenes were annotated into 56 GO terms using Blast2GO. Biological processes made up the majority of the GO annotations with 24 GO terms (62.5%), including 2691 unigenes. Among the biological processes, the most dominant subcategories were cellular process (63.5%) and cellular metabolic process (53.2%). The molecular function included 14 GO terms and 3300 unigenes were annotated into this category, among which the most dominant subcategories were binding (60.8%) and catalytic activity (54.7%). For the cellular component category, 1693 unigenes were annotated into 18 terms, and the most highly enriched subcategories were cell (79.7%) and cell part (79.7%). We also identified the GO functions for the DEGs and the results were shown in fig. S1. A total of 261 DEGs were classified into GO categories, including 166 up-regulated and 95 down-regulated genes.
KEGG pathway analysis
A total of 15,546 unigenes were assigned to 255 reference canonical pathways in the KEGG database (fig. S2). The results of pathway enrichment showed that there were 16 significantly enriched pathways (P < 0.05) (table S4). Obviously, ‘metabolic pathways' was the most dominant pathway in the midgut (2719 unigenes). Furthermore, nine ABC transporter DEGs were identified in the KEGG pathways. Among them, five genes (Unigene1088_All, CL1385.Contig1_All, CL1008.Contig1_All, CL138.Contig3_All, Unigene1164_All) had higher and three (CL138.Contig3_All, CL3225.Contig1_All, Unigene1213_All) had lower expression levels in the insecticide-resistant strains.
Identification of unigenes that are putatively involved in insecticide resistance
To investigated metabolic enzymes and putative insecticide target site resistance related genes, a total of 33 P450, 27 GST, and 23 CCE related unigenes were identified in the midgut of G. pyloalis (table 1). In addition, genes that were potentially involved in insecticide target-site resistance were also identified, including 1 AChE, 5 ryanodine receptors and 2 sodium channel-related unigenes (table 1).
Unigenes potentially involved in insecticide resistance.
SNPs identification in P450s of G. pyloalis midguts
The sequenced fragments of the 33 P450s from the insecticide-resistant strain were compared with those of the insecticide-susceptible strain to identify sequence variation involved in insecticide resistance. As a result, 21 SNPs were identified in three P450-related unigenes, CL1429.Contig1_All CYP9A20, Unigene2326_All CYP6AB47, and CL2956.Contig1_All CYP6AW1. The CYP9A20 gene has 15 polymorphic sites at locations 476, 572, 599, 620, 650, 850, 857, 890, 900, 919, 1082, 1142, 1253, 1311, and 1670 (table 2). Among these SNPs, only four were non-synonymous SNPs resulting from nucleotide substitutions. The substitutions led to conversion of amino acids where valine (V) changes into isoleucine (I), lysine (K) into arginine (R), asparagine (N) into lysine (K), and asparagine (N) into serine (S). CYP6AB47 has only one SNP at location 1405, which was a synonymous SNP. There were also five SNPs identified in CYP6AW1 at location 536, 1196, 1322, 1583, and 3137, among which five SNPs were all synonymous SNPs.
The identification of SNPs in P450s.
DEGs putatively related to insecticide metabolism
One of the most commonly used methods for obtaining differential expressed gene data is cluster analysis. The results of cluster analysis of the metabolic enzyme-related DEGs are shown in fig. 2. 19 DEGs were identified between susceptible and resistant G. pyloalis, including nine P450s, five GSTs, and five CCEs. To further elucidate the orthologous relationships of the DEGs that belong to the three conventional detoxification enzymes, we constructed a neighbor-joining tree of the P450s, GSTs and CCEs with a dataset of four other insect species (B. mori, Drosophila melanogaster, Cnaphalocrocis medinalis, and Papilio xuthus) and G. pyloalis.
Comparative distribution of transcripts associated with detoxification enzymes of insecticide-susceptible and resistant G. pyloalis strains. RPKM values were represented as a scale of colors ranging from green to red, which encompass the lowest (green) to highest (red) RPKM values.
To date, there are four main clans of insect P450s, CYP2, CYP3, CYP4, and the mitochondrial CYP clan, that have been identified (Scott, Reference Scott1999; Tompkins & Wallace, Reference Tompkins and Wallace2007; Feyereisen, Reference Feyereisen2011; Jean-Philippe et al., Reference Jean-Philippe, Mahmoud, Alexia and Ingraham2013). Here, nine DEGs with significant similarity with P450s were found (table 3). Phylogenetic analysis showed that the 33 midgut GpCYPs clearly segregated into four clans: CYP2, CYP3, CYP4, and Mito CYP clan (fig. 3). Out of the 33 GpCYPs, five were from the mitochondrial CYP clan, 4 from the CYP2 clan, and 4 from the CYP4 clan (fig. 3). Most of the P450 related unigenes in G. pyloalis midgut were segregated to CYP3 clan, which has 20 unigenes. Combined with the results shown in table 3 and fig. 3, among the nine DEGs, eight (CL911.Contig1_All CYP9A78, CL2253.Contig2_All CYP324A19, CL80.Contig1_All CYP6AB10, CL835.Contig2_All CYP6AB51, CL2956.Contig1_All CYP6AW1, CL129.Contig1_All CYP6AE76, CL274.Contig3_All CYP9A19, and Unigene6135_All CYP6AB49) were clustered into the CYP3 clan, and only one (CL2076.Contig2 CYP304F17) was clustered into the CYP2 clan.
The detoxification-related genes were identified from Glyphodes pyloalis.
Neighbor-joining tree of the different CYP clans of G. pyloalis, B. mori, D. melanogaster, C. medinalis and P. xuthus. The bootstrap value as percentage for 1000 replication was indicated for each branch (only values >70% are shown). The up-regulated unigenes are shown in black asterisks, while the down-regulated unigenes are shown in white asterisks.
GSTs catalyze the coordination of released GSH with xenobiotics and toxic compounds which can make them less toxic (Motoyama & Dauterman, Reference Motoyama and Dauterman1980). For insects, GSTs are mainly divided into seven classes: Delta, Epsilon, Omega, Sigma, Theta, Zeta and microsomal classes. Among these classes, the Delta and Epsilon classes are insect specific and contributed to detoxification metabolism (Chelvanayagama et al., Reference Chelvanayagama, Parker and Board2001; Ranson et al., Reference Ranson, Rossiter, Ortelli, Jensen, Wang, Roth, Collins and Hemingway2001; Yu et al., Reference Yu, Lu, Li, Fang, Zuo, Dai, Zhang and Xiang2008). As shown in fig. 4, the phylogenetic tree assigned the 27 GSTs into the six expected GST groups: Sigma, Epsilon, Theta, Omega, Zeta and Delta. The Sigma class was the dominant class, having eight GST-related unigenes. Here, we identified five DEGs that had significant similarity with GSTs. Among the five DEGs, two up-regulated unigenes were assigned to the Delta and Sigma classes (GSTd2, GSTs5), whereas the down-regulated unigenes were assigned to the Epsilon and Sigma classes (GSTe3, GSTe4, GSTs2) (table 3, fig. 4). The results suggest that these genes may participate in insecticide resistance in the midgut.
Neighbor-joining tree of the different GST classes of G. pyloalis, B. mori, D. melanogaster, C. medinalis and P. xuthus. The bootstrap value as percentage for 1000 replication was indicated for each branch (only values >70% are shown). The up-regulated unigenes are shown in black asterisks, while the down-regulated unigenes are shown in white asterisks.
CCEs constitute a large multifunctional superfamily of polyproteins and play a crucial role in insecticide resistance (Satoh & Hosokawa, Reference Satoh and Hosokawa1998). Insect CCEs can be subdivided into three major classes: intracellular catalytic, secreted catalytic, and neurodevelopmental classes (Ranson et al., Reference Ranson, Claudianos, Ortelli, Abgrall, Hemingway, Sharakhova, Unger, Collins and Feyereisen2002). As shown in fig. 5, only one unigene belongs to the neurodevelopmental class, and four to the catalytic class. Most of the CCE-related unigenes belong to the intracellular catalytic class. The intracellular catalytic class of CCEs plays a vital role in detoxifying xenobiotics, and some are involved in organophosphorus insecticide resistance (Claudianos et al., Reference Claudianos, Ranson, Johnson, Biswas, Schuler, Berenbaum, Feyereisen and Oakeshott2006). Here, we found that all the DEGs belonged to the intracellular catalytic class.
Neighbor-joining tree of the different CCE classes of G. pyloalis, B. mori, D. melanogaster, C. medinalis and P. xuthus. The bootstrap value as percentage for 1000 replication was indicated for each branch (only values >70% are shown). The up-regulated unigenes are shown in black asterisks, while the down-regulated unigenes are shown in white asterisks.
qPCR analysis
The differential expression levels between susceptible and resistant G. pyloalis strains were validated by qRT-PCR. Six DEGs (CL911.Contig1_All CYP9A78, CL2076.Contig2 CYP304F17, CL835.Contig2_All CYP6AB51, CL2956.Contig1_All CYP6AW1, Unigene16962_All GSTs5, Unigene755_All GSTd2) were selected to confirm the reliability of transcriptome data that we obtained. These unigenes had higher expression level in the insecticide-resistant larva midgut (fig. 6a). The changing trend of the qRT-PCR results were consistent with the different expression gene profiles obtained by DEG expression profiling (table S6).
(a) qRT-PCR analysis of the expression levels of six unigenes. The X axis represents the six genes selected for qRT-PCR validation, and the Y axis represents the relative expression level. (b) qRT-PCR analysis of seven unigenes after propoxur or phoxim treatment for 24 h. The X axis represents the seven genes for qRT-PCR validation, and the Y axis represents the relative expression level. Ribosomal protein L32 (Rpl32) was used as an internal control. The details of the primers are shown in table S5. Data are expressed as mean ± SE (n = 3). * P < 0.05, ** P < 0.01.
Expression of the detoxification enzyme genes with propoxur and phoxim exposure
For further exploring the expression level of the detoxification enzyme genes after propoxur or phoxim exposure, qRT-PCR was performed to identify the expression levels of the 11 DEGs, which were up-regulated in the insecticide-resistant strains. We have identified and deposited the sequences of these genes in GenBank (table S5). Furthermore, we found that seven unigenes were induced by the two insecticides. As shown in fig. 6b, the expression levels of CYP324A19, CYP304F17, CYP6AW1, CYP6AB10, GSTs5, and AChE-like were significantly up-regulated after the propoxur exposure for 24 h, while the transcripts of CYP324A19, CCE001c, and AChE-like were significantly induced in response to phoxim treatment. This result showed that these detoxification enzyme genes, especially CYP324A19 and AChE-like, were potentially involved in the detoxification of phoxim and propoxur in G. pyloalis.
Discussion
G. pyloalis is one of the most destructive pests in mulberry field that has developed a high level of insecticide resistance. In this study, we performed a transcriptome investigation using the midguts of the insecticide-susceptible and insecticide-resistant G. pyloalis larvae. A total of 1740 DEGs were identified, including 643 up- and 1097 down-regulated genes. The abundance of detoxification metabolism-related genes indicated that the insecticide metabolic resistance level might increase in the insecticide-resistant G. pyloalis strain. However, most of these unigenes are not directly associated with the molecular mechanisms of insecticides resistance. The results are consistent with those of insecticide-susceptible and resistant strains of Spodoptera frugiperda (Lepidoptera: Noctuidae) (Nascimento et al., Reference Nascimento, Fresia, Cônsoli and Omoto2015) and Plutella xylostella (Lepidoptera: Plutellidae) (Lin et al., Reference Lin, Jin, Hu, Chen, Yin, Li, Dong, Zhang, Ren and Feng2013). It implies that studies on insecticide resistance mechanisms are still insufficient, and many genes related to insecticide resistance have not yet been discovered. To the best of our knowledge, our study is the first transcriptomic analysis of G. pyloalis involved in insecticide resistance and the results offer novel insights into transcriptomics adaptations to insecticide in G. pyloalis.
Insecticide resistance occurs in many insects. It is an adaptive process that results from the selection pressure from the continuous exposure to insecticides (Ffrench-Constant, Reference Ffrench-Constant2007). Many studies on the resistance of chemical insecticides focused on three main detoxification enzymes: CYPs, GSTs, and CCEs. These enzymes are responsible for insect metabolic resistance to insecticides. Our differential expression results between the susceptible and resistant G. pyloalis strains have found unigenes that are homologous to CYPs, GSTs, and CCEs. Here, we identified that 33 unigenes had significant similarity with P450s. Among these 33 P450s, nine were DEGs, of which six were up-regulated and three were down-regulated. Five up-regulated genes were assigned to the CYP3 clan. In general, the clan 3 P450 genes in insects can be further subdivided into CYP6, CYP9 and CYP321 families, which encode enzymes involved in xenobiotic detoxification (Pittendrigh et al., Reference Pittendrigh, Aronstein, Zinkovsky, Andreev, Campbell, Daly, Trowell and Ffrench-Constant1997; Yang et al., Reference Yang, Chen, Wu, Yue and Wu2006; Li et al., Reference Li, Schuler and Berenbaum2007). Earlier reports have identified that the CYP6 subfamily is linked to the detoxification of deltamethrin resistance in Tribolium castaneum (Zhu et al., Reference Zhu, Parthasarathy, Bai, Woithe, Kaussmann, Nauen, Harrison and Palli2010). In other Lepidoptera, Amyelois transitella CYP6AB11 enzymes have been identified and implicated in detoxification of several insecticides (Niu et al., Reference Niu, Rupasinghe, Zangerl, Siegel, Schuler and Berenbaum2011). Here, CYP6AB51 and CYP6AB10 were found to be overexpressed in the insecticide-resistant G. pyloalis strain. In addition, enzymes from the CYP9A subfamily are associated with resistance to pyrethroid and midacloprid in Helicoverpa armigera (Yang et al., Reference Yang, Chen, Wu, Yue and Wu2006) (CYP9A12A) and B. mori (Yamamoto et al., Reference Yamamoto, Ichinose, Aso and Fujii2010) (CYP9A20). Similarly, we identified that CYP9A78, clustered with the CYP9A subfamily, was significantly overexpressed in the insecticide-resistant G. pyloalis strain. In addition, a recent report showed that CYP304F-1 had higher expression levels in insecticide-resistant P. xylostella strains (Yu et al., Reference Yu, Tang, He, Ma, Vasseur, Baxter, Yang, Huang, Song and You2015). Likewise, we found that CYP304F17 was highly expressed in the insecticide-resistant strain. It implies that the five overexpressed CYP3 clan genes might involve in the insecticide resistance of G. pyloalis. Furthermore, the expression levels of four P450 genes (CYP324A19, CYP304F17, CYP6AW1, and CYP6AB10) were significantly increased after propoxur treatment for 24 h, while CYP324A19 was significantly induced with phoxim exposure for 24 h. The results showed that the resistance-related genes have a different response to different insecticides, suggesting that these P450 genes probably have their own specific substrates. Further genotype-phenotype association studies using field populations are needed to identify the functional polymorphism responsible for the increased expression of these P450 genes.
P450s are dominant regulators of insecticide metabolism. Certain SNPs have been shown to have a tremendous effect on P450 activity (Preissner et al., Reference Preissner, Hoffmann, Preissner, Dunkel, Gewiess and Preissner2013). To identify sequence variations involved in insecticide resistance, we performed SNP analysis of G. pyloalis midgut P450s. Twenty-one SNPs were found in three P450-related unigene sequences, which were due to nucleotide substitution. A previous study identified 12 SNPs, eight of which were found in CYP6D1 in a pyrethroid-resistant housefly strain (Mahmood et al., Reference Mahmood, Højland, Asp and Kristensen2016). A recent report identified SNPs in CYP6P9 and CYP6P4 in Anopheles funestus and validated them to be resistance markers in laboratory pyrethroid-resistant strains (Wondji et al., Reference Wondji, Irving, Morgan, Lobo, Collins, Hunt, Coetzee, Hemingway and Ranson2009). Similarly, our findings showed that CYP9A20, CYP6AB47, and CYP6AW1 had 15, 1, and 5 SNPs, respectively. Earlier reports identified BmCYP9A20 as being highly overexpressed in B. mori after 24 h phoxim exposure (Gu et al., Reference Gu, Zhou, Xie, Li, Ma, Sun, Wu, Wang, Wang, Hong, Shen and Li2014). Our results showed that the 15 SNPs in CYP9A20 were non-synonymous, which led to conversion of amino acids. CYP9A20 was more metabolically active against insecticides compared with other genes. Further analysis of these SNPs associated with insecticide resistance will be of great interest.
GSTs play a critical role in the detoxification metabolism not only of endogenous xenobiotics such as insecticides but also of endogenous compounds such as lipid peroxidation (Enayati et al., Reference Enayati, Ranson and Hemingway2005; Huang et al., Reference Huang, Huang, Sun, Huang and Wang2017). Here, we have identified 27 GSTs from G. pyloalis larval midguts, with five DEGs in the insecticide-resistant strains. The overexpression of GSTs in the insecticide-resistant G. pyloalis strain indicates their involvement in detoxification of insecticides. Several reports have showed that insect GSTs might be responsible for resistance to organophosphate (Huang et al., Reference Huang, Hu, Yao, Wu, Chiang and Sun1998) and pyrethroids (Vontas et al., Reference Vontas, Small and Hemingway2001, Reference Vontas, Small, Nikou, Ranson and Hemingway2002; Enayati et al., Reference Enayati, Ranson and Hemingway2005). Wu et al., showed that GST activities were significantly higher in phoxim-resistant strain compared with the susceptible strain in H. armigera (Wu et al., Reference Wu, Liang and Guo1997). DmGSTD1 is overexpressed in DDT-resistant D. melanogaster (Tang & Tu, Reference Tang and Tu1994). Similarly, high expression levels of GSTs5 and GSTd2 were observed in insecticide-resistant G. pyloalis. Moreover, qRT-PCR analysis showed that the expression level of GSTs5 (CL612.Contig2_All) was significantly up-regulated after propoxur exposure for 24 h, suggesting that it may be involved in insecticide resistance.
For insects, another important metabolic enzyme family involved in insecticide resistance is CCE. Several studies on insect CCEs have shown that the enzymes may be involved in the detoxification of organophosphates, pyrethroids, and carbamates through gene upregulation, amplification, or mutation (Field & Devonshire, Reference Field and Devonshire1998; Claudianos et al., Reference Claudianos, Russell and Oakeshott1999). Previous studies have described that the resistance mechanisms of CCEs to insecticide are associated with specific site-mutations. These mutations result in the overexpression of CCEs, which decreases insecticide metabolism (Cui et al., Reference Cui, Qu, Cong, Liu and Qiao2007). Notably, AChE which belongs to the neurodevelopmental class, J clade CCEs, is the target site of OP and carbamate insecticides (Oakeshott et al., Reference Oakeshott, Claudianos, Campbell, Newcomb and Russell2005). The expression levels of four DEGs were higher in the insecticide-resistant than in the insecticide-susceptible G. pyloalis strain. Furthermore, we showed that CCE001C and AChE-like were significantly induced in G. pyloalis larvae exposed to insecticide. Taken together, we have described the important relationships between the three detoxification enzymes and metabolic process of insecticides in G. pyloalis. However, the mechanism of insecticide resistance in insects is a sequence of complex metabolic processes. The ABC transporter is a solute transport system that plays a role in the interaction between organisms and the environment (Lei et al., Reference Lei, Zhu, Xie, Wu, Wang, Guo, Xu, Li, Zhou and Zhang2014). To date, the involvement of ABC transporters in insecticide resistance for arthropods has not been well documented. However, an increasing number of reports have showed the important roles of ABC transporters in insecticide resistance (Dermauw & Van Leeuwen, Reference Dermauw and Van Leeuwen2014). He et al. (Reference He, You, Vasseur, Yang, Xie, Cui, Bai, Liu, Li and Xu2012), have identified ABC transporters as important detoxification enzymes in P. xylostella. Epis et al. (Reference Epis, Porretta, Mastrantonio, Urbanelli, Sassera, Marco, Mereghetti, Montagna, Ricci and Favia2014), found that AnstABCG4 expression is up-regulated in permethrin-treated Anopheles stephensi. Here, we identified the up-regulated expression of ABC transporters in insecticide-resistant G. pyloalis. Likewise, a recent study revealed that the transcription of CmABCG1, CmABCG4 and CmABCC2 were significantly induced by chlorpyrifos in C. medinalis (Yu et al., Reference Yu, Xu, Wang, Ma, Yu, Fei, Zhang and Wang2017). In this study, the identification of transcripts in these classes of ABC transporters promotes further studies into their involvement in insecticide resistance.
Conclusion
In summary, we present the first whole-transcriptome analysis for differential expression genes in insecticide-resistant and susceptible G. pyloalis strains. Genes involved in insecticide metabolic processes, including P450, GST, and CCE, were identified in the larval midguts. Seven of these detoxification enzyme genes were significantly induced by two insecticides. In addition, 15 SNPs in CYP9A20 were non-synonymous and led to the conversion of amino acids in different locations. Overall, the transcriptome results revealed several genes that are potentially relevant to insecticide-resistant of G. pyloalis. Many of which need be understood with deep functional studies. Nonetheless, our study presents a new layer of complexity in studies of insecticide resistance and facilitates the development of insecticides with greater specificity to G. pyloalis.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485319000257.
Acknowledgments
We thank Pengfei Zhan and Qiujie Qian for technical assistance. This work was funded by the National Natural Science Foundation of China (Grant No. 31572321; 31602010). The work was also supported by the public project grant (2016C32005) from the Science and Technology Department of Zhejiang Province and the grant (No.2015GZ06) from the Science and Technology Department of Huzhou.
Conflict of interest
The authors declare no competing financial interest.
