Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Transcriptomic and Proteomic Responses of Sweetpotato Whitefly, Bemisia tabaci, to Thiamethoxam

  • Nina Yang,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Wen Xie,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Xin Yang,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Shaoli Wang,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Qingjun Wu,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Rumei Li,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Huipeng Pan,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Baiming Liu,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Xiaobin Shi,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Yong Fang,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Baoyun Xu,

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

  • Xuguo Zhou ,

    xuguozhou@uky.edu (XGZ); zhangyoujun@caas.cn (YJZ)

    Affiliation Department of Entomology, University of Kentucky, Lexington, Kentucky, United States of America

  • Youjun Zhang

    xuguozhou@uky.edu (XGZ); zhangyoujun@caas.cn (YJZ)

    Affiliation Department of Plant Protection, Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences, Beijing, P. R. China

Abstract

Background

The sweetpotato whitefly, Bemisia tabaci (Hemiptera: Aleyrodidae), is one of the most widely distributed agricultural pests. Although it has developed resistance to many registered insecticides including the neonicotinoid insecticide thiamethoxam, the mechanisms that regulate the resistance are poorly understood. To understand the molecular basis of thiamethoxam resistance, “omics” analyses were carried out to examine differences between resistant and susceptible B. tabaci at both transcriptional and translational levels.

Results

A total of 1,338 mRNAs and 52 proteins were differentially expressed between resistant and susceptible B. tabaci. Among them, 11 transcripts had concurrent transcription and translation profiles. KEGG analysis mapped 318 and 35 differentially expressed genes and proteins, respectively, to 160 and 59 pathways (p<0.05). Thiamethoxam treatment activated metabolic pathways (e.g., drug metabolism), in which 118 transcripts were putatively linked to insecticide resistance, including up-regulated glutathione-S-transferase, UDP glucuronosyltransferase, glucosyl/glucuronosyl transferase, and cytochrome P450. Gene Ontology analysis placed these genes and proteins into protein complex, metabolic process, cellular process, signaling, and response to stimulus categories. Quantitative real-time PCR analysis validated “omics” response, and suggested a highly overexpressed P450, CYP6CX1, as a candidate molecular basis for the mechanistic study of thiamethoxam resistance in whiteflies. Finally, enzymatic activity assays showed elevated detoxification activities in the resistant B. tabaci.

Conclusions

This study demonstrates the applicability of high-throughput omics tools for identifying molecular candidates related to thiamethoxam resistance in an agricultural important insect pest. In addition, transcriptomic and proteomic analyses provide a solid foundation for future functional investigations into the complex molecular mechanisms governing the neonicotinoid resistance in whiteflies.

Background

Neonicotinoids, a relatively new class of synthetic insecticides, have been used primarily to control aphids, leafhoppers, whiteflies, and other sap-sucking pests [1], [2]. Imidacloprid and thiamethoxam, two major neonicotinoids on current market, specifically act on the insect's nicotinic acetylcholine receptors (nAChR) within the central nervous system [2], [3]. Extensive and repetitive use of neonicotinoids have led to the development of resistance in the fruit fly Drosophila melanogaster, the green peach aphid Myzus persicae, the house fly Musca domestica L, the brown planthopper Nilaparvata lugens [4], [5], [6], [7], and the sweetpotato whitefly Bemisia tabaci [8], [9].

Bemisia tabaci is an important pest of arable and horticultural crops in the temperate regions of the world. It has a broad host range and can cause tremendous damages directly by feeding and indirectly by transmitting 115 species of begomoviruses [10], [11]. Bemisia tabaci complex contains more than 24 morphologically indistinguishable biotypes, although recent studies have suggested these are genetically distinct cryptic species [12], [13], [14]. During the past two decades, B. tabaci biotype B, one of the most invasive and destructive species/biotypes, has been introduced into at least 54 countries from its origin in the Middle East-Asia Minor region and become a world-wide invasive whitefly species [12], [13]. In China, B. tabaci was first recorded in the late 1940s [15]. However, the crop damages and economic losses caused by this phloem-feeding insect had not been severe until the introduction of B. tabaci B biotype in the mid-1990s [16]. Since then, B biotype has quickly displaced the indigenous B. tabaci populations, rapidly invaded the entire country, and has led to serious yield losses in many crops [17].

The management of B. tabaci has been relied heavily on synthetic insecticides. As a result, pesticide resistance has been developed in B. tabaci in many parts of the world. In Israel and Spain, for example, B. tabaci field populations were found highly resistant to thiamethoxam [8], [18], [19]. In Crete, B. tabaci developed over 1,000-fold resistance to imidacloprid in comparison to its susceptible counterpart in the field [20]. In China, field collected B. tabaci has developed high level of resistance to both imidacloprid and thiamethoxam [21].

Study of insecticide resistance relies heavily on detailed biochemical, genetic, and molecular analyses. In general, the development of insecticide resistance involves one of the following mechanisms: 1) the over-expression of enzymes that break down or bind to (sequester) the pesticide; 2) target-site modifications (mutations) that reduce sensitivity to the insecticide; or 3) reduced penetration of the pesticide through the insect cuticle [7], [22], [23]. For example, a point mutation (Y151S) in two nAChR subunits led to the development of neonicotinoid resistance in N. lugens [7]. Other studies have suggested that the resistance to neonicotinoid pesticides was associated with over-expression of detoxification enzyme cytochrome P450 CYP6G1 in D. melanogaster [4]. Similarly, resistance to neonicotinoid insecticides in the green peach aphid M. persicae has been linked with the over-expression of CYP6CY3 [5].

In B. tabaci, neonicotinoid resistance has yet to be correlated with target-site modification, but has been associated with elevated expression of specific genes. For example, over-expression of a cytochrome P450, CYP6CM1, has been linked to imidacloprid resistance in B. tabaci [24]. In addition, B. tabaci resistance to neonicotinoids, especially to thiamethoxam, has been associated with elevated activities of detoxification enzymes [25], [26]. Most recently, the molecular basis of thiamethoxam resistance in B. tabaci was investigated using the suppression subtractive hybridization (SSH) cDNA library approach [27]. Based on the results of the differential screening, 298 and 209 clones were picked and sequenced, respectively, from the forward and reverse cDNA libraries, representing up- and down-regulated genes between the thiamethoxam-resistant and -susceptible B. tabaci. BLASTX analysis matched 24.5% of these differentially expressed transcripts to genes with known functions. Among them, a NAD-dependent methanol dehydrogenase-like EST from B.tabaci was substantially overexpressed in the resistant whiteflies (∼12-fold). Despite recent progresses, molecular mechanisms underlying B. tabaci resistance to neonicotinoids remain poorly understood.

Functional genomics and proteomics provide an unprecedented opportunity for scientific community to gain a global understanding of the spatial and temporal dynamics of molecular and cellular processes in a living organism and notably facilitate the analysis of genetics and metabolic pathways governing these processes [28], [29]. Also, the application of genomic tools to previously intractable cases of insecticide resistance should greatly expand our understanding of how insecticide resistance evolves and can be avoided or managed [30], [31]. In this study, we compared thiamethoxam resistant and susceptible B. tabaci at transcriptome (RNA-seq) and proteome (iTRAQ) level to enhance our genetic and molecular understanding of the thiamethoxam resistance in an agriculturally important insect pest.

Results

Omics analyses

Thiamethoxam susceptible and resistant B. tabaci transcriptomes were sequenced individually, generating approximately 12 million raw reads for each library (Table 1). After removal of the low quality reads, the total number of clean reads per library ranged from 11–12 million. To reveal the molecular events underlying transcriptomic profiles, sequence reads were mapped to a reference transcriptome containing both B and Q biotypes [32], [33]. Among the two RNA-seq libraries, 77.69% and 79.03% of reads were mapped to a gene in the reference database with a perfect match ratio of 64.55% and 65.45%, respectively (Table 1). The percentage of reads mapped to unique genes were approximately 40%; representing a critical subset of RNA-seq libraries.

thumbnail
Table 1. Summary of RNA-seq metrics from B. tabaci transcriptomes.

https://doi.org/10.1371/journal.pone.0061820.t001

The total number of mass spectrometry detected in B. tabaci proteomes was 39316, representing 5765 peptide spectra and 2226 distinct peptides (Table 2, Dataset S1). Of the 1005 peptides identified, more than 70% (711) were assigned to a putative protein by homology search against non-redundant (NR) database; leaving approximately 30% (297) of the peptides unidentified. Among these annotated proteins, 372 were either hypothetical, putative, or predicted.

thumbnail
Table 2. Summary of iTRAQ metrics from B. tabaci proteomes.

https://doi.org/10.1371/journal.pone.0061820.t002

Differentially expressed genes between susceptible and resistant B. tabaci

Following thiamethoxam treatment, a total of 664 and 674 up- and down-regulated transcripts, respectively, were differentially expressed (FDR≤0.001 and |log2Ratio|≥1) between susceptible (TH-S) and resistant (TH-2000) whiteflies (Figure 1A). Majority of these transcripts (400, ∼66%), however, expressed within 1–2 fold differences (Figure 1B).

thumbnail
Figure 1. Differentially expressed genes between thiamethoxam resistance and susceptible B. tabaci.

(A) All 338 differentially expressed genes between the two B. tabaci strains were selected with a cutoff value of FDR≤0.001 and |log 2 Ratio|≥1. (B) The distribution of differentially expressed genes based on fold of changes.

https://doi.org/10.1371/journal.pone.0061820.g001

Table S1 showed the GO classification of 1338 differentially expressed transcripts between TH-S and TH-2000 whiteflies (2∼fold change, FDR≤0.001). Using Blast2GO, 186 differentially expressed transcripts were able to assign to 37 GO classes (Figure 2A). Majority of these genes were assigned to metabolic process, cellular process, multi-organismal process, response to stimulus, binding, catalytic, organelle, cell, and cell part (Table S1). In the up- regulated group, most of genes were assigned to the cathepsin b-like proteinase, NADH-dehydrogenase, glutathione-S-transferase (GST) genes. To investigate their biological functions, 318 differentially expressed genes were mapped to 160 pathways in the KEGG database. Among them, 40 pathways were substantially enriched (p-value≤0.05), such as “Metabolic pathway” and “Drug metabolism pathway” (Table 3). Specifically, 21 genes encoding enzymes in drug metabolism pathway were highly enriched, including 4 cytochrome P450s and 5 GSTs (Table S2). Interestingly, we also found 53 up-regulated and 65 down-regulated genes in metabolic pathways (Table S3). Up-regulated transcripts included GST, cytochrome P450, glucosyl/glucuronosyl transferases, udp glucuronosyl transferases, nadh dehydrogenase, arginine kinase and cytochrome c oxidase, whereas down-regulated ones were 1-acyl-sn-glycerol-3-phosphate acyltransferase, hemocyanin subunit, NADH- dehydrogenase, aconitase and cytochrome p450. It is worth noting that some of these differentially expressed genes were linked to central nervous system (CNS) diseases pathways, such as Parkinson's disease and Alzheimer's disease. A plausible explanation is that nAChRs, the target site for neonicotinoid insecticides, reside in the CNS.

thumbnail
Figure 2. Gene Ontology classification of differentially expressed genes and proteins between thiamethoxam resistance and susceptible B.

tabaci. The differentially expressed genes or proteins are grouped into three hierarchically-structured GO terms, biological process, cellular component, and molecular function. The y-axis indicates the number of genes or proteins in each GO term. (A) Differentially expressed genes identified by RNA-seq. (B) Differentially expressed proteins identified by iTRAQ.

https://doi.org/10.1371/journal.pone.0061820.g002

thumbnail
Table 3. Significantly enriched KEGG pathways in B. tabaci transcriptome.

https://doi.org/10.1371/journal.pone.0061820.t003

Thiamethoxam-induced differential protein expression between susceptible and resistant B. tabaci

After challenged with thiamethoxam, 52 differentially expressed proteins (p-value≤0.05) were identified between susceptible (TH-S) and resistant (TH-2000) B. tabaci. Among them, 38 proteins were up-regulated (≥1.2∼Fold, p-value≤0.05) and 14 proteins were down-regulated (≤0.8∼Fold, p-value≤0.05) (Table S4). Following in-gel digestion by trypsin, peptides were identified by liquid chromatography-electrospray ionization/multistage mass spectrometry (LC-ESI-MS/MS; Table S4). Glutathione- S-transferase and glucosyl/glucuronosyl transferase, which are involved in xenobiotics detoxification, were up-regulated by 1.96 and 1.56-fold, respectively, in the resistant TH-2000 B. tabaci relative to the susceptible TH-S strain. Other up-regulated peptides in the resistant B. tabaci included UDP glucuronosyl transferase (1.73-fold), implicated in the inactivation and excretion of both endogenous and exogenous compounds [34]; luciferin regenerating enzyme (1.71-fold), playing an important role in the recycling of oxyluciferin into luciferin [35]; eukaryotic initiation factor (1.23-fold), associated with protein translation initiation and elongation processes [36]; glycyl-tRNA synthetase (1.38-fold), involved in RNA modification, RNA transportation, and amino acid-tRNA synthesis [37]; and ADP/ATP translocase proteins (1.49-fold), a group of enzymes catalyzing the exchange of ADP and ATP across the mitochondrial inner membrane [38]. In addition, proteins related to energy regulation, protein transportation and binding were also differentially expressed between the resistant and susceptible B. tabaci (Table S4).

To correlate protein with mRNA expression profiles, accession numbers from the proteomic dataset was extracted and compared to the annotated RNA-seq libraries. Differentially expressed peptides were compared to the nucleotide sequences using BLASTp. Differentially expressed peptides were compared to the nucleotide sequences using BLASTp as follows: (i) E-value is less than 1e-15, (ii) the number of mismatches is no more than 6%, and (iii) the alignment is at least 30 amino acids. Table S5 shows the directional correlation between mRNAs and proteins, and the correlation coefficient between the proteins and gene expression profiles was 0.6643 (Figure 3, Table S5).

thumbnail
Figure 3. Correlation between the differently expressed proteins and genes.

Scatter plots illustrates the distribution of differentially expressed proteins and related genes. The Pearson correlation coefficient between proteins and mRNA expression profiles is shown in the upper left corner of the plot.

https://doi.org/10.1371/journal.pone.0061820.g003

Gene ontology and pathway analysis

Among the 52 differentially expressed proteins, 28 were sub-categorized into 31 hierarchically-structured GO classes including 18 Biological Process, 8 Cellular Component, and 5 Molecular Function (Figure 2B). Specifically, “biological regulation” (9, 32.1%), “cellular process” (17, 60.7%), and “metabolic process” (14, 50%) were highly represented in “Biological Process”, whereas, “cell” (19, 67.9%), “cell part” (19, 67.9%), and “organelle” (12, 43%) were the most common categories in “Cellular Component”. In “Molecular Function”, the top three categories are “catalytic activity” (18, 64.3%), “binding” (10, 35.7%) and “transporter activity” (5, 17.9%) (Figure 2B, Table S6). Upon further investigation by GO enrichment analysis (p-value≤0.05), sodium/potassium-transporting ATPase subunit beta-2, mitochondrial import inner membrane translocase subunit TIM44, thioredoxin, serpin, cathepsin L, ruvB-like protein 1 isoform 1, transport protein sec23, mitochondrial ATP synthase coupling factor 6 precursor, 14-3-3 protein epsilon, peroxiredoxin-like protein, ribosomal protein S15Aa-isoform D, mitochondrial ATP synthase F chain, p270 proteins, GST, and cytochrome P450 CYP3A were significantly impacted by thiamethoxam treatment. To investigate which biological pathways were active when exposed to thiamethoxam, 35 differently expressed proteins were assigned to the reference pathways in KEGG. As a result, 14 pathways were enriched (p-value≤0.05), including “Metabolism of xenobiotics by cytochrome P450”, “Drug metabolism- cytochrome P450”, and “Drug metabolism-other enzymes”. Among them, “Metabolism of xenobiotics by cytochrome P450”, and “Drug metabolism- cytochrome P450” had the lowest p-value (Table 4). KEGG pathway analysis also revealed that the most-enriched peptides, including GSTs, UDP-glucuronosyltransferases, glucosyl/glucuronosyl transferases, and cytochrome P450s, were involved in xenobiotic metabolism (Table S7 and Table S8).

thumbnail
Table 4. Significantly enriched KEGG pathways in B. tabaci proteome.

https://doi.org/10.1371/journal.pone.0061820.t004

qRT-PCR validation study

To validate results from transcriptomic and proteomic analyses, genes encoding differentially expressed proteins between thiamethoxam-resistant and -susceptible B. tabaci were subjected to the qRT-PCR analysis. Among the genes tested (41), 75% (31) were in agreement with the omics results (Table S9). Expression profiles of representative transcripts encoding UDP-glucuronosyltransferase 2B10-like isoform 1 (XP_002704642.1), glucosyl/glucuronosyl transferases (XP_969321.2), GST (ACH90394.1), and cytochrome P450 CYP6CX1 (ACT78507.2), respectively, were shown in Figure 4. All four transcripts showed significantly higher expressions (p-value<0.05) in the resistant TH-2000 strain in comparison to the susceptible TH-S strain, in which CYP6CX1 had the highest differential expression (∼11-fold).

thumbnail
Figure 4. Quantitative real-time PCR analysis.

The relative gene expression of selected target genes was normalized to a BestKeeper composed of endogenous reference genes EF-1a and β-actin. Standard errors were generated from the three biological replicates. Asterisks denote significant gene expression differences between resistant and susceptible B. tabaci, as determined by a paired t-tests (* p<0.05, ** p<0.01).

https://doi.org/10.1371/journal.pone.0061820.g004

Enzymatic activity assay

Resistant B. tabaci strains (TH-R and TH-2000) were derived from a field-collected susceptible strain (TH-S), which had no previous exposure to insecticides. Both of these resistant strains have been exposed to thiamethoxam continuously for over 60 generations. Samples collected from TH-2000 were LC80 survivors, whereas TH-R were LC50 survivors. Theoretically, individuals from TH-2000 have a higher level of tolerance to thiamethoxam treatment, and this difference was reflected in the metabolic enzyme activity assays. For example, GST activity of TH-R strain was slightly higher than that of the susceptible TH-S strain, however, the difference was not substantial (one-way ANOVA, P<0.05; LSD test). In contrast, GST activity was significantly elevated in TH-2000 strain, which exhibited a 1.79-fold higher GST activity than TH-S when using CDNB and GSH as substrates (Table 5). The same trend was observed in P450 monooxygenase activities as well (Table 5). The PNOD activities of the two resistant TH-R and TH-2000 strains are 1.67 and 2.82-fold, respectively, greater than that of the susceptible TH-S strain. Furthermore, both GST and P450 activities between the two resistant strains (TH-R and TH-2000) were also significantly different (one-way ANOVA, P<0.05; LSD test) (Table 5).

thumbnail
Table 5. Metabolic enzyme activity among thiamethoxam resistance and susceptible B. tabaci.

https://doi.org/10.1371/journal.pone.0061820.t005

Discussion

Insects are surrounded with an environment filled with toxic compounds, and constantly expose to xenobiotics through physical contact, ingestion, and inhalation. The repetitive application of synthetic insecticides in an artificially simplified (e.g., monoculture) agricultural ecosystem, insects develop an intricate detoxification system to maximize their survivorship by metabolizing insecticides and plant secondary compounds into less toxic substances and/or increasing their solubility to facilitate their sequestration [39]. Typically, the metabolic detoxification system in insects consists of three major groups of enzymes. The phase I detoxification enzymes, acting on a broad range of substrates directly to reduce their toxicity, are represented by the cytochrome P450 monooxygenases (P450s). The phase II enzymes, including glutathione S-transferases (GSTs), UDP-glucuronosyltransferases(UGTs), and carboxylesterases (COE), facilitate the excretion of hydrophobic toxic compounds by improving their hydrophilicity. The phase III enzymes are mainly composed of ATP binding cassette (ABC) and other transmembrane transporters that actively pump the conjugated xenobiotics out of the cell. A total of 196 genes encoding Phase I and II detoxification enzymes have been identified in D. melanogaster including 89 P450s, 39 GSTs, 35 COEs, and 33UGT [40], [41], [42]. In this study, the differentially expressed genes identified from transcriptome and proteome are predominantly Phase I and II detoxification enzymes, except one overexpressed Phase III enzyme, ABC sub-family G member 4-like isoform 1. Enzymatic activity assay results from this study and a previous report [26] clearly demonstrate the elevated metabolic enzyme responses toward thiamethoxam in B. tabaci.

Phase I enzyme involved in thiamethoxam resistance

Cytochrome P450s (CYP, EC 1.14.14.1) are an extremely important metabolic system involved in catabolism of various classes of insecticides, and the CYP6 family has been implicated in insecticide resistance more than any other P450 families [43], [44]. Increased expression of Cyp6g1 led to resistance to DDT, neonicotinoids, and insect growth regulators in D. melanogaster [4], [45]. Elevated expression of CYP6D1, CYP6A, and CYP6D3 were associated with pyrethroid and neonicotinoids resistance in M. domestica [6], [46]. A brain-specific cytochrome P450, CYP6BQ9, was determined to be responsible for the deltamethrin resistance in Tribolium castaneum [47]. Two duplicated P450 genes, CYP6P4 and CYP6P9, were over-expressed in pyrethroid resistant Anopheles funestus [48]. In Anopheles gambiae, CYP6P3 was significantly over-expressed in field populations and this P450 was linked directly to the permethrin resistance [49].

In B. tabaci, increased detoxification by cytochrome P450s has been associated with neonicotinoid resistance. CYP6cm1, which was over expressed in several laboratory-selected B. tabaci strains, was considered as a legitimate molecular marker for screening imidacloprid resistance in the field [24]. The expression levels of CYP6C and CYP9F, respectively, were elevated 5–7 fold in imidacloprid resistant B. tabaci [50]. After more than 30 generations of selection, thiamethoxam resistance reached 66-fold in resistant B. tabaci [25], [26]. Synergism studies revealed that piperonyl butoxide, a cytochrome P450 inhibitor, had significant synergistic effects toward thiamethoxam, in which resistant ratio dropped to 3-fold in resistant B. tabaci [26]. In this study, combined data from omics analyses and enzymatic activity assays are in accordance with previous reports, strongly suggesting the involvement of P450 monooxygenases in thiamethoxam resistance in B. tabaci. Specifically, CYP6CX1 exhibited significantly higher mRNA expression in a laboratory-selected thiamethoxam resistant strain. In addition, CYP6CX1 has recently been implicated in B. tabaci resistance to fenvalerate, chlopyrifos, and avermectin in the field [51]. To elucidate its role in thiamethoxam resistance, molecular cloning and RNAi-based functional characterization are currently in progress.

Phase II enzyme involved in thiamethoxam resistance

Similar to Phase I P450s, Phase II detoxification enzymes were also elevated in response to insecticide exposures. Elevated production of Phase II enzymes, GSTs (EC 2.5.1.18), has been documented as a mechanism of resistance to organochlorines, organophosphates (OP), and pyrethroids [52]. Over-expression of GSTe2 led to resistance to both DDT and permethrin in Aedes Aegypti [53]. Elevated GST activity was also linked with acaricide resistance in Tetranychus urticae [54]. MdGST-3 was over expressed in OP-resistant M. domestica [55]. In a lambda-cyhalothrin resistant N. lugens, the expression of nlgst1-1 was highly elevated in comparison to a susceptible strain [56]. In the decamethrin resistant Tenebrio molitor, GST provided a passive protection mechanism by binding to this pyrethroid insecticide [57]. In this study, GSTs in a resistance TH-2000 strain was significantly over-expressed in comparison to the susceptible TH-S strain, and it is consistent with a previous study utilizing the same thiamethoxam resistant and susceptible B. tabaci strains [26]. This result is supported by the differentially expressed GSTs identified from both transcriptomic and proteomic analyses, indicating that GSTs play a central role in the detoxification of thiamethoxam in B. tabaci.

UDP-glucuronosyl transferases (EC 2.4.1.17; UGTs) are a major class of phase II drug metabolizing enzymes that play an important role in detoxifying a large number of xenobiotics through conjugation reactions [58]. Specially, the biological process predominantly involved UGTs is glucuronidation, in which a glucuronic acid moiety is added to a suite of endogenous and/or exogenous substrates to facilitate the removal of some of these compounds [59]. UGTs are highly expressed in liver and they are responsible for the metabolism of over 200 drugs in human body [58], [60]. Insect UGTs, on the other hand, are implicated in the detoxifications of insecticides [34]. UGTs activity was highly expressed in the larval and adult stages in D. melanogaster, and correlated with increased detoxification capability [61]. In Bombyx mori, expressed BmUGT1 can conjugate a wide range of substrates including flavonoids, coumarins and other phenolic compounds, suggesting a role of UGTs in detoxification processes [62]. Besides detoxification, UGTs are also involved in olfaction, endobiotic modulation, UV protection, and sequestration [63]. In this study, UGT transcripts and proteins were increased 3.43 and 1.73-fold, respectively, in response to thiamethoxam exposure. Although there is no direct evidence linking UGTs with neonicotinoid resistance, the apparent over-expression of UGT at both transcription and translation levels in the resistant B. tabaci suggests the potential involvement of UGTs in thiamethoxam resistance.

Mitochondria is an important target site for toxic compounds and its involvement in apoptosis controls life and death decisions in the cell [64], [65]. Mitochondria cytochrome c oxidase (COX) (EC 1.9.3.1) is a complex metalloproteinase that provides a critical function in cellular respiration in both prokaryotes and eukaryotes [66]. In mammalian brain, reduced COX activity is an important factor for graded neuronal depolarization (or hyperexcitability) and neuronal death [67]. Elevated expression of COX was involved in the development of resistance to methotrexate in Chinese hamster ovary cells [68]. In Blattella germanica, COX1 was over-expressed in a permethrin resistant strain [69]. Also, increased expression of CO3 was tightly regulated in response to permethrin treatment in A. Aegypti [70]. In Schistosoma mansoni, SCOX1 was substantially increased in a praziquantel resistant strain [71]. In this study, the COX is over-expressed in the TH-2000 strains, suggesting that cytochrome c oxidase may play a role in thiamethoxam resistant in whiteflies.

Enhanced production of carboxylesterases through gene amplification and/or up-regulation has been demonstrated in resistance to organophosphates, carbamates and pyrethroids in M. persicae, Schizaphis graminum, and N. lugens [72], [73], [74]. Although our transcriptomic and proteomic results did not identify differentially expressed carboxylesterases between resistant and susceptible B. tabaci showed specific activity of carboxylesterases was significantly higher in thiamethoxam resistance strain [26]. This is consistent with a previous study focusing on neonicotinoid cross-resistance in a B. tabaci B biotype in Israel [19], and indicates the potential involvement of carboxylesterases in thiamethoxam resistance. In a previous effort, a total of 72 and 52 up- and down-regulated transcripts were identified from the forward and reverse SSH libraries, respectively [27]. Similar to this study, differentially expressed genes between the thiamethoxam-resistant and -susceptible B. tabaci include, but not limit to, cell communication, response to abiotic stimulus, response to stress, lipid particle, nuclear envelope, cell proliferation, and nutrient reservoir activity. qRT-PCR analysis, however, validated 50% of the randomly selected transcripts showed significant differences. With the advent of Omics Era, transcriptomic and proteomic responses of B. tabaci to thiamethoxam were studied here using RNA-seq and iTRAQ, respectively. Ultra high-throughput omics analyses not only provided 10-fold more differentially expressed transcripts (1338 in comparison to 124), but also identified 52 differentially expressed proteins. Besides drastically improved sequencing output, qRT-PCR validation study indicates that coordinated omics approach (75%) is more accurate than SSH method (50%).

Correlations between mRNA and protein expression

In this study, only a modest set of genes and proteins were concordantly differentially expressed (Figure 5). Expression profiles are highly dynamic, and mRNA expression and protein levels do not always correlate [75], [76]. Discrepancies in the directional changes between transcriptome and proteome are likely due to the single sampling time-point, as the temporal changes in transcript versus protein in vivo are poorly understood [77]. In addition, the regulatory mechanism of transcriptome and proteome are complex, and both turnover and mRNA stability are important factors contributing to translation efficiency [76], [77]. Furthermore, without a fully sequenced B. tabaci genome, difference between the expressed transcripts and the translated proteins will most likely be the norm rather than the exception.

thumbnail
Figure 5. Venn diagram summarizing the proportion of proteins and genes significantly expressed between thiamethoxam resistance and susceptible B. tabaci.

Differentially expressed transcripts and peptides are represented in respective circles. The overlapping region denotes specific transcripts with their corresponding peptides.

https://doi.org/10.1371/journal.pone.0061820.g005

Summary

The exposure to xenobiotics typically induces a global transcriptional response that leads to the over expression of the detoxification machinery in insects. Our previous and current efforts using conventional SSH approach, coordinated omics analyses, and metabolic pathway analysis identified a core set of differentially expressed genes and proteins between the thiamethoxam-resistant and -susceptible B. tabaci. Specifically, the expression of a suite of Phase I and Phase II detoxification enzymes, including cytochrome P450s, UDP-glucuronosyltransferases, and glutathione S-transferases, were substantially elevated, suggesting a metabolic detoxification mechanism underlying the thiamethoxam resistance. Taken collectively, our results point the way to multiple avenues of “omics” approaches into the molecular understanding of insecticide resistance in an agriculturally important insect pest.

Materials and Methods

Ethics statement

Bemisia tabaci B biotype strains used in this study were initially collected in the field at Beijing in 2000, and have been maintained in a greenhouse at the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Sciences. The species in the genus of Aleyrodidae are common agricultural pests and are not included in the “List of Protected Animals in China”. No specific permits were required for the described field studies.

Sample preparation

Thiamethoxam susceptible (TH-S) and resistant (TH-R) B. tabaci strains were maintained on cabbage, Brassica oleracea L., cv. Jingfeng 1 [25]. TH-S strain has been cultured without exposure to any insecticides [25]. In contrast, TH-R has been exposed to thiamethoxam for more than 60 generations, and exhibited >70-fold resistance to thiamethoxam in comparison to TH-S strain [26], [27]. In addition, a total of 10,000 TH-R adults were treated with 2,000 mg/L of thiamethoxam to eliminate the heterozygous individuals. After 48 hours, survivors were collected and designed as the TH-2000. This concentration (LC80) was determined based on a dose-response curve generated from an in vivo toxicity bioassay, in which 80% individuals from the resistance TH-R strain were killed by 2,000 mg/L of thiamethoxam insecticide. Approximately 1,500 TH-S and TH-2000 adults were collected, respectively, and snap-frozen immediately in liquid nitrogen and stored at −80°C.

RNA-seq library preparation and sequencing

Total RNA was isolated from TH-S and TH-2000 whitefly samples, respectively, with Trizol (Invitrogen) and according to the manufacturer's protocol. The quantity and quality of RNA were determined with Nanodrop ND-1000. To remove residual DNA contamination, total RNA was treated with RNase-free DNase I (New England BioLabs). mRNA was purified from 6 ug of total RNA from each sample with Dynal oligo (dT) beads (Invitrogen), and then fragmented using RNA fragmentation kit (Ambion). The first cDNA strand was synthesized using random hexamer primers. The double-stranded cDNA fragments were processed by an end repair using T4 DNA polymerase, Klenow DNA polymerase, and T4 polynucleotide kinase (NEB) followed by a single Adenine base addition using Klenow 3′ to 5′ exo-polymerase, and concluded by ligation with Illumina's adaptor. The products were purified with QiaQuick PCR extraction kit (QiaGen) and enriched by PCR amplification. Finally, the library products were then subjected to sequencing analysis on the flow cell via Illumina HiSeq™ 2000.

Annotation and de novo gene expression

Raw reads were transformed into clean reads by removing adaptor sequences, empty sequences (sequences with only adaptor but no reads), low-quality sequences (reads with unknown sequences [“N”]), and reads with a copy number of 1 (likely a sequencing error). All clean reads were then mapped to references sequences [32], [33] and only allowed no more than 1 nucleotide mismatch. Clean reads that map to multiple genes were filtered, and the remaining clean reads were designated as unambiguous reads. We, therefore, quantified transcript levels in reads per kilobase per million mapped reads (RPKM) [78]. For functional annotation, distinct sequences were BLAST against the NCBI NR database with a cut-off E-value of 10−5. In addition, Blast2GO (http://www.blast2go.org) was used to assign Gene Ontology terms (http://www.geneontology.org), while Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/ or http://www.kegg.jp/), a database resource that integrates genomic, chemical, and systemic functional information, was adopted to annotate molecular networks (pathways).

Screening for differentially expressed genes

The mapped reads were assembled using a software package SOAPdenovo (short oligonucleotide alignment program) [79]. Then, the RPKM value for each transcript was measured in reads per kilobase of transcript sequence per million mapped reads. The significance of digital gene expression profiles were analyzed as described previously [80] with minor modifications. This algorithm, similar to GFOLD [81], denotes the number of unambiguous clean reads from gene A as x, as every gene's expression occupies only a small part of the library. P(x) yields to the Poisson distribution.The fold change of each transcript was then calculated by the formula of log2 (TH-2000_RPKM/TH-S_RPKM). Where N1 and N2 indicate the total number of clean reads in TH-S and TH-2000, respectively, and x and y denote the mapped clean reads in each sample. FDR (False Discovery Rate) method was used to determine the threshold of p-value in differential gene expression tests. The probability that gene A expresses equally between two samples can be calculated with the following formula:In this study, “FDR≤0.001 and the absolute value of log2Ratio≥1” was the threshold to evaluate the significance level of differentiated gene expression.

Protein quantification and database search using iTRAQ labeling

1 mg of whiteflies (approximately 100 individuals) were dissected in lysis buffer (7M urea, 2M thiourea, 4% CHAPS, 40 mM Tris-HCl, pH 8.5), and 1 mM PMSF (phenylmethanesulfonyl fluoride) and 2 mM EDTA (ethylene diamine tetraacetic acid) were added later, respectively. After 5 min, 10 mM DTT was added to the lysis solution, and then it was centrifuged at 4°C, 30,000×g for 15 min. The supernatant was collected, and the concentration of total proteins was determined using a 2D Quantification Kit (GE Healthcare). For quality check, 30 µg of total protein from each sample was subjected to the SDS-PAGE analysis. After that, 100 µl protein from each sample was digested with trypsin gold (Promega) (protein∶ trypsin = 30∶1) at 37°C for 16 h, and resultant peptides were dried by vacuum centrifugation. The peptides were reconstituted in 0.5 M TEAB and processed according to the manufacturer's protocol for 8-plex iTRAQ (Applied Biosystems, Inc). Samples (100 mg total protein/sample) from TH-S and TH-2000 strains were labeled with 113 and 114 iTRAQ tags, respectively. Then pooled mixtures of iTRAQ-labeled peptides were fractionated by SCX chromatography (Phenomenex, Inc, USA) using a Shimadzu LC-20AB HPLC Pump system. Collected fractions were pooled into 10 final fractions and analyzed by nano LC-MS/MS analysis after desalting by Strata XC18 column (Phenomenex) and vacuum-dried. Nano LC-MS/MS analysis on each of these fractions was performed using a LTQ-Orbitrap Velos mass spectrometer (Thermo Fisher Scientific Inc. Rockford, IL., USA) equipped with nanoelectrospray ionization.

Peptides were identified by searching against a combined database containing B. tabaci B and Q biotypes (reference transcriptomes included both 454 and illumina sequencing data) [32], [33] using a MS/MS data interpretation algorithm within Mascot. Searches were restricted to whitefly database with carbamidomethyl cysteine as a fixed modification, and oxidized methionine as a variable modification. A peptide mass tolerance of 3 ppm and fragment mass tolerance of ±0.05 Da with one missed cleavage site of trypsin were allowed. When the Mascot software was used to search the database, 1005 proteins were identified with a false discovery rate (FDR) of less than 1%. Differential expression ratios for proteins were obtained from Mascot software (http://www.matrixscience.com), which calculates protein ratios using only ratios from the spectra that are distinct for each protein and excluding the shared peptides of protein isoforms. To calculate differential expression ratios, all identified spectra from a protein were used to obtain an average protein ratio relative to the control label (i.e. fold change). Student t-test was used to analyze the differential expression of proteins between resistant and susceptible B. tabaci. The p-value was calculated using the confidence intervals from the error factor generated in Mascot aswhere, N is the number of peptide ratios, s is the standard deviation, and x represents the mean of peptide ratios. In this study, we used P≤0.05 and the fold change >1.2-fold or <0.8-fold as the thresholds to judge the significance level of differentiated protein expression.

Functional and pathway enrichment analyses

Functional annotation of transcripts and proteins identified in B. tabaci was carried out using Blast2GO, an integrated GO annotation and data mining tool that assigns gene ontology through BLAST searches against nucleotide and/or protein databases [82]. GO enrichment analysis provides all GO terms that significantly enriched in differentially expressed genes/proteins in comparison to susceptible whiteflies. This method firstly maps all differentially expressed genes/proteins to GO terms in the database (http://www.geneontology.org/), calculating gene numbers for every term, then using hypergeometric test to find significantly enriched GO terms in differentially expressed genes/proteins comparing to B. tabaci transcriptome/proteome background. The calculating formula isWhere N is the number of all genes/or proteins with GO annotation; n is the number of differentially expressed genes/or proteins in N; M is the number of all genes/or proteins that are annotated to the certain GO terms; m is the number of differentially expressed genes/or proteins in M. The calculated p value was subjected to a Bonferroni Correction, taking a corrected p-value of 0.05 as a threshold. GO terms fulfilling this condition were defined as significantly enriched GO terms in differentially expressed genes/or proteins.

All identified transcripts and proteins were mapped to pathway in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Significantly enriched metabolic pathways or signal transduction pathways in differentially expression genes and proteins were identified using the same calculating formula as in GO analysis. Here N is the number of all genes/or proteins that with KEGG annotation, n is the number of differentially expressed genes/or proteins in N, M is the number of all genes/or proteins annotated to specific pathways, and m is number of differentially expressed genes/or proteins in M.

Correlation between protein and mRNA expression

To assess the correlation between transcriptomic and proteomic platforms, we first set cutoff values to select subsets of genes and proteins with distinctive expression signals. All the protein sequences identified by iTRAQ were analyzed and converted into a searchable database. For each protein, we queried the RNA-seq data from the same expression pattern of a matching transcript (P-value<0.05). The significance level of the overlap between detected proteins and transcripts was determined using Pearson's Chi-squared test with Yates' continuity correction.

Quantitative real time PCR (qRT-PCR) analysis

Total RNA was extracted from TH-2000 and TH-S adults, respectively, using Trizol (Invitrogen) following the manufacturer's protocols. The total RNA obtained was resuspended in nuclease-free water and the concentration was measured using Nanodrop (Thermo Scientific Nanodrop 2000). cDNA was synthesized using the SYBR PrimeScript reverse transcription-PCR (RT-PCR) kit (Takara). qRT-PCRs were carried out on an ABI 7500 real-time PCR system (Applied Biosystems) with SYBR Green Real-time PCR Master Mix (TaKaRa) following a cycling regime of 95°C for 3 min, 40 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 35 s. The primers used for real-time PCR are listed in Table S10. Three biological replicates for each sample were used for qRT-PCR analysis. The relative gene expression of selected target genes was normalized to the BestKeeper [83], a composite internal reference gene composed of endogenous EF-1a and β-actin to eliminate sample-to-sample variations. The relative expression levels were calculated using the 2−ΔΔCt method [84].

Enzymatic activity assays

GST activity was measured using the 1-chloro-2,4-dinitrobenzene (CDNB) and reduced GSH as substrates according to Feng [29]. The total reaction volume per well of a 96-well microtiter plate total volume was 300 µL, consisting of 100 µL crude enzyme,100 µL 1-chloro-2,4-dinitrobenzene (CDNB) (1.2 mM), and 100 µL glutathione reduced (GSH) (12 mM). Change in absorbance was monitored continuously at 340 nm for 10 min at 25°C using a thermomax kinetic microplate reader (Molecular Devices). Cytochrome-P450-dependent monooxygenase activity was determined using p-nitroanisole (PNA) as substrates with minor modifications [29]. Incubation mixture contained 375 µL PNA (2 µM), 37.5 µL NADPH (Nicotinamide adenine dinucleotide phosphate) (9.6 mM), and 337.5 µL crude enzyme. The mixture was incubated for 30 min at 34°C in one air atmosphere. The reaction volume per well of a 96-well microtiter plate was 200 µL of incubation mixture. The P450 monooxygenase activity is determined spectrophotometrically by monitoring the accumulation of metabolic products according to a p-nitrophenol standard curve at OD405.

Supporting Information

Dataset S1.

Peptide spectra in B. tabaci proteome.

https://doi.org/10.1371/journal.pone.0061820.s001

(XLSX)

Table S1.

Gene ontology analysis of B. tabaci transcriptome.

https://doi.org/10.1371/journal.pone.0061820.s002

(XLSX)

Table S2.

Differentially expressed genes within drug metabolism pathway in response to thiamethoxam exposure.

https://doi.org/10.1371/journal.pone.0061820.s003

(XLSX)

Table S3.

Differentially expressed genes within metabolism pathways in response to thiamethoxam exposure.

https://doi.org/10.1371/journal.pone.0061820.s004

(XLSX)

Table S4.

Differentially expressed peptides identified by LC-ESI-MS/MS.

https://doi.org/10.1371/journal.pone.0061820.s005

(XLSX)

Table S5.

Differentially expressed peptides and their corresponding mRNAs in response to thiamethoxam exposure.

https://doi.org/10.1371/journal.pone.0061820.s006

(DOCX)

Table S6.

Gene ontology analysis of B. tabaci proteome.

https://doi.org/10.1371/journal.pone.0061820.s007

(XLSX)

Table S7.

Differentially expressed peptides within drug metabolism pathway in response to thiamethoxam exposure.

https://doi.org/10.1371/journal.pone.0061820.s008

(XLSX)

Table S8.

Differentially expressed peptides within metabolism pathways in response to thiamethoxam exposure.

https://doi.org/10.1371/journal.pone.0061820.s009

(XLSX)

Table S9.

qRT-PCR verified differentially expressed genes.

https://doi.org/10.1371/journal.pone.0061820.s010

(XLSX)

Table S10.

Primers used for qRT-PCR analyses.

https://doi.org/10.1371/journal.pone.0061820.s011

(XLSX)

Acknowledgments

Authors are grateful to editor and anonymous reviewers for their critical comments and suggestions. Special thanks go to Dr. John Obrycki (Department of Entomology, University of Kentucky) for his comments on an earlier draft.

Data Deposition

Bemisia tabaci transcriptome datasets are available at NCBI GEO with an accession number GSE42102.

Bemisia tabaci proteome datasets are available at PeptideAtlas under a submission number PASS00135.

Author Contributions

Conceived and designed the experiments: NNY WX XGZ YJZ. Performed the experiments: NNY XBS YF XY. Analyzed the data: NNY. Contributed reagents/materials/analysis tools: NNY RML HPP BML SLW QJW BYX XGZ YJZ. Wrote the paper: NNY XGZ.

References

  1. 1. James DG, Price TS (2002) Fecundity in two-spottd spider mite (Acari:Tetranychidae) is increased by direct and systemic exposure to imidacloprid. J Econ Entomol 95: 729–732.
  2. 2. Tomizawa M, Casida JE (2005) Neonicotinoid insecticide toxicology: Mechanisms of Selective Action. Annu Rev Pharmacol Toxicol 45: 247–268.
  3. 3. Tomizawa M, Casida JE (2003) Selective toxicity of neonicotinoids attributable to specificity of insect and mammalian nicotinic receptors. Annu Rev Entomol 48: 339–364.
  4. 4. Daborn PJ, Yen JL, Bogwitz MR, Le Goff G, Feil E, et al. (2002) A single P450 allele associated with insecticide resistance in Drosophila. Science 297: 2253–2256.
  5. 5. Puinean AM, Foster SP, Oliphant, Denholm I, Field LM, et al. (2010) Amplification of a Cytochrome P450 gene is associated with resistance to neonicotinoid insecticides in the aphid Myzus persicae. PLoS Genet 6: e1000999.
  6. 6. Markussen MDK, Kristensen M (2010) Cytochrome P450 monooxygenase- mediated neonicotinoid resistance in the house fly Musca domestica L. Pestici Biochem Physiol 98: 50–58.
  7. 7. Liu ZW, Williamson MS, Lansdell SJ, Denholm I, Han ZJ, et al. (2005) A nicotinic acetylcholine receptor mutation conferring target-site resistance to imidacloprid in Nilaparvata lugens (brown planthopper). Proc Natl Acad Sci 102: 8420–8425.
  8. 8. Elbert A, Nauen R (2000) Resistance of Bemisia tabaci (Homoptera: Aleyrodidae) to insecticides in southern Spain with special reference to neonicotinoids. Pest Manag Sci 56: 60–64.
  9. 9. Wang ZY, Yao MD, Wu YD (2009) Cross-resistance, inheritance and biochemical mechanisms of imidacloprid resistance in B-biotype Bemisia tabaci. Pest Manag Sci 65: 1189–1194.
  10. 10. Stansly PA, Naranjo SE (2010) Bemisia: Bionomics and management of a global pest. Springer Science+Business Media B V 545.
  11. 11. Hogenhout SA, Ammar El-D, Whitfield AE, Redinbaugh MG (2008) Insect vector interactions with persistently transmitted viruses. Annu Rev Phytopathol 46: 327–359.
  12. 12. Dinsdale A, Cook L, Riginos C, Buckley YM, De Barro P (2010) Refined global analysis of Bemisia tabaci (Hemiptera: Sternorrhyncha: Aleyrodoidea: Aleyrodidae) mitochondrial cytochrome oxidase 1 to identify species level genetic boundaries. Anna Entomol Soci Amer 103: 196–208.
  13. 13. De Barro PJ, Liu SS, Boykin LM, Dinsdale AB (2011a) Bemisia tabaci: A statement of species status. Annu Rev Entomol 56: 1–19.
  14. 14. De Barro PJ, Ahmed MZ (2011b) Genetic networking of the Bemisia tabaci cryptic species complex reveals pattern of biological invasions. PLoS ONE 6: e25579.
  15. 15. Zhou Y (1949) A list of Aleyrodidae from China. Entomol Sinica 3: 1–18.
  16. 16. Luo C, Yao Y, Wang RJ, Yan FM, Hu DX, et al. (2002) The use of mitochondrial cytochrome oxidase mtCOI gene sequences for the identification of biotypes of Bemisia tabaci (Gennadius) in China. Acta Ent omologica Sinica 45: 759–763.
  17. 17. Liu SS, De Barro PJ, Xu J, Luan JB, Zang LS, et al. (2007) Asymmetric mating interactions drive widespread invasion and displacement in a whitefly. Science 318: 1769–1772.
  18. 18. Horowitz AR, Kontsedalov S, Ishaaya I (2004) Dynamics of resistance to the neonicotinoids acetamiprid and thiamethoxam in Bemisia tabaci (Homoptera: Aleyrodidae). J Econ Entomol 97: 2051–2056.
  19. 19. Rauch N, Nauen R (2003) Biochemical markers linked to neonicotinoid cross-resistance in Bemisia tabaci (Hemiptera: Aleyrodidae). Arch Ins Biochem Physiol 54: 165–176.
  20. 20. Roditakis E, Grispou M, Morou E, Kristoffersen J B, Roditakis N, et al. (2008) Current status of insecticide resistance in Q biotype Bemisia tabaci populations from Crete. Pest Manag Sci 65: 313–322.
  21. 21. Wang ZY, Yan HF, Yang YH, Wu YD (2010) Biotype and insecticide resistance status of the whitefly Bemisia tabaci from China. Pest Manag Sci 66: 1360–1366.
  22. 22. Feyereisen R (1995) Molecular biology of insecticide resistance. Toxicol Lett 82/83: 83–90.
  23. 23. Perry T, Batterham P, Daborn PJ (2011) The biology of insecticidal activity and resistance. Insect Biochem Mol Biol 41: 411–422.
  24. 24. Karunker I, Benting J, Lueke B, Ponge T, Morin S, et al. (2008) Over-expression of cytochrome P450 CYP6CM1 is associated with high resistance to imidacloprid in the B and Q biotypes of Bemisia tabaci (Hemiptera: Aleyrodidae). Insect Biochem Mol Biol 38: 634–644.
  25. 25. Feng YT, Wu QJ, Xu BY, Wang SL, Zhang YJ, et al. (2008) Fitness costs and morphological change of laboratory-selected thiamethoxam resistance in the B-type Bemisia tabaci (Hemiptera: Aleyrodidae). J Appl Entomol 133: 466–472.
  26. 26. Feng YT, Wu QJ, Wang SL, Chang XL, Zhang YJ, et al. (2009) Cross-resistance study and biochemical mechanisms of thiamethoxam resistance in B-biotype Bemisia tabaci (Hemiptera: Aleyrodidae). Pest Manag Sci 66: 313–318.
  27. 27. Xie W, Yang X, Wang SL, Yang NN, Zhang YJ, et al. (2012a) Gene expression profiling in the thiamethoxam resistant and susceptible B-biotype sweetpotato whitefly, Bemisia tabaci. J Insect Sci 12: 46–60.
  28. 28. Kislinger T, Cox B, Kannan A, Chung C, Emili A, et al. (2006) Global survey of organ and organelle protein expression in mouse: combined proteomic and transcriptomic profiling. Cell 125: 173–186.
  29. 29. Carolanv JC, Caragea D, Reardon KT, Mutti NS, Edwards OR, et al. (2011) Predicted effector molecules in the salivary secretome of the pea aphid (Acyrthosiphon pisum): a dual transcriptomic/proteomic approach. J Proteome Res 4: 1505–1518.
  30. 30. Oakeshott JG, Horne I, Sutherland TD, Russell RJ (2003) The genomics of insecticide resistance. Genome Biol 4: 202–205.
  31. 31. Liu B, Jiang GF, Zhang YF, Li JL, Ran C, et al. (2011) Analysis of transcriptome differences between resistant and susceptible strains of the citrus red mite panonychus citri (Acari: Tetranychidae). PLoS ONE 6: 1–9.
  32. 32. Wang XW, Luan JB, Li JM, Bao YY, Zhang CX, et al. (2010) De novo characterization of a whitefly transcriptome and analysis of its gene expression during development. BMC Genomics 11: 400–410.
  33. 33. Xie W, Meng QS, Wu QJ, Wang SL, Zhang YJ, et al. (2012b) Pyrosequencing the Bemisia tabaci transcriptome reveals a highly diverse bacterial community and a robust system for insecticide resistance. PLoS ONE DOI:https://doi.org/10.1371/journal.pone.0035181.
  34. 34. Huang FF, Chai CL, Zhang Z, Xiang ZH (2008) The UDP-glucosyltransferase multigene family in Bombyx mori. BMC Genomics 9: 563–577.
  35. 35. Emamzadeh R, Emamzadeh S, Hemati R, Sadeghizadeh M (2010) RACE-based amplification of cDNA and expression of a luciferin-regenerating enzyme (LRE): An attempt towards persistent bioluminescent signal. Enzyme Microb Tech 47: 159–165.
  36. 36. Chaudhuri J, Si K, Maitra U (2012) Function of eukaryotic translation initiation factor 1A (eIF1A) (formerly called eIF-4C) in initiation of protein synthesis. J Biol Chem 272: 7883–7891.
  37. 37. Sang LJ, Gyu PS, Park H, Seol W, Kim S, et al. (2002) Interaction network of human aminoacyl-tRNA synthetases and subunits of elongation factor 1 complex. Biochem Biophys Res Commun 291: 158–64.
  38. 38. Kokoszka JE, Waymire KG, Levy SE, Sligh JE, Wallace DC, et al. (2004) The ADP/ATP translocator is not essential for the mitochondrial permeability transition pore. Nature 427: 461–465.
  39. 39. Misra JR, Horner MA, Lam G, Thummel CS (2011) Transcriptional regulation of xenobiotic detoxification in Drosophila. Genes Dev 25: 1796–1806.
  40. 40. Low WY, Ng HL, Morton CJ, Parker MW, Robin C, et al. (2007) Molecular evolution of glutathione s-transferases in the genus Drosophila. Genetics 177: 1363–1375.
  41. 41. Ranson H, Claudianos C, Ortelli F, Abgrall C, Feyereisen R, et al. (2002) Evolution of supergene families associated with insecticide resistance. Science 298: 179–181.
  42. 42. Tijet N, Helvig C, Feyereisen R (2001) The cytochrome P450 gene superfamily in Drosophila melanogaster: annotation, intron-exon organization and phylogeny. Gene 262: 189–198.
  43. 43. Scott JG (1999) Cytochromes P450 and insecticide resistance. Insect Biochem Mol Biol 29: 757–777.
  44. 44. Puinean AM, Foster SP, Oliphant L, Field LM, Bass C, et al. (2011) Amplification of a cytochrome p450 gene is associated with resistance to neonicotinoid insecticides in the aphid myzus persicae. PLoS Genet 6 (6) e1000999.
  45. 45. Joußen N, Schuphan I, Schmidt B (2010) Metabolism of methoxychlor by the P450-Monooxygenase CYP6G1 involved in insecticide resistance of drosophila melanogaster after expression in cell cultures of nicotiana tabacum. Chem Biodivers 7: 722–735.
  46. 46. Liu NN, Scott JG (1998) Increased transcription of CYP6D1 causes cytochrome P450-mediated insecticide resistance in house fly. Insect Biochem Mol Biol 28: 531–535.
  47. 47. Zhu F, Parthasarathy R, Bai H, Woithe K, Palli SR, et al. (2010) A brain-specific cytochrome P450 responsible for the majority of deltamethrin resistance in the QTC279 strain of Tribolium castaneum. Proc Natl Acad Sci 107: 8557–8562.
  48. 48. Wondji CS, Irving H, Morgan J, Hunt RH, Ranson H, et al. (2009) Two duplicated P450 genes are associated with pyrethroid resistance in Anopheles funestus, a major malaria vector. Genome Res 19: 452–459.
  49. 49. Müller P, Warr E, Stevenson BJ, Pignatelli PM, Donnelly MJ, et al. (2008) Field-caught permethrin-resistant Anopheles gambiae overexpress CYP6P3, a P450 that metabolises pyrethroids. PLoS Genet 4: e1000286.
  50. 50. Qiu BL, Liu L, Li XX, Mathur V, Ren SX (2009) Genetic mutations associated with chemical resistance in the cytochrome P450 genes of invasive and native Bemisia tabaci (Hemiptera: Aleyrodidae) populations in China. Insect Sci 16: 237–245.
  51. 51. Zhuang HM, Wang KF, Zheng L, Wu ZJ, Wu G, et al. (2011) Identification and characterization of a cytochrome P450 CYP6CX1 putatively associated with insecticide resistance in Bemisia tabaci. Insect Sci 18: 484–494.
  52. 52. Li XC, Schuler MA, Berenbaum MR (2007) Molecular mechanisms of metabolic resistance to synthetic and natural xenobiotics. Annu Rev Entomol 52: 231–253.
  53. 53. Lumjuan N, McCarroll L, La-aied Prapanthadara, Hemingway J, Ranson H (2005) Elevated activity of an Epsilon class glutathione transferase confers DDT resistance in the dengue vector, Aedes aegypti. Insect Biochem Mol Biol 35: 861–871.
  54. 54. Stumpf N, Nauen R (2002) Biochemical markers linked to abamectin resistance in Tetranychus urticae (Acari: Tetranychidae). Pest Biochem Physiol 72: 111–121.
  55. 55. Syvanen M, Zhou Z, Wharton J, Goldsbury C, Clark A (1996) Heterogeneity of the glutathione transferase genes encoding enzymes responsible for insecticide degradation in the housefly. J Mol Evol 43: 236–240.
  56. 56. Vontas JG, Small GJ, Dimitra DC, Ranson H, Hemingway J (2002) Purification, molecular cloning and heterologous expression of a glutathione S-transferase involved in insecticide resistance from the rice brown planthopper, Nilaparvatalugens. J Biochem 362: 329–337.
  57. 57. Kostaropoulos I, Papadopoulos AI, Metaxakis A, Boukouvala E, Papadopoulou-Mourkidou E (2001) Glutathione S-transferase in the defence against pyrethroids in insects. Insect Biochem Mol Biol 31: 313–319.
  58. 58. Miley MJ, Zielinska AK, Keenan JE, Bratton SM, Redinbo MR (2007) Crystal structure of the cofactor-binding domain of the human phase II drug-metabolism enzyme UDP -glucuronosyltransferase 2B7. J Mol Biol 369: 498–511.
  59. 59. Hundle BS, O'Baien DA, Alberti M, Beyer P, Hearst P (1992) Functional expression of zeaxanthin glucosyltransferase from Erwinia herbicola and a proposed uridine diphosphate binding site. Proc NatI Acad Sci 89: 9321–9325.
  60. 60. Samokyszyn VM, Gall WE, Zawada G, Freyaldenhoven A, Radominska-Pandya A (2000) 4-Hydroxyretinoic acid, a novel substrate for human liver microsomal UDP-glucuronosyltransferase(s) and recombinant UGT2B7. J Biol Chem 275: 6908–6914.
  61. 61. Rausell C, Llorca J, Real MD (1997) Separation by FPLC chromatofocusing of UDP-glucosyltransferases from three developmental stages of drosophila melanogaster. Arch Insect Biochem 34: 347–358.
  62. 62. Luque T, O'Reilly DR (2002) Functional and phylogenetic analyses of a putative Drosophila melanogaster UDP-glycosyltransferase gene. Insect Biochem Mol Biol 32: 1597–1604.
  63. 63. Ahn SJ, Vogel H, Heckel DG (2012) Comparative analysis of the UDP-glycosyltransferase multigene family in insects. Insect Biochem Mol Biol 42: 133–147.
  64. 64. Hand SC, Menze MA (2008) Commentary mitochondria in energy-limited states: mechanisms that blunt the signaling of cell death. J Exp Biol 211: 1829–1840.
  65. 65. Tournier C, Hess P, Yang DD, Xu J, Davis RJ (2000) Requirement of JNK for stress induced activation of the cytochrome c–mediated death pathway. Science 288: 870–874.
  66. 66. Capaldi RA (1990) Structure and function of cytochrome c oxidase. Annu Rev Biochem 59: 569–96.
  67. 67. Greenamyre JT, MacKenzie G, Peng TI, Stephans SE (1999) Mitochondrial dysfunction in Parkinson's disease. Biochem Soc Symp 66: 85–97.
  68. 68. Alemany C, Noe V, Ciudad CJ (2000) Identification by RNA-based arbitrarily primed PCR of the involvement of cytochrome c oxidase in the development of resistance to methotrexate. Biochem Biophys Acta 1495: 319–326.
  69. 69. Pridgeon JW, Liu NN (2003) Overexpression of the cytochrome c oxidase subunit I gene associated with a pyrethroid resistant strain of German cockroaches, Blattella germanica (L.). Insect Biochem Mol Biol 33: 1043–1048.
  70. 70. Pridgeon JW, Becnel JJ, Clark GG, Linthicum KJ (2009) Permethrin induces overexpression of cytochrome c oxidase subunit 3 in Aedes aegypti. J Med Entomol 46: 810–819.
  71. 71. Pereira C, Fallon PG, Cornette J, Capron A, Pierce RJ, et al. (1998) Alterations in cytochrome-c oxidase expression between praziquantel-resistant and susceptible strains of Schistosoma mansoni. Parasitology 117: 63–73.
  72. 72. Hemingway J (2000) The molecular basis of two contrasting metabolic mechanisms of insecticide resistance. Insect Biochem Mol Biol 30: 1009–1015.
  73. 73. Ono M, Swanson JJ, Field LM, Devonshire AL, Siegfried BD (1999) Amplification and methylation of an esterase gene associated with insecticide-resistance in greenbugs, Schizaphis graminum (Rondani) (Homoptera: Aphididae). Insect Biochem Mol Biol 29: 1065–1073.
  74. 74. Small GJ, Hemingway J (2000) Molecular characterization of the amplified carboxylesterase gene associated with organophosphorus insecticide resistance in the brown planthopper, Nilaparvata lugens. Insect Mol Biol 9: 647–653.
  75. 75. Nie L, Wu G, Zhang W (2006) Correlation between mRNA and protein abundance in desulfovibrio vulgaris: A multiple regression to identify sources of variations. Biochem Biophys Res Commun 339: 603–610.
  76. 76. Bolognani F, Perrone-Bizzozero NI (2008) RNA-protein interactions and control of mRNA stability in neurons. J Neurosci Res 86: 481–489.
  77. 77. Popesku JT, Martyniuk CJ, Denslow ND, Trudeau VL (2010) Rapid dopaminergic modulation of the fish hypothalamic transcriptome and proteome. PLoS ONE 5 (8) e12338 .
  78. 78. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B (2008) Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods 5: 621–644.
  79. 79. Li R, Yu C, Li Y, Lam TW, Yiu SM, et al. (2009) SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25: 1966–1967.
  80. 80. Audic S, Claverie JM (1997) The significance of digital gene expression profiles. Genome Res 7: 986–995.
  81. 81. Feng J, Meyer CA, Wang Q, Liu JS, Zhang Y, et al. (2012) GFOLD: a generalized fold change for ranking differentially expressed genes from RNA-seq data. Bioinformatics .
  82. 82. Conesa A, Götz S, García-Gómez JM, Terol J, Robles M, et al. (2005) Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21: 3674–3676.
  83. 83. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP (2004) Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper – Excel-based tool using pair-wise correlations. Biotechnol Lett 26: 509–515.
  84. 84. Pfaffl MW (2001) A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res 29: e45.