Transcriptome Analysis and Identification of Major Detoxification Gene Families and Insecticide Targets in Grapholita Molesta (Busck) (Lepidoptera: Tortricidae).

The oriental fruit moth, Grapholita molesta (Busck) (Lepidoptera: Tortricidae), is an important pest of most stone and pome fruits and causes serious damage to the fruit industry worldwide. This insect pest has been primarily controlled through the application of insecticides; as a result, G. molesta has developed resistance to many different types of insecticides. To identify detoxification genes, we have, de novo, sequenced the transcriptome of G. molesta (Lepidoptera: Tortricidae) and yielded 58,970 unigenes of which 26,985 unigenes matched to known proteins. In total, 2,040 simple sequence repeats have been identified. The comprehensive transcriptome data set has permitted us to identify members of important gene families related to detoxification in G. molesta, including 77 unigenes of putative cytochrome P450s, 28 of glutathione S-transferases, 46 of Carboxylesterases, and 31 of insecticide targets. Orthologs of some of these unigenes have shown to play a pivotal role in insecticide resistance in other insect species and those unigenes likely have similar functions in G. molesta.

The oriental fruit moth, Grapholita molesta (Busck) (Lepidoptera: Tortricidae), is an important pest of most stone and pome fruits primarily from the family Rosaceae (Rothschild and Vickers 1991, Kirk et al. 2013). G. molesta is invasive and capable of adapting to various environmental conditions (Hughes et al. 2004). It attacks both vegetative tissues and fruits and can cause significant loss to the fruit industry (Du et al. 2015). Adults lay eggs adjacent to young shoots and successive larval instars continue to feed on the shoots causing the dieback or distortion of the growing tips. Serious damage occurs when their larvae bore into the fruit and feed around the seed. The final-instar larvae leave the feeding plant materials to search for split subcutaneous site for pupation (Rothschild and Vickers 1991). Historically, this insect pest is primarily controlled by the frequent application of insecticides which has led to insecticide resistance in its field populations (Kanga et al. 2003). The insecticide resistance in G. molesta was first reported in 1960. Since then, it has developed resistance to several classes of insecticides, such as azinphosmethyl and phosmet (organophosphorus), methomyl (carbamates), endosulfan (organochloride), and cypermethrin (pyrethroids) (Pree et al. 1998; Kanga et al. 1999, 2003) (www.pesticideresistance.com).

In insect, there are two major biochemical resistance mechanisms: increased metabolism and target site insensitivity to insecticides. The increased metabolism of insecticides is mainly caused by the elevated activities of detoxifying enzymes, such as cytochrome P450s (P450s), glutathione S-transferases (GSTs) and carboxylesterases (CarEs) (Ranson et al. 2002, Enayati et al. 2005, Feyereisen 2005, Oakeshott et al. 2005). There are several excellent reviews of these enzymes in insecticide resistance (McKenzie and Batterham 1994, Scott et al. 1998, Ranson et al. 2002, Daborn and Le Goff 2004, Feyereisen 2005, Oakeshott et al. 2005, Li et al. 2007, Yan et al. 2009, Liu 2015). The cytochrome P450 genes distribute in 265 different families and subfamilies (Nelson 2009) and the number of these genes is growing rapidly (Feyereisen 2006). Majority of P450 genes linked to insecticide resistance belong to CYP4, 6, 9, and 12 families (Feyereisen 2005, Li et al. 2007, Guo et al. 2012a,b). GSTs are important phase II metabolic enzymes and widespread in most organisms (Hayes and Pulford 1995). Insect GSTs can be divided into seven classes: Delta, Epsilon, Omega, Sigma, Theta, Zeta, and Microsomal. The Delta and Epsilon classes are frequently associated with the insecticide resistance (Alias and Clark 2007, Lumjuan et al. 2011, Riveron et al. 2014). The third major family of the detoxifying enzymes is CarEs. The CarE genes are frequently implicated in the insect resistance to organophosphates, carbamates, and pyrethroids through gene amplification, upregulation, coding sequence mutations, or a combination of these mechanisms (Li et al. 2007, Zhang et al. 2011).

The major targets of the insecticides include acetylcholinesterase (AchE), acetyl-CoA carboxylase (ACCase), nicotinic acetylcholine receptor (nACHR), voltage-gated sodium channel (VGSC), γ-aminobutyric acid receptor, glutamate-gated chloride channel (GluCl), and ryanodine receptor (RyR) (Wei et al. 2013, Sparks and Nauen 2015). In insect, acetylcholine (ACh) esterase has shown to be the target of two types of insecticides organophosphates and carbamates, while the nACHR is the target of neonicotinoids, GABA receptors/GluCl of organochlorides and fipronil, and VGSC of pyrethroids (Ichikawa 2015). In addtition, sodium and GABA-activated channels are the targets of pyrethroids and cyclodienes. When exposed to pyrethroids and 2,2-bis(4-Chlorophenyl)-1,1,1-trichloroethane (DDT), mammals, and insects show severe neurological symptoms with modified sodium channel activity. These insecticides have also been indicated to suppress GABA and glutamate receptor-channel complexes (Narahashi 1992, 1996). Furthermore, ligand-gated chloride channels have been shown to act as the target sites for several groups of insecticides (Bloomquist 2003, Insecticide Resistance Action Committee).

Till now, very little study has been done in G. molesta with regards to the gene families associated with insecticide resistance. Despite the economic importance of G. molesta, at the time of this research, only a limited number of studies had been conducted at the molecular level, including the sequencing of the mitochondrial genome (Gong et al. 2012), characterization of two general odorant binding protein genes (Zhang et al. 2012), population genetic structure (Kirk et al. 2013), transcriptome analysis of sex pheromone glands (Jung and Kim 2014), and characterization of two heat shock proteins (Chen et al. 2014). Although these studies generated rich genomic and EST sequence data for G. molesta, it lacked sequence information related to detoxification or other insecticide related genes. To identify genes of P450s, GSTs, CarEs, and insecticide targets in G. molesta, we have first de novo sequenced the transcriptome of G. molesta. The comprehensive transcriptome data set has permitted the identification and phylogenetic classification of important gene families that are potentially associated with insecticide resistance.

Materials and Methods

Insect Samples

A laboratory colony of G. molesta was derived from a field collection in 2012 from a peach orchard in Taigu, Shanxi Province, China, and had been reared indoors with artificial diet without exposure to any insecticides for four years. The insect was maintained in a growth chamber with 26 ± 1 °C, 70–80% RH, and a photoperiod of 15:9 (L:D) h (Wang et al. 2011). To generate a comprehensive transcriptome dataset, insect samples were collected from various developmental stages of G. molesta. To collect the eggs, male and female adults of G. molesta were placed in a glass tank. The bottom of the tank was covered by a waxed paper, and the top opening of the tank was covered with a cotton ball spread with honey water to supply nutrients for the moths. Once the females deposited their eggs onto the waxed paper, the eggs were then collected. To collect the samples of other developmental stages, the neonates were transferred into individual tubes immediately after hatching and reared on the artificial diet. When larvae pupated, the pupae were transferred into individual pupation boxes. Unmated male and female moths were collected once individual insect had emerged. Due to the difference in the insect size from the developmental stages, different numbers of the insects were chosen. The final samples included 1,000 eggs, 500 first instar larvae, 100 second instar larvae, 40 third instar larvae, 40 pupae, and 40 adults (20 males and 20 females), and all samples were kept in liquid nitrogen until further process.

RNA Isolation, cDNA Library Construction, and Illumina Sequencing

Total RNA was extracted separately from insect samples of each developmental stage using Trizol reagent (Invitrogen, Carlsbad, CA) with RNase-free DNase I (Takara, China) treatment following the manufacturer's instructions, and then purified using RNeasy spin columns (Qiagen, Valencia, CA, USA). The RNA was quantified using a NanoVue UV-Vis spectrophotometer (GE Healthcare Bio-Science, Uppsala, Sweden) and examined for its integrity by 1% agarose gel electrophoresis.

RNA samples with equal amount from each developmental stage were combined and used for the construction of the cDNA library and sequencing following the standard protocols from the BGI (Shenzhen, China). Briefly, cDNAs were synthesized from the mRNA fragments with random oligonucleotides primers using RNase H (Invitrogen, Carlsbad, CA) and DNA polymerase I (New England BioLabs). The cDNA fragments were then ligated to the adapter containing primer sequences, amplified through PCR and sequenced using paired-end Illumina HiSeq 2000 technology (Illumina Inc., San Diego, CA, USA).

De Novo Assembly and Bioinformatics Analysis

Before sequence analysis, sequencing adaptors and low-quality (defined as with 5% or more Ns) sequence reads were removed. The clean reads were de novo assembled using Trinity (Grabherr et al. 2011). Unigenes <200 bp were excluded from further analysis. Assembled unigenes (>200 bp) were divided into either clusters or singletons. The clusters were prefixed with “CL” by the cluster id and contig number. One cluster includes more than one contig with sequence similarity >70%. Singletons were prefixed with “Unigene.”

Functional Annotation of Unigenes and Marker Discovery

To annotate the unigenes, distinct sequences were first searched against the National Center for Biotechnology Information (NCBI). Unigene sequences were searched by BLASTX to NCBI nucleotide databases, and aligned with Swiss-Prot knowledgebase (UniProtKB)/Swiss-Prot, KEGG (http://www.genome.ad.jp/kegg/) and COG databases (http://www.ncbi.nlm.nih.gov/COG/) (with a cutoff e-value≤ 1e − 10). Unigenes were accepted if more than 60% of the sequences (not including Ns) matched to known genes. Individual unigenes if matched to more than five different known genes have been excluded according to the criteria established by Birzele et al. (2010). Sequences with BLASTX hits were annotated according to GO terms using Blast2GO (http://www.blast2GO.com) and displayed by WEGO (http://wego.genomics.org.cn/cgi-bin/wego/index.pl). Unigenes that had no match in the databases were further analyzed using the ESTScan software to identify potential open reading frames (ORFs), translation directions, and putative amino acid sequences (Iseli et al. 1999). Potential simple sequence repeat (SSR) markers were identified with the software MicroSAtellite (http://pgrc.ipk-gatersleben.de/misa/) (Sharma et al. 2007) with unigenes as references. All these data were submitted to NCBI by Bankit.

Identification of Major Detoxification Gene Families

To identify the unigene sequences encoding P450s, GSTs, CarEs, as well the target proteins of six different classes of commonly used insecticides, the gene annotation list of the G. molesta Nr database was searched using the full names, abbreviated names and alternative names of these gene families. The candidate unigenes were blasted (BLASTX) against NCBI database to further validate annotation predictions with cutoff E-value < 1.0E  10. The prediction of the ORFs) and amino acid sequence translation were carried out using ExPaSy (http://www.expasy.org/tools/dna.html). The major detoxification gene families in G. molesta were compared to different model insect species, Bombyx mori, Anopheles gambiae, Drosophila melanogaster, and Tribolium castaneum. Phylogenetic trees of P450s, GSTs, and CarEs gene families were constructed based on neighbor-joining using MEGA 5 (Tamura et al. 2011).

Results and Discussion

Transcriptome Sequencing and De Novo Assembly

A comprehensive transcriptome dataset (GenBank accession numbers: SRP068815: PRJNA309585) was generated by Illumina sequencing of combined different developmental stages of G. molesta. The statistics of the sequencing are summarized in Table 1. The sequencing of the G. molesta transcriptome yielded a total of 56,506,100 raw reads and 4.69 Gb nucleotides. In total, 52,162,048 (92%) clean reads were assembled into 125,246 contigs and 58,970 unigenes that consist of 17,620 clusters and 41,350 singletons. Of these unigenes, 9,012 (15.3%) were longer than 1,000 bp (but shorter than 2,000 bp), 2,376 (4.0%) were between 2,000 and 3,000 bp, and 1,275 (2.2%) were longer than 3,000 bp. The size distributions of all the contigs and unigenes are shown in Supp Figure 1 [online only]. The number of unigenes in G. molesta is comparable with other insect species, such as Liposcelis entomophila (54,220 unigenes) (Wei et al. 2013), Bactrocera dorsalis (49,804 unigenes) (Shen et al. 2011). Up to today, the only published transcriptome data set from G. molesta was generated from female antennae with 5.2 million clean reads (one tenth of the read number in this present article) (Li et al. 2015). The significantly large sequence data set generated from the various stages of G. molesta in our study provides valuable resources for the scientific community in the molecular and genetic studies of many aspects of insect biology not only in this important insect pest but also in other related insect pests.

Table 1.

Summary statistics for the Illumina sequencing from Grapholita molesta transcriptome

Sequencing
Total number of raw reads	56,506,100
Total number of clean reads	52,162,048
Total clean nucleotides (bp)	4,694,584,320
Q20 percentage (%)	97.99
N percentage (%)	0.00
Number of contigs	125,246
Mean length of contigs (bp)	346
Number of unigenes	58,970
Mean length of unigenes (bp)	738
Unigenes annotations against nr	26,985
Unigenes annotations against KEGG	17,779
Unigenes annotations against COG	9,328
Unigenes annotations against GO	12,188

SSR Discovery

As currently there is no genome sequence accessible for G. molesta, the transcriptome provided an excellent opportunity to identify the SSRs distribution in a high throughput manner for this insect species. To evaluate the quality of assembly and facilitate genetic studies, the distribution of SSRs in G. molesta transcriptome was analyzed. In total, 2,040 SSRs were identified. The distribution of SSRs is shown in Supp Figure 2 [online only]. Fifty-one motifs of SSRs were found in the transcriptome of G. molesta. The main repeats of the identified SSR motif frequencies were A (16.37%), C (14.17%), AT (3.87%), GC (2.79%), TA (2.75%), CGC (2.70%), GAT (2.30), GCC (2.79%), GCG (2.84%), and GGC (3.73%). Using de novo transcriptome sequencing approach to discover SSRs has been performed in many other insect species, such as Diaphorina citri (Chen et al. 2014), L. entomophila (Wei et al. 2013), and Maruca vitrata (Margam et al. 2011).

SSRs have also been widely used in the construction of linkage maps, map-based cloning, and marker-assisted selection and genome-wide association studies in many insect species (Cao et al. 2015, Schemerhorn et al. 2015, Wu et al. 2016). The SSRs identified in this research will be valuable information for the genetic studies of insecticide-related issues and many other aspects in G. molesta.

Annotation of Unigenes

In total, 26,985 unigenes (45.8% of all the unigenes) matched to genes with known functional proteins. The distribution of E-values, percentages of unigenes with matches, and matched insect species are shown in Figure 1. The percentage of unigenes matching known functional proteins was equivalent or much higher than those reported for D. citri (44.8%) (Chen et al. 2014), Nilaparvata lugens (40.0%) (Bao et al. 2012), Bemisia tabaci (16.2%) (Wang et al. 2010), and Plutella xylostella (22.3%) (He et al. 2012), suggesting the high quality of the de novo assembly of G. molesta transcriptome sequences. The advancement of overall gene annotation in the public domains should help improve the gene annotation for G. molesta and other insect species. Among the 58,970 unigenes, 31,985 (54.2%) showed no significant similarities to known function genes (E-value > 10 − 10) and some of these unigenes may be unique to G. molesta or not characterized yet.

Fig. 1.

Homology analysis of unigenes for Grapholita molesta. (A) E-value distribution of BLAST hits for each unigenes with a cut-off E-value of 1.0E −5. (B) Similarity distribution of the top BLAST hits for each sequence. (C) Species distribution.

For the predicted proteins of the 26,985 unigenes, 16,155 (59.9%) proteins had over 60% similarities with other insects in the NR database (Fig. 1C). The majority of these sequences (66.0%) matched to sequences from the monarch butterfly Danaus plexippus while a small portion of the sequences matched to other insect species, such as B. mori, the mode of Lepidoptera. Although they all belong to the Lepidoptera insects, B. mori for genomic sequencing has been domesticated >10 thousand years without exposure to any insecticides (Xia et al 2009), whereas D. plexippus was captured in the wild were used as starting material for genomic sequencing (Zhang et al. 2011). The moth for sequencing in our research has been captured in 2012 and reared in the laboratory with artificial diet without exposure to any insecticides.

Further Functional Annotation and Metabolic Pathway of Unigenes

To further annotate the function of the unigenes, gene ontology (GO) assignment and cluster of orthologous groups (COG) of proteins were first used to categorize the functional groups of the predicted proteins. In total 12,188 (45.17%), unigenes were annotated against GO. These assigned unigenes were distributed into 58 functional groups of three ontologies; biological process was the largest category (55.59%) followed by cellular component (28.66%) and the molecular function (15.75%) (Fig. 2; Supp Table 1 [online only]). The unigene distribution patterns in the three categories were similar to those found in Adelphocoris suturalis (Tian et al. 2015), Bactrocera oleae (Pavlidi et al. 2013), Tomicus yunnanensis (Zhu et al. 2012), but differs from those of L. entomophila (Enderlein) (Wei et al. 2013), Cimex lectularius (Shen et al. 2011), and B. dorsalis (Bai et al. 2011). Within the biological process category, the metabolic process was among the predominant groups along with cellular process, signaling and single-organism process.

Fig. 2.

Classification of the gene ontology (GO) for the transcriptome of Grapholita molesta.

The COG of proteins has also been used to predict and classify the possible functions of the unigenes. Of these 25 COG categories, “General function prediction only” represents the most common category, followed by “translation, ribosomal structure and biogenesis,” and “replication, recombination and repair.” By contrast, the “RNA processing and modification,” “Extracellular structures” and “Nuclear structure” fell into the smallest groups (Fig. 3). The lack of genomic information in the genus Grapholita may contribute to the relatively low number of unigenes having COG founctional classification. However, the results were similar to the transcriptome studies of B. dorsalis (Bai et al. 2011), A. suturalis (Tian et al. 2015), and T. yunnanensis (Zhu et al. 2012).

Fig. 3.

COG classification of the Grapholita molesta transcriptome.

KEGG can be used to analyze biochemical pathways including metabolic pathways and regulatory pathways. In total, 17,779 unigenes were mapped to 258 KEGG pathways (Supp Fig. 3 [online only]). The top one pathway is “metabolic pathways,” which include 14.69% unigenes. The similar patterns were also found in the transcriptome studies of B. dorsalis (Bai et al. 2011) and T. yunnanensis (Zhu et al. 2012). The high representation of unigenes involved in the metabolic process identified through GO, COG, or KEGG analysis provides valuable resources for the biological and molecular study of G. molesta, particularly in relation to detoxification.

Identification of the Major Detoxification Gene Families

To identify unigenes associated with insecticide resistance, we have focused on the three major detoxification enzyme gene families P450s, GSTs, and CarEs, as well unigenes encoding insecticide target proteins. The unigenes of G. molesta which fell into these gene families were identified and further classified in comparison with four model insect species, B. mori, A. gambiae, D. melanogaster, and T. castaneum (Table 2).

Table 2.

Classification of the cytochrome P450s, GSTs) and CarE (CarEs) in different insect species

Enzymes/class	G. molesta	B. mori	A. gambiae	D. melanogaster	T. castaneum
P450s
CYP2	14	7	10	7	8
CYP3	30	28	40	36	79
CYP4	18	33	46	32	47
Mitochondrial	15	10	9	12	9
Subtotal	77	78	105	87	143
GSTs
Delta	5	4	12	11	3
Epsilon	9	8	8	14	19
Omega	4	4	1	5	4
Sigma	4	2	1	1	7
Theta	1	1	2	4	1
Zeta	2	2	1	2	1
Microsomal	3	0	3	1	0
Other	0	2	3	0	0
Subtotal	28	23	31	38	35
CarEs
Dietary class
A clade	7	42	0	0	0
B clade	0	13	14	2	14
C clade	3	0	2	11	12
Hormones and pheromone
D clade	0	4	0	3	2
E clade	9	2	4	3	7
F clade	0	0	4	2	2
G clade	0	2	4	0	0
H clade	1	1	9	2	1
I clade	1	0	2	3	1
Neurodevelopmental
J clade	2	2	2	0	2
K clade	0	1	1	0	1
L clade	0	1	5	0	5
M clade	2	6	2	0	2
N clade	1	2	0	0	0
Other	20
Subtotal	46	76	49	26	49

Cytochrome P450s

In total, 138 unigenes encoding putative cytochrome P450s (P450s) were identified from the G. molesta transcriptome dataset (Supp Table 2 [online only]). For the phylogenic analysis, we removed the unigene sequences <500 bp, which are <150 residues after protein translation. Such short amino acids were not sufficient for the phylogenic analysis. After removing the short unigene sequences, 77 P450 unigenes were identified in the G. molesta. The number is comparable to B. mori (78) (Xia et al. 2004) and D. melanogaster (87) (Tijet et al. 2001), but lower than T. castaneum (143) (Zhu et al. 2013) and A. gambiae (105) (Pavlidi et al. 2013) (Table 2). Phylogenic analysis revealed the distribution of the 77 P450 unigenes in G. molesta that was 30 in CYP clade 3, 18 in CYP clade 4, 14 in CYP clade 2, and 15 in the Mitochondrial clade, the four clades of P450 genes in insect (Fig. 4) (Nelson 1998, Feyereisen 2006). Like other species, it appeared that G. molesta has relatively more members in CYP3 and CYP4 clades (Table 2). CYP3 clade is a large group of insect CYPs and comprises CYP6, CYP9, CYP28, and CYP308-310 families. CYP4 clade includes >20 families, but mainly CYP4 family. The genes of these P450 families participate primarily in xenobiotics metabolism and insecticide resistance (Feyereisen 2006; Guo et al. 2012a,b). It is possible that during the evolution, insects have developed detoxification genes to adapt to the toxins present in their diets and defense themselves from toxins in their environment, such as insecticides (Schuler 2011).

Fig. 4.

Phylogenetic analysis of cytochome P450s in Grapholita molesta (Gm) with other insects. Sequences were aligned using ClustalW. Nodes with >50% bootstrap support (1,000 pseudoreplicates) are indicated. Bm, Bombyx mori; Dm, Drosophila melanogaster; Ms, Manduca sexta; Ha, Helicovera armigera; Sl, Spodoptera littoralis; Tc, Tribolium castaneum.

Twenty-four P450 unigenes were mapped in the KEGG pathways, such as “metabolic pathways” (Ko01100), “metabolism of xenbiotics by cytochrome P450” (Ko00980), “insect hormone biosynthesis” (Ko00981), “Drug metabolism-cytochrome P450” (Ko00982), and “Drug metabolism-other enzymes” (Ko00983) (Supp Table 3 [online only]). Sixteen P450 unigenes, which were potentially involved in metabolism related pathways, were further analyzed. Majority of them were grouped into CYP6 families and two into CYP9 family (Table 3). The CYP6 family is unique to the class Insecta, and its biochemical function has essentially been associated with the metabolism of xenobiotics (Bergé et al. 1998). For example, in D. melanogaster, knockdown of CYP4g1, CYP6g1 and CYP12d1 by UAS-RNAi system increased its susceptibility to DDT (Gellatly et al. 2015, Sun et al. 2015); in Aedes aegypti, several CYP6 and CYP9 genes showed to degrade deltamethrin in CAL pyrethroid-resistant strains (Dusfour et al. 2015); in Anopheles funestus the CYP6Z1 gene conferred cross-resistance to carbamate/pyrethroid (Ibrahim et al. 2016). Moreover, CYP6B subfamily could be induced by certain monoterpenses and high concentrations monoterpenses caused pyrethroid resistance in Helicoverpa armigera (Liu et al. 2013). Overexpression of CYP6BQ9, a brain-specific cytochrome P450 gene, resulted in the increased the metabolism of deltamethrin in T. castaneum (Zhu et al. 2010). Several CYP6F subfamily genes and CYP9AQ2 in Locusta migratoria were associated with the detoxification of deltamethrin or carbaryl (Guo et al. 2015, 2016). The sixteen P450s potentially play similar roles in insecticide detoxification.

Table 3.

Putative cytochrome P450s revealed from Grapholita molesta transcriptome database

Unigene ID	Accession No.	Length (bp)	Identity (%)	E-value	Similar subfamily by BlastX
CL9294.Contig1	KY421206	1709	64	0.0	CYP6AB (D. pastinacella)
Unigene34892	KY421207	1767	55	0.0	CYP6AB (B. mandarina)
Unigene35102	KY421208	1957	69	0.0	CYP6AB (S. littoralis)
Unigene35103	KY421209	1930	69	0.0	CYP6AB (C. suppressalis)
CL4219	KY421210	1914	53	0.0	CYP6AE (H. armigera)
Unigene31911	KY421211	1210	51	4e−134	CYP6AE (M. sexta)
CL8871.Contig1	KY421212	1658	57	0.0	CYP6B (P. machaon)
Unigene2501	KY421213	1566	55	0.0	CYP6B (M. sexta)
Unigene34714	KY421214	1829	58	0.0	CYP6B (P. machaon)
Unigene25288	KY421215	1042	59	6e−149	CYP6B (A. transitella)
CL8472.Contig1	KY421216	1621	37	4e−103	CYP6K (A. rosae)
Unigene28130	KY421217	1169	71	5e−168	CYP6U (O.brumata)
Unigene32755	KY421218	1752	58	0.0	CYP6 (P. xylostella)
CL2501.Contig1	KY421219	1283	62	2e−158	CYP9A (B. mori)
CL6760.Contig1	KY421220	1430	67	0.0	CYP9E (A. transitella)
Unigene29849	KY421221	1586	53	2e−180	CYP338A (C. medinalis)

Glutathione S-transferases

The analysis of the G. molesta transcriptome data set led to the identification of 28 putative GSTs (Supp Table 3 [online only]). The lengths of the identified GST unigene sequences range from 370 to 5,048 bp encoding putative protein residues from 74 to 702 amino acids. In comparison, the total number of GST unigenes in G. molesta is more than B. mori and less than other three model insect species (Table 2). Phylogenetic analysis revealed the distribution of these GST unigenes in the seven main classes: Delta, Epsilon, Omega, Sigma, Theta, Zeta, and Microsomal GSTs (Fig. 5).

Fig. 5.

Phylogenetic analysis of GSTs in Grapholita molesta (Gm) with other insects. Bm, Bombyx mori; Dm, Drosophila melanogaster; Sl, Spodoptera littoralis; Cq, Culex quinquefasciatus; Px, Papilio xuthus.

All the GST unigenes were mapped to the KEGG pathways (Supp Table 3 [online only]). They fell into the pathways “metabolic pathways (Ko01100)”, “drug metabolism-Cytochrome P450 (Ko00982)” and “flutathione metabolism (Ko00480).” Twenty-four GST unigenes mapped to the KEGG pathways were further annotated using Blast2go and were assigned to five classes involved in xenobiotics metabolism (Table 4). Twelve GST unigenes belonged to Delta and Epsilon classes. These two GST classes have been implicated in resistance to diverse insecticides, such as organophosphates, organochlorines, and pyrethroid insecticides (Ranson et al. 2001, Wang et al. 2008). For instance, in A. aegypi, overexpression of one Delta GST gene and several Epsilon GST genes caused the insect resistance to DDT and pyrethroid (Lumjuan et al. 2005); In B. dorsalis, three Epsilon class GSTs have been shown to contribute to the resistance to malathion. Overexpression of the genes showed increased malathion detoxification capabilities in B. dorsalis and repression of their RNA level led to increased susceptibility to malathion (Lu et al. 2016). L. magratoria LmGSTd1 gene might be involved in the chlorpyrifos resistance, treatment with chlorpyrifos induced the mRNA and protein expression of the LmGSTd1 (Qin et al. 2014).

Table 4.

Putative glutathione S-transferases revealed from Grapholita molesta transcriptome database

Unigene ID	Accession No.	Length (bp)	Identity (%)	E-value	Similar subfamily by BlastX
CL810.Contig1	KY421222	1242	59	2e−96	GST Epsilon (C. medinalis)
Unigene22800	KY421223	843	60	3e−101	GST Epsilon (C. medinalis)
CL5817.Contig1	KY421224	988	61	2e−92	GST Epsilon (C. medinalis)
CL2089.Contig3	KY421225	847	48	3e−62	GST Epsilon (C. medinalis)
Unigene1668	KY421226	487	49	1e−53	GST Epsilon (D. plexippus)
Unigene31881	KY421227	763	55	9e−87	GST Epsilon (D. plexippus)
Unigene18830	KY421228	728	57	3e−88	GST Epsilon (B. mori)
CL3473.Contig2	KY421229	697	59	2e−94	GST Epsilon (B. mori)
CL1195.Contig2	KY421230	1293	83	4e−143	GST (B. mori)
Unigene18898	KY421231	651	72	8e−50	GST Delta (B. mori)
Unigene25519	KY421232	715	80	3e−102	GST Delta (S. litura)
CL330.Contig1	KY421233	934	56	2e−76	GST Delta (C. suppressalis)
Unigene33518	KY421234	2083	65	6e−89	GST Delta (P. xuthus)
Unigene25468	KY421235	847	84	2e−128	GST theta (A. pernyi)
Unigene27154	KY421236	899	71	7e−118	GST theta (A. transitella)
CL3096.Contig1	KY421237	1051	98	4e−153	GST Zeta (C. suppressalis)
Unigene21539	KY421238	720	66	2e−96	GST Zeta (S. litura)
Unigene31172	KY421239	1009	89	4e−158	GST Omega (S. litura)
Unigene3868	KY421240	2923	72	1e−129	GST Omega (B. mori)
Unigene35251	KY421241	958	72	2e−138	GST Omega (B. mori)
CL1012.Contig1	KY421242	958	64	7e−115	GST Omega (D. plexippus)
Unigene17494	KY421243	370	34	4e−12	GST Sigma (D. plexippus)
Unigene19895	KY421244	674	81	3e−124	GST Sigma (C. fumiferana)
CL9209.Contig8	KY421245	1517	65	2e−93	GST Sigma (C. medinalis)

Carboxylesterases

In this study, 44 unigenes were found encoding putative CarEs. The length of the CarE unigenes, ranged from 207 to 2,652 bp, encoding putative proteins with 56 to 697 amino acid residues (Supp Table 4 [online only]). Based on phylogenetic analysis with other known insect CarEs, 24 CarE unigenes identified in G. molesta were assigned to seven classes: seven in class A, three in clade C, nine in clade E, two in clade M, one in each H, I, and N clade. Different from other insect species, no CarE unigenes in G. molesta belonged to B, D, F, G, and J-L clades (Fig. 6). Clades A-C are among the dietary class, which appears to have broad substrate specificities and general dietary or detoxification functions. This class is mainly composed of α-esterases. Some of α-esterases in the dipteran insects are involved in organophosphate resistance (DeSilva et al. 1997, Guerrero 2000). Clades D-H contain secreted and catalytically active esterases, whereas the third class including clades I-N, which contains neuro/developmental esterases (Oakeshott et al. 2010). The study on B. mori have found that there were more CarE genes in the clade A than in other clades where they were also suggested to be involved in coping with xenobiotics in the food (Yu et al 2009). Our sequencing data in G. molesta are consistent with this notion for Clade A. It is possible that the clade A CarEs in G. molesta are the major CarEs involved in the metabolism of foreign chemicals of natural and synthetic origin.

Fig. 6.

Phylogenetic analysis of CarEs from Grapholita molesta (Gm) with other insects. Bm, Bombyx mori; Dp, Danaus plexxippus; Ha, Helicovera armigera; Nv, Nasonia vitripennis; Ms, Manduca sexta; Ar, Athalia rosae; Ag, Anopheles gambiae; Sl, Spodoptera littoralis.

Analysis of KEGG pathways revealed that these CarE unigenes are putatively involved in “metabolic pathways,” “insect hormone biosynthesis,” “drug metabolism-other enzymes” (Supp Table 4 [online only]). The CarE unigenes that potentially involved in metabolism-related pathways were further searched against the NCBI databases by BlastX. Most of these CarE unigenes showed similarity to esterase FE4 and E4 in Myzus persicae (Sulzer) (Table 5). The overproduction of CarE E4 and FE4 was one of the major insecticide resistance mechanisms in M. persicae (Field and Devonshire 1998, Guillemaud et al. 2003, Rivi et al. 2013). Biochemical assays suggested that temephos resistance was associated with elevated CarEs in Aedes albopictus with the significant upregulation of CCEae6a and CCEae3a genes in the temephos resistance (Poupardin et al. 2014, Grigoraki et al. 2015).

Table 5.

Putative CarEs revealed from Grapholita molesta transcriptome database

Unigene ID	Accession No.	Length (bp)	Identity (%)	E- value	Similar subfamily by BlastX
CL4488.Contig1	KY421246	1818	82	0.0	Antennal esterase CXE6 (C. pomonella)
CL976.Contig1	KY421247	2652	66	0.0	Antennal esterase CXE7 (C. pomonella)
Unigene24026	KY421248	622	91	2e−84	CarE CXE10 (C. pomonella)
Unigene24027	KY421249	1238	82	0.0	CarE CXE10 (C. pomonella)
Unigene25181	KY421250	784	77	3e−134	CarE type B (C. pomonella)
Unigene35543	KY421251	1754	80	0.0	CarE CXE8 (C.pomonella)
Unigene34503	KY421252	366	71	4e−51	CarE CXE21, partial (S.littoralis)
Unigene32047	KY421253	2110	66	0.0	Carboxyl/choline CCE025a (H. armiqera)
Unigene35794	KY421254	1199	90	0.0	Esterase B1 (P. xuthus)
CL2290.Contig1	KY421255	2564	55	0.0	Esterase B1 (P. xuthus)
CL888.Contig2	KY421256	1352	51	2e−148	Esterase B1 (P. xuthus)
Unigene34501	KY421257	1219	57	1e−137	Esterase FE4 (P.xuthus)
Unigene33813	KY421258	704	65	1e−108	Esterase E4-like (P. xuthus)
Unigene21050	KY421259	853	42	2e−74	Esterase E4 (P. xylostella)
Unigene1259	KY421260	340	61	3e−141	Esterase E4-like (A. transitella)
CL1669.Contig1	KY421261	2398	72	0.0	Juvenile hormone esterase-like isoform X1 (A. transitella)
Unigene5391	KY421262	414	53	2e−25	Venom CarE-6 (P. xylostella)
Unigene7697	KY421263	928	71	2e−141	Venom CarE-6-like isoform X2 (A. transitella)

Identification of Insecticide Targets

Target-site insensitivity was another insecticide resistance mechanism, which is often caused by point mutations in target protein (Puggioni et al. 2016). Unigenes with high sequence similarity to genes encoding insecticide targets or having functions in insecticide metabolism were also identified in G. molesta (Table 6). These unigenes included two encoding putative AchE, 10 encoding nACHR, one encoding putative ACCase, four encoding putative VGSC, six encoding putative GABA receptor, two encoding putative GluCl, and six encoding putative RyR. These unigenes were further annotated using Blast2go and were assigned to their corresponding subfamily separately (Table 6).

Table 6.

Putatively identified insecticide targets in Grapholita molesta

Unigene ID	Accession No.	Length (bp)	Target type	Insecticide class	Similar subfamily by BlastX
Unigene14021	KY421264	744	Acetylcholines	Organophosphates,	AChE2
Unigene9773	KY421265	612	terase (AchE)	Carbamates	AChE1
CL6055.Contig1	KY421266	1031	nACHR	neonicotinoids	AChR Subunit alpha3
CL6338.Contig2	KY421267	2707			AChR Subunit alpha 7
Unigene13470	KY421268	338			AChR Subunit alpha 9
Unigene18329	KY421269	391			AChR Subunit beta 1
Unigene2511	KY421270	1091			AChR Subunit beta 1
Unigene26844	KY421271	819			AChR Subunit alpha 9
Unigene31962	KY421272	1486			AChR Subunit alpha 9
Unigene32054	KY421273	1186			AChR Subunit alpha 9
Unigene39496	KY421274	254			AChR Subunit alpha 3
Unigene7358	KY421275	366			AChR Subunit alpha 3
CL560.Contig1	KY421276	7277	ACCase	Tetronic & Tetramic acid derivatives	ACCase X4
Unigene38784	KY421277	387	VGSC	Pyrethroids, Pyrethrins	VGSC Subunit alpha
Unigene40058	KY421278	329			VGSC Subunit alpha
Unigene5001	KY421279	412			VGSC Subunit alpha
Unigene9007	KY421280	673			VGSC Subunit alpha
Unigene3694	KY421281	319	GABA receptor	Organochlorines, Phenylprazoles (Fiproles)	GABA Subunit alpha
Unigene38533	KY421282	280			GABA Subunit beta
Unigene39781	KY421283	277			GABA Subunit alpha
Unigene492	KY421284	2266			GABA
Unigene5732	KY421285	230			GABA Subunit beta
Unigene7667	KY421286	267			GABA Subunit beta
CL3467.Contig1	KY421287	1199	Glutamate-	Organochlorines,	GluCl
Unigene23454	KY421288	262	gated chloride channel (GluCl)	Phenylprazoles (Fiproles)	GluCl
Unigene9483	KY421289	398	RyR	Diamides (chlorantraniliprole, cyanthraniliprole, flubendiamide)	RyR 44F
CL1038.Contig1	KY421290	528			RyR 44F
CL2588.Contig1	KY421291	10452			RyR 44F
CL4350.Contig3	KY421292	1196			RyR 44F
Unigene28442	KY421293	1394			RyR 44F
Unigene37417	KY421294	508			RyR 44F

Insecticide resistance has frequently associated with point mutations of select positions at different target receptors. Spinosad resistance in Drosophila and spinosyn resistance in Tuta absoluta (Meyrick) were associated with a G275E mutation in nAChR α6 subunit (Zimmer et al. 2016, Silva et al. 2016). The former was proved by using the gene editing CRISPR/Cas9 system. Imidacloprid resistance in a rice crop pest, N. lugens, was found to be caused by decreased nAChR subunit NIα8 expression levels (Zhang et al. 2016). Similarly in N. lugens, fipronil resistance was affected by two target-site mutations which were synergistic and compensatory in the GABA receptor RDL (Zhang et al. 2016). In A. albopictus, F1534S mutation in VGSG has been reported to be associated with deltamethrin resistance (Xu et al. 2016). Multi-country survey revealed the relationship between the G4946E mutation in RyR) and diamide resistance in P. xylostella (Guo et al. 2014, Steinbach et al. 2015, Troczka et al. 2015). In the present study, we have identified a list of unigenes encoding insecticide targets in G. molesta. How these unigenes function with regards to resistance to insecticides is yet to be investigated.

Summary

To our knowledge, this is the most comprehensive transcriptome study in G. molesta reported so far. It provides valuable resources for genomic and molecular study of G. molesta, an economically important insect pest. The large quantity of SSRs identified in this research will be valuable information for the genetic studies of insecticide related issues and many other aspects in G. molesta. Analysis of the transcriptome data has uncovered the major detoxification related gene families in G. molesta. G. molesta has relatively more P450s in CYP3 and CYP4 clades than other insect species compared. These G. molesta P450s in the clade CYP3 are likely to play an important role in the detoxification of insecticides. In comparison, the total number of GSTs in G. molesta is more than B. mori and less than other three model insect species. These unigenes in the sigma GST subfamily are likely to be involved in insecticide detoxification. There were much more CarEs in the clade A than other clades in G. molesta. It is possible that the clade A CarEs in G. molesta are the major CarEs involved in the metabolism of foreign chemicals of natural and synthetic origin. These major detoxification enzymes, P450s, GSTs, and CarEs, identified in G. molesta potentially have similar functions to other insect species and will be further investigated in the role of insecticide metabolism. Our results may provide a rich transcriptome platform for this fruit-boring insect pest and important information for researchers to understand the molecular mechanism of insecticide resistance in G. molesta.

Click here for additional data file.