Skip to main content

Transcriptional analysis of Bemisia tabaci MEAM1 cryptic species under the selection pressure of neonicotinoids imidacloprid, acetamiprid and thiamethoxam

Abstract

Background

Neonicotinoids are widely applied in the control of the destructive agricultural pest Bemisia tabaci, and resistance against these chemicals has become a common, severe problem in the control of whiteflies. To investigate the molecular mechanism underlying resistance against nenonicotinoids in whiteflies, RNA-seq technology was applied, and the variation in the transcriptomic profiles of susceptible whiteflies and whiteflies selected by imidacloprid, acetamiprid and thiamethoxam treatment was characterized.

Results

A total of 90.86 GB of clean sequence data were obtained from the 4 transcriptomes. Among the 16,069 assembled genes, 584, 110 and 147 genes were upregulated in the imidacloprid-selected strain (IMI), acetamiprid-selected strain (ACE), and thiamethoxam (THI)-selected strain, respectively, relative to the susceptible strain. Detoxification-related genes including P450s, cuticle protein genes, GSTs, UGTs and molecular chaperone HSP70s were overexpressed in the selected resistant strains, especially in the IMI strain. Five genes were downregulated in all three selected resistant strains, including 2 UDP-glucuronosyltransferase 2B18-like genes (LOC 109030370 and LOC 109032577).

Conclusions

Ten generations of selection with the three neonicotinoids induced different resistance levels and gene expression profiles, mainly involving cuticle protein and P450 genes, in the three selected resistant whitefly strains. The results provide a reference for research on resistance and cross-resistance against neonicotinoids in B. tabaci.

Peer Review reports

Background

The cotton whitefly Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae) is one of the most invasive and destructive pests worldwide. It has a wide host range including more than 600 species of agriculture and horticulture plants. B. tabaci is a complex of at least 41 discrete species that belong to 11 major genetic groups [1,2,3,4,5,6,7,8,9]. The B. tabaci genome is highly divergent from other sequenced hemipteran genomes, showing no detectable synteny [10]. MEAM1 (Middle East-Asia Minor 1, B biotype) and MED (Mediterranean, Q biotype) are the two most widespread and damaging cryptic species in this complex, and MEAM1 has been called a “superbug” because of its highly destructive power in local agriculture [11]. The first global invasion of whiteflies was caused by the transportation of MEAM1 cryptic species on ornamental plants through trade since the late 1980s [12]. MEAM1 whitefly spread to Xinjiang, China, at the end of the 1990s, and its status has been elevated from an accidental pest to a major pest in a short time; these whiteflies remain in some cotton fields in the southern part of Xinjiang at present. Toxicology tests have shown that the resistance of this pest to several pesticides has increased with intensive chemical application in local areas, especially the application of neonicotinoids [13].

Neonicotinoids are commonly used pesticides that target the nicotinic acetylcholine receptor (nAChR) of the central neural system, resulting in excitation and the convulsion of insects, leading to their to death [14,15,16,17]. Imidacloprid, acetamiprid and thiamethoxam are three typical neonicotinoids that were launched in 1991, 1995 and 1998, respectively, and accounted for over 75% of neonicotinoid sales worldwide in 2012 [18]. With the worldwide use of neonicotinoids, whiteflies, especially the MEAM1 and MED cryptic species, have developed different levels of resistance to these insecticides [19,20,21,22,23,24]. Imidacloprid resistance in whiteflies was first reported in 1996, followed by thiamethoxam and acetamiprid resistance [18, 20, 21, 25].

Increased detoxification enzyme activity and the overexpression of corresponding genes are correlated with the resistance of B. tabaci to neonicotine insecticides [26,27,28,29]. CYP6CM1 was the first cytochrome P450 reported to be involved in the imidacloprid resistance of the MEAM1 and MED cryptic species based on high overexpression, and it was also found to participate in the metabolism of thiamethoxam, pymetrozine and pyriproxyfen in B. tabaci [26, 30,31,32]. Further research revealed that a mitogen-activated protein kinase (MAPK) signaling pathway mediates the activation of cAMP response element binding protein (CREB) and that CREB mediates the increase in CYP6CM1 expression via phosphorylation-mediated signal transduction [33]. Resistance to imidacloprid in field populations of B. tabaci is also associated with the increased expression of CYP6CM1 and another cytochrome p450 gene, CYP4C64 [34]. In addition to cytochrome P450 enzymes, carboxylesterases (COEs), NAD-dependent methanol dehydrogenase and glutathione S-transferases (GSTs) have been suggested to be involved in the thiamethoxam resistance of whiteflies [35,36,37,38,39]. RNA sequencing (RNA-seq) and proteomic research have suggested that the overexpression of GST, UDP glucuronosyltransferase (UGT), glucosyl/glucuronosyl transferase and cytochrome P450 proteins plays a role in resistance against thiamethoxam in B. tabaci [40]. CYP303 and CYP6CX3 might be related to the imidacloprid and acetamiprid resistance of B. tabaci, and no evidence of nAChR mutation related to the neonicotinoid resistance of B. tabaci has been found to date [41].

In this study, we established a susceptible strain of the B. tabaci MEAM1 cryptic species by raising the population under controlled indoor conditions without exposure to any pesticide. An imidacloprid-selected strain, acetamiprid-selected strain and thiamethoxam-selected strain were established from the parental susceptible strain by continuous selection with the three neonicotinoids. To determine how this whitefly reacts to long-term selection, the transcriptome profiles of the four strains were sequenced and compared via RNA-seq. It was shown that long-term selection by neonicotinoids led to the upregulation of detoxification-correlated genes belonging to various gene superfamilies (mainly P450 and cuticle proteins), especially in the imidacloprid-selected strain.

Results

Resistance levels of the three neonicotinoid-selected strains

After selection by imidacloprid, acetamiprid and thiamethoxam for 10 generations, the three selected neonicotinoid strains originating from the susceptible strain had developed various levels of resistance against the corresponding neonicotinoids in the following order: IMI strain (43.42-fold) > ACE strain (13.63-fold) > THI strain (2.57-fold) (Table 1).

Table 1 Resistance level of three resistance-selected strains of MEAM1 cryptic species at the 10th generation

Transcriptome data analysis

RNA-seq was utilized to quantify B. tabaci gene expression in the IMI strain, ACE strain, THI strain and SUS strain (control), and a total of 90.86 GB of clean sequence data were obtained. The GC content was 40.07–43.41%, and the Q20 and Q30 values were ≥ 95.81% and ≥ 89.64%, respectively. The average total mapping ratio and average unique mapping ratio were 85.97 and 84.63%, respectively (Table 2).

Table 2 Summary of the RNA-seq data

A total of 16,069 genes were merged and assembled based on all the mapped reads from the four strains using Stringtie (1.3.3b) [42]. The distributions of the expression levels of all the genes were similar among the three selected resistant strains and the susceptible strain (Additional file 1). More than 68% of the genes were overexpressed (FPKM ≥1), and approximately 4.3% of the genes were extremely highly expressed (FPKM > 100) in each strain (Additional file 2).

A total of 1079 new transcripts were detected, and 208 genes corresponded to transcripts in the NCBI RefSeq of B. tabaci (ASM185493v1) (Additional file 3). According to Gene Ontology (GO) [43] classifications, 111 genes were assigned to corresponding GO terms (Additional file 4). The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment results revealed 43 genes assigned to corresponding pathways, including phagosome, endocytosis, lysosome and metabolic pathways (Additional file 5).

Differentially expressed genes between the susceptible and selected resistant strains

Relative to the SUS strain, 584 genes were upregulated and 114 genes were downregulated in the IMI strain, 110 genes were upregulated and 150 genes were downregulated in the ACE strain, 147 genes were upregulated and 50 genes were downregulated in the THI strain. Thirteen genes were upregulated in all 3 selected resistant strains, and 8 of the 13 genes were “protein coding” biotypes. While 4 of these genes were annotated as E3 ubiquitin-protein ligase BRE1-like transcript variant X1, protein D1-like transcript variant X5, maltase 2-like and zingipain-2-like transcript variant X2, the annotation of the remaining genes will require more database information for reference. Two phase II enzymes, UDP-glucuronosyltransferase 2B18-like genes (LOC109030370 and LOC109032577), were obviously downregulated in all three selected resistant strains relative to the SUS strain (Fig. 1) (Additional files 6, 7, 8).

Fig. 1
figure 1

Differentially expressed genes between the whitefly susceptible strain and three resistance selected strains. A: Venn plot of up-regulated DEG numbers between the whitefly susceptible strain and three resistance selected strains. B: Venn plot of down-regulated DEG numbers between the whitefly susceptible strain and three resistance selected strains. Abbreviations: SUS, susceptible strain; IMI, imidacloprid selected strain; ACE, acetamiprid selected strain; THI, thiamethoxam selected strain

The IMI strain exhibited the largest number of upregulated detoxification genes, including 17 cuticle proteins, 8 P450s, 1 carboxylesterase, 1 UDP-glucuronosyltransferase, 1 ABC (ATP-binding cassette), 3 HSPs (heat shock proteins) and 3 potassium channel subfamily K members. Among these genes, venom carboxylesterase-6-like (COE) (LOC109038718) was most highly overexpressed (log2FC = 3.77), followed by cytochrome P450 4C1-like (LOC 109043232) and 10 cuticle proteins (log2FC > 2.00). Cytochrome P450 4C1-like (LOC109043232) and glutathione S-transferase-like transcript variant X1 (LOC109029898) were overexpressed in the ACE strain. Cytochrome P450 4c3-like transcript variant X1 (LOC109043950) and HSP 68-like (LOC109035112) were upregulated in the THI strain (Fig. 2A) (Additional file 9).

Fig. 2
figure 2

Heatmaps of differentially expressed genes correlated with pesticide resistance. Log2foldchange was calculated by DESeq2 method; Log2foldchange value > 1 or < − 1 with padj value < 0.05 was considered as significant. A: Differentially expressed cytochrome P450 genes among the four strains B: Differentially expressed cuticle protein genes among the four strains. C: Differentially expressed acetylcholine receptor genes among the four strains. Abbreviations: SUS, susceptible strain; IMI, imidacloprid selected strain; ACE, acetamiprid selected strain; THI, thiamethoxam selected strain

A total of 17 cuticle protein genes were overexpressed in the IMI strain, including 4 larval cuticle protein A2B-like genes (LOC 109040103, LOC109040136, LOC 109041764, LOC 109040104) and 1 pupal cuticle protein G1A-like (LOC109033087). Ten of the cuticle protein genes were highly expressed (log2FC > 2). No cuticle protein genes were upregulated in the ACE and THI strains. In contrast, 6 cuticle protein genes, including a pupal cuticle protein G1A-like gene (LOC 109033087, log2FC = − 1.19), were downregulated in the ACE strain. Cuticle protein 8-like (LOC109040130) was inhibited in the THI strain (log2FC = − 6.47) and was also the most obviously reduced gene in the strain (Fig. 2B) (Additional file 10).

A total of 13 acetylcholine receptors were annotated in all of the transcriptome profiles, 2 of which were muscarinic acetylcholine receptors (mAChRs), and 11 of which were nicotinic acetylcholine receptors (nAChRs). Nine of the nAChRs belonged to the α subunit group, and the other 2 belonged to the β subunit group. mAChR (LOC109042619) was upregulated in the IMI strain (log2FC = 1.16) and THI strain (log2FC = 1.06) relative to the SUS strain, and nAChR subunit alpha-L1 (LOC109040817, log2FC = 1.10) was upregulated in the IMI strain. No acetylcholine receptors were found to be differentially expressed between the ACE and SUS groups (Fig. 2C) (Additional file 11).

GO and KEGG enrichment analysis

GO enrichment was carried out according to a threshold of padj < 0.05 to analyze the upregulated differentially expressed genes (DEGs) between the SUS strain and the IMI, ACE and THI strains. Compared with the SUS strain, upregulated DEGs in the IMI strain were mainly enriched in biological process (BP) categories including the chitin metabolic, amino sugar metabolic, and glucosamine-containing compound metabolic categories molecular function (MF) categories including the structural constituent of cuticle, chitin binding, and structural molecule activity categories; only 19 genes were enriched in cellular component (CC) categories, among which extracellular region was the predominant in the category. The upregulated genes in the ACE strain were enriched in BP and MF categories, among which DNA integration and O-acyltransferase activity were the two predominant categories, and no genes were enriched in any CC category. The upregulated genes in the THI strain were also enriched in the BP and MF categories, among which DNA integration and acid phosphatase activity were the two predominant categories, and no genes were enriched in any CC category (Fig. 3).

Fig. 3
figure 3

GO enriched terms of DEGs of the three resistance selected strains versus the susceptible strain. The x-axis lists the sub-Go terms under categories of biological process, cellular component, and molecular function. The y-axis is the significant level of GO term enrichment, and a larger value represents a more significant level. A: imidacloprid selected strain versus susceptible strain. B: acetamiprid selected strain versus susceptible strain. C: thiamethoxam selected strain versus susceptible strain. Abbreviations: BP, biological process; CC, cellular component; MF, molecular function

KEGG pathway enrichment analysis was performed using the Cluster Profiler R package. Neuroactive ligand-receptor interaction and the Wnt signaling pathway were the two most significantly DEG-enriched pathways (false discovery rate < 0.05) between the IMI and SUS strains (Fig. 4A). Two overrepresented pathways, lysosome and pentose and glucuronate interconversions, were identified among the DEGs between the ACE strain and SUS strain (Fig. 4B). Fatty acid elongation and protein processing in the endoplasmic reticulum were the two most significantly enriched DEGs between the THI strain and the SUS strain (Fig. 4C). The lysosome pathway was found to be significantly enriched in all 3 selected resistant strains.

Fig. 4
figure 4

KEGG enriched pathways in the resistance selected strains versus the susceptible strain. The pathways were deduced by pairwise comparison of DEGs between the resistance selected strains versus the susceptible strain. Gene ratio represents the ratio of DEG numbers versus the number of genes annotated in the pathway. Larger gene ratio indicates a greater level of enrichment. The padj values ranging from 0 to 1 are P values corrected by multiple hypothesis tests, with lower values indicating greater enrichment. A: imidacloprid selected strain versus susceptible strain. B: acetamiprid selected strain versus susceptible strain. C: thiamethoxam selected strain versus susceptible strain

Validation of DEGs by qRT-PCR

The expression levels of nine genes in the three resistant selected strains were compared by real-time quantitative PCR, and the expression patterns were consistent with the RNA-seq results (Fig. 5). Three kinds of genes were selected to check the consistency of the gene expression trends (selected resistant strains versus susceptible strains) between the RNA-seq and qPCR results. The first group of genes were upregulated (log2FC > 1) in all the selected resistant strains and included genes such as zingipain-2-like transcript variant X2 and maltase 2-like. The second group of genes were upregulated in one or two selected resistant strains relative to the susceptible strain, which included genes such as P450 4C1-like, probable muscarinic acetylcholine receptor gar-2 transcript variant X1, heat shock protein 68-like (LOC109040186), heat shock protein 68-like (LOC109035112), glutathione S-transferase-like transcript variant X1, and cathepsin B-like. Cytochrome P450 4C3-like transcript variant X1 belonged to the third group of genes, which were not upregulated in any of the selected resistant strains relative to the susceptible strain.

Fig. 5
figure 5

Validation of gene expression patterns by RT-qPCR. Means (± SE) were used to determine transcript levels with the 2-ΔΔCt method. One-way analysis of variance (ANOVA) was used to analyze the relative expression levels of selected genes. Panels A - I: cytochrome P450 4C1-like; cytochrome P450 4C3-like transcript variant X1; probable muscarinic acetylcholine receptor gar-2 transcript variant X1; heat shock protein 68-like (LOC109040186); heat shock protein 68-like (LOC109035112); glutathione S-transferase-like transcript variant X1; cathepsin B-like; zingipain-2-like transcript variant X2; maltase 2-like. Columns labeled with different lowercase letters indicate significant differences among the three resistance selected strains. Turkey’s multiple range test was used for pairwise comparison for mean separation (P < 0.05). The RNA-seq results were represented by broken lines and RT-qPCR results were represented by columns

Discussion

Resistance against imidacloprid has become a common and severe problem among whitefly populations worldwide. In this study, the susceptible MEAM1 whitefly population developed the highest resistance level against imidacloprid (RR = 43.42) after 10 generations of selection relative to acetamiprid selection (RR = 13.63) and thiamethoxam selection (RR = 2.57) under laboratory conditions. Whitefly seems to be able to develop a high level of resistance against imidacloprid within several generations, which might lead to the variation in cross-resistance among neonicotinoids. After continuous selection with thiamethoxam for 30 generations (B biotype), the laboratory strain TH-R (resistance ratio = 25.6) developed a higher level of cross-resistance against imidacloprid (cross-resistance ratio = 47.3) than acetamiprid (cross-resistance ratio = 35.8) [35]. Q biotype whiteflies collected in 2000 (Spain) and 2001 (Germany) showed the highest resistance against imidacloprid, and neonicotinoid cross-resistance decreased in the order of imidacloprid > thiamethoxam > acetamiprid [44]. The physiological and molecular mechanisms underlying the ability to rapidly establish resistance against imidacloprid are worth considering.

Physiological research has revealed that pesticide resistance is largely associated with increased detoxification enzyme activity, and the P450 enzyme is crucial for the establishment of resistance against neonicotinoids [35, 44, 45]. Important molecular evidence of neonicotinoid resistance was also obtained, and CYP6CM1 was suggested to be one of the most important genes in imidacloprid resistance in whiteflies [26, 46, 47]. In this study, transcriptomic comparisons between the three neonicotinoid-selected MEAM1 whitefly strains and the susceptible strain revealed that continuous selection by imidacloprid led to the overexpression of a greater number of genes (especially P450 and cuticle protein genes) than acetamiprid selection or thiamethoxam selection. Despite breakthroughs in imidacloprid resistance research in whiteflies, key genes involved in the metabolism of thiamethoxam and acetamiprid remain unclear. In the transcriptome of the IMI strain, eight P450s were obviously upregulated, including CYP6CM1 and four P450 4C1-like genes, while the P450 4C1-like gene (LOC109043232) was the only overexpressed P450 superfamily member in the ACE strain transcriptome, and the P450 4C3-like gene (LOC 109043950) was the only overexpressed P450 in the THI strain transcriptome. Further study of these P450s might provide us with a comprehensive view of the survival strategy and neonicotinoid resistance formation in the whitefly MEAM1 cryptic species.

Insect cuticle is a crucial determinant of insecticide resistance, and cuticle modifications are mainly attributed to the overexpression of cuticular protein genes [48]. In vivo penetration assays using [3 H] imidacloprid revealed a significant reduction in the penetration of the insecticide through the cuticle in a resistant clone of the peach-potato aphid Myzus persicae [49]. In the current study, the large number and high levels of overexpressed cuticle proteins in the IMI strain obviously increased resistance against the pesticide. Cuticle protein 8, together with 12 other cuticle protein genes, was significantly downregulated in a thiamethoxam-resistant strain of the cotton aphid Aphis gossypii Glover [50]; similarly, cuticle protein 8 was downregulated in the THI strain of whiteflies in our study. How the up- or downregulation pattern of cuticle protein genes affects resistance against the three tested neonicotinoids warrants further investigation.

The toxicity of neonicotinoids is related to the binding affinity of neonicotinoids and relative nAChRs to some extent [14, 51, 52]. An RNA-seq analysis of MED whiteflies revealed no nAChR polymorphism potentially related to resistant phenotypes, but nAChR subunits β2 and α7 were upregulated in an acetamiprid-selected strain and imidacloprid-selected strain relative to a susceptible strain [41]. In the current study, nAChRα-L1 was obviously upregulated in the IMI strain and slightly upregulated in the THI strain (log2FC = 0.88), and nAChR α3 (log2FC = 0.93) and nAChR β1 (log2FC = 0.72) were also upregulated in the IMI strain. Unlike the common link between the overexpression of dextoxification enzyme genes and pesticide resistance, decreased expression levels of nAChR subunits β1 and β2 were observed after imidacloprid exposure in Apis cerana cerana [53], and reduced expression of nAChR subunit α2 was correlated with the neonicotinoid resistance phenotype in Musca domestica L. [54]. The whiteflies used in both the research by Illias et al. (2015) and our study included newly emerged adult whiteflies without further pesticide treatment, which showed the upregulation of some nAChRs, while the expression of nAChRs in insects that had experienced recent neonicotinoid treatment was obviously reduced, as mentioned above. The exposure of the larval aphid Acyrthosiphon pisum to imidacloprid and thiamethoxam was shown to result in the up- or downregulation of various nAChRs [55]. It can still be argued that reduced expression of nAChRs is more likely to be correlated with instant neonicotinoid stimuli, and the upregulation of nAChRs is probably a survival strategy implemented by insect pests when experiencing long-term selection by neonicotinoids. A comparison of the expression patterns of nAChRs in selected resistant and susceptible strains of B. tabaci before and after neonicotinoid stimulation might provide a more vivid illustration of pesticide resistance mechanisms.

The muscarinic acetylcholine receptor (mAChR) is a metabotropic G-protein-coupled receptor that initiates continuous intracellular signaling events [56]. Five mAChRs exist in mammals (named m1 – m5), which can be fully activated by the classical agonist muscarine and blocked by the classical mAChR antagonist atropine [57]. GAR-1, GAR-2 and GAR-3 are three new mAChRs in Caenorhabditis elegans that are similar to but pharmacologically distinct from mAChRs m1- m5 [58, 59]. GAR-2 is most similar to GAR-1 and closely related to GAR-3, but GAR-2 was not inhibited at any of three tested muscarinic antagonists (atropine, scopolamine, and pirenzepine) and was found to mainly be expressed in different types of neuronal cells [58]. In our study, mAChR gar-2 was upregulated in all three neonicotinoid-selected strains, and mAChR DM1 was slightly upregulated in the IMI strain. Since no insect mAChR has been reported to participate in the establishment of insecticide resistance, whether the observed overexpression is related to the resistance of the whitefly MEAM1 cryptic species needs further examination.

GSTs and UGTs are thought to be involved in the metabolism of xenobiotics in phase II by increasing the polarity of xenobiotics and facilitating their excretion. GSTd7 has been related to resistance to imidacloprid in the MEAM1 and MED cryptic species [60], and GST14 is thought to be related to the resistance of thiamethoxam of the Q-biotype [61]. UGTs have been connected with resistance against imidacloprid, thiamethoxam, and chlorantraniliprole in Leptinotarsa decemlineata, Aphis gossypii Glover, and Plutella xylostella (L.), respectively [62,63,64,65,66]. In the current study, long-term selection with the three neonicotinoids led to the overexpression of one GST gene (glutathione S-transferase-like variant X1, Sigma class) in the ACE strain and one UGT gene (UGT 2B7-like, 2B class) in the IMI strain. No GSTs were found to be up- or downregulated in the IMI strain or the THI strain. Since no GSTs or UGTs of the same classes as the two differentially expressed GST-like and UGT 2B-like genes have been associated with pesticide resistance in whiteflies or other insect pests, how these two genes are related to imidacloprid resistance or acetamiprid resistance needs further investigation.

Conclusions

Whiteflies of the MEAM1 cryptic species originating from the susceptible parental strain were selected continuously for 10 generations with imidacloprid, acetamiprid or thiamethoxam. The transcriptomes of the four examined strains were established by RNA-seq to investigate the resistance mechanism and survival strategy under the selection pressure imposed by neonicotinoids. The results demonstrated that different resistance levels developed in the whitefly MEAM1 cryptic species in the order of IMI > ACE > THI. In accord with the significantly higher resistance level, more DEGs were found in the IMI strain, mainly P450s and cuticle proteins. DEGs were also found in the ACE strain and THI strain, including P450s, GSTs, nAChRs, mAChRs, HSPs, and COEs. The DEGs between the susceptible strain and the selected resistant strains will provide a basis for further understanding neonicotinoid resistance and cross-resistance in B. tabaci.

Materials and methods

Insects

The susceptible MEAM1 (SUS) strain was collected at Xinjiang Agricultural University in 2010 and reared indoors on cotton plants without any exposure to insecticides thereafter. The identity of the strain was determined by mitochondrial cytochrome oxidase I gene sequencing (mtCOI) [67].

Resistance selection

Adults of the susceptible strain were treated with imidacloprid, acetamiprid and thiamethoxam for 10 generations to establish an imidacloprid-selected (IMI) strain, acetamiprid-selected (ACE) strain and thiamethoxam-selected (THI) strain as follows. The following technical products were used for resistance selection and bioassays: imidacloprid (700 g/kg, Bayer CropScience), acetamiprid (200 g/kg, Nippon Soda), and thiamethoxam (250 g/kg, Syngenta). For resistance selection, the three pesticides were diluted with 0.01% Tween 80 to test the LC50 value against the corresponding strain in each generation, and the population was then separately selected with the three pesticides based on the LC50 values. The LC50 value test was carried out following methods described previously with slight modification [13, 68]. In brief, cotton leaves were cut into discs (38 mm diameter), which were then dipped into pesticide solutions for 20 s and allowed to dry in air. The leaf discs were then transferred to agar (1.7%) plates prepared in advance. Female adults that emerged after 2–3 days were collected and transferred to the treated cotton leaf discs and kept for 48 h; 30–40 individuals were used for each of the treatments, which were replicated 3 times. The LC50 values were calculated with the PoloPlus Software Program (Leora Software, Berkeley, CA). The three pesticides were then diluted to the proper concentration based on the LC50 value, sprayed on cotton leaves and allowed to dry. Adult whiteflies were then transferred to the treated cotton plants and selected for one generation. The procedure was repeated in each generation until the 10th generation for further testing. All strains were kept under controlled environmental conditions with a temperature of 25 ± 1 °C, relative humidity of 70–880% and a photoperiod of 16 L/8 D.

RNA extraction and quality control

Total RNA was extracted from susceptible and resistant MEAM1 whitefly adults (mixture of 50 male individuals and 50 female individuals). Three biological replicates were included for each of the 4 strains, and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library preparation for sequencing

A total amount of 1 μg RNA per sample was used for RNA sample preparation. Messenger RNA (mRNA) was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in First Strand Synthesis Reaction Buffer (5X). First-strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase (RNase H). Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exon-nuclease/polymerase activities. After the adenylation of the 3′ ends of the DNA fragments, adaptors with hairpin loop structures were ligated to prepare for hybridization. To preferentially select cDNA fragments of 370–420 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and an index (X) Primer. Finally, the PCR products were purified (AMPure XP system), and library quality was assessed on an Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150 bp paired-end reads were generated [69].

Transcriptome data analysis

Raw data (raw reads) in FASTQ format were first processed with in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N sequences and low-quality reads from the raw data. At the same time, the Q20, Q30, and GC contents of the clean data were calculated. All downstream analyses were based on clean data with high quality. Reference genome and gene model annotation files were downloaded from the B. tabaci reference genome (assembly ASM185493v1, http://www.whiteflygenomics.org/cgi-bin/bta/index.cgi) [10]. The index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5 [70]. We selected Hisat2 as the mapping tool because Hisat2 can generate a database of splice junctions based on the gene model annotation file and thus provides better mapping results than other non-splice mapping tools.

Differential gene expression analysis and functional enrichment

Feature Counts v1.5.0-p3 [71] was used to count the read numbers mapped to each gene. Then, the FPKM value (the expected number of fragments per kilobase of transcript sequence per million base pairs sequenced) of each gene was calculated based on the length of the gene and read count mapped to the gene. The differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.20.0) [72, 73]. A corrected P-value of 0.05 and a log2 foldchange of ±1 were set as the thresholds for significant differential expression.

GO annotation and GO/KEGG enrichment

The GO [43] enrichment analysis of DEGs was implemented with the cluster Profiler R package, in which gene length bias was corrected. GO terms with corrected P values < 0.05 were considered significantly enriched by DEGs. The Cluster Profiler R package was used to test the statistical enrichment of DEGs in KEGG pathways (http://www.genome.jp/kegg/) [74, 75]. Gene set enrichment analysis (GSEA) is a computational approach for determining whether a predefined gene set shows a significant consistent difference between two biological states. The genes were ranked according to the degree of differential expression in the two samples, and the predefined gene set was then tested to determine whether the enriched genes were at the top or bottom of the list. GSEA can include subtle expression changes. We used the local version of the GSEA tool (http://www.broadinstitute.org/gsea/index.jsp), and the GO and KEGG datasets were used for GSEA independently.

Validation of DEGs by quantitative real-time PCR (qRT-PCR)

Primers for nine genes were designed and validated by qRT-PCR. The primers (Additional file 12) were designed with Primer3plus (http://www.primer3plus.com/cgi-bin/dev/primer3plus.cgi). First-strand cDNA was synthesized using Total RNA Extractor (TRIzol) (Sangon Biotech, China) and an M-MuLV First Strand cDNA Synthesis Kit (Sangon Biotech, China). Total RNA (0.5 μg) was reverse-transcribed, and qRT-PCR was conducted with a 7500 Fast Real-time PCR System (ABI) in a 20 μl reaction volume using 2X SG Fast qPCR Master Mix (Low Rox) (Sangon Biotech, China) according to the following protocol: 3 min of activation at 95 °C, followed by 40 cycles of 3 s at 95 °C and 25 s at 60 °C. Relative changes in gene expression were assessed using the 2−ΔΔCt method [76], and elongation factor 1 alpha (EF1α) was used as a reference gene [77, 78]. Samples were assessed in triplicate. Data were analyzed by one-way ANOVA, followed by Tukey’s multiple comparisons and analysis with SPSS v. 19.0 (SPSS, Chicago, IL, USA). For ANOVA, the data were transformed for homogeneity of variance tests. Differences were considered statistically significant when P < 0.05.

Availability of data and materials

The RNA-Seq data are available from the NCBI/GenBank BioProject database under the identifier PRJNA752712.

References

  1. Kunz D, Tay WT, Court LN, Elfekih S, Gordon KHJ, Evans GA, et al. Draft mitochondrial DNA genome of a 1920 Barbados cryptic Bemisia tabaci ‘New World’ species (Hemiptera: Aleyrodidae). Mitochondrial DNA B. 2019;4(1):1183–4.

    Article  Google Scholar 

  2. Chu D, Qu W-M, Guo L. Invasion genetics of alien insect pests in China: research progress and future prospects. J Integr Agric. 2019;18(4):748–57.

    Article  Google Scholar 

  3. Vyskocilova S, Tay WT, van Brunschot S, Seal S, Colvin J. An integrative approach to discovering cryptic species within the Bemisia tabaci whitefly species complex. Sci Rep. 2018;8(1):10886–98.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Mugerwa H, Seal S, Wang HL, Patel MV, Kabaalu R, Omongo CA, et al. African ancestry of New World, Bemisia tabaci-whitefly species. Sci Rep. 2018;8(1):2734–44.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Liu SS, De Barro PJ, Xu J, Luan JB, Zang LS, Ruan YM, et al. Asymmetric mating interactions drive widespread invasion and displacement in a whitefly. Science. 2007;318(5857):1769–72.

    Article  CAS  PubMed  Google Scholar 

  6. Boykin LM, Bell CD, Evans G, Small I, De Barro PJ. Is agriculture driving the diversification of the Bemisia tabaci species complex (Hemiptera: Sternorrhyncha: Aleyrodidae)?: dating, diversification and biogeographic evidence revealed. BMC Evol Biol. 2013;13(228):1–10.

    Google Scholar 

  7. Basit M. Status of insecticide resistance in Bemisia tabaci: resistance, cross-resistance, stability of resistance, genetics and fitness costs. Phytoparasitica. 2019;47(2):207–25.

    Article  CAS  Google Scholar 

  8. Liu B, Yan F, Chu D, Pan H, Jiao X, Xie W, et al. Difference in feeding behaviors of two invasive whiteflies on host plants with different suitability: implication for competitive displacement. Int J Biol Sci. 2012;8(5):697–706.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Wosula EN, Chen W, Fei Z, Legg JP: Unravelling the genetic diversity among cassava Bemisia tabaci whiteflies using NextRAD sequencing. Genome Biol Evol 2017; (No.11):2958–2973.

  10. Chen W, Hasegawa DK, Kaur N, Kliot A, Pinheiro PV, Luan J, et al. The draft genome of whitefly Bemisia tabaci MEAM1, a global crop pest, provides novel insights into virus transmission, host adaptation, and insecticide resistance. BMC Biol. 2016;14(1):110.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Barinaga M. Is devastating whitefly invader really a new species? Science. 1993;259(5091):30.

    Article  CAS  PubMed  Google Scholar 

  12. De Barro PJ, Liu SS, Boykin LM, Dinsdale AB. Bemisia tabaci: a statement of species status. Annu Rev Entomol. 2011;56:1–19.

    Article  PubMed  Google Scholar 

  13. Ma DY, Gorman K, Devine G, Luo WC, Denholm I. The biotype and insecticide-resistance status of whiteflies, Bemisia tabaci (Hemiptera: Aleyrodidae), invading cropping systems in Xinjiang Uygur autonomous region, northwestern China. Crop Prot. 2007;26(4):612–7.

    Article  CAS  Google Scholar 

  14. Tomizawa M, Casida JE. Selective toxicity of neonicotinoids attributable to specificity of insect and mammalian nicotinic receptors. Annu Rev Entomol. 2003;48(1):339–64.

    Article  CAS  PubMed  Google Scholar 

  15. Casida JE, Durkin KA. Neuroactive insecticides: targets, selectivity, resistance, and secondary effects. Annu Rev Entomol. 2013;58(1):99–117.

    Article  CAS  PubMed  Google Scholar 

  16. Simon-Delso N, Amaral-Rogers V, Belzunces LP, Bonmatin JM, Chagnon M, Downs C, et al. Systemic insecticides (neonicotinoids and fipronil): trends, uses, mode of action and metabolites. Environ Sci Pollut Res. 2015;22(1):5–34.

    Article  CAS  Google Scholar 

  17. Casida JE. Neonicotinoids and other insect nicotinic receptor competitive modulators: Progress and prospects. Annu Rev Entomol. 2018;63:125–44.

    Article  CAS  PubMed  Google Scholar 

  18. Bass C, Denholm I, Williamson MS, Nauen R. The global status of insect resistance to neonicotinoid insecticides. Pestic Biochem Physiol. 2015;121:78–87.

    Article  CAS  PubMed  Google Scholar 

  19. Prabhaker N, Toscano NC, Castle SJ, Henneberry TJ. Selection for imidacloprid resistance in silverleaf whiteflies from the Imperial Valley and development of a hydroponic bioassay for resistance monitoring. Pest Manag Sci. 1997;51(4):419–28.

    Article  CAS  Google Scholar 

  20. Elbert A, Nauen R. Resistance of Bemisia tabaci (Homoptera: Aleyrodidae) to insecticides in southern Spain with special reference to neonicotinoids. Pest Manag Sci. 2000;56(1):60–4.

    Article  CAS  Google Scholar 

  21. Horowitz AR, Kontsedalov S, Ishaaya I. Dynamics of resistance to the neonicotinoids acetamiprid and thiamethoxam in Bemisia tabaci (Homoptera: Aleyrodidae). J Econ Entomol. 2004;97(6):2051–6.

    Article  CAS  PubMed  Google Scholar 

  22. Schuster DJ, Mann RS, Toapanta M, Cordero R, Thompson S, Cyman S, et al. Monitoring neonicotinoid resistance in biotype B of Bemisia tabaci in Florida. Pest Manag Sci. 2010;66(2):186–95.

    Article  CAS  PubMed  Google Scholar 

  23. Wang Z, Yan H, Yang Y, Wu Y. Biotype and insecticide resistance status of the whitefly Bemisia tabaci from China. Pest Manag Sci. 2010;66(12):1360–6.

    Article  CAS  PubMed  Google Scholar 

  24. Vassiliou V, Emmanouilidou M, Perrakis A, Morou E, Vontas J, Tsagkarakou A, et al. Insecticide resistance in Bemisia tabaci from Cyprus. Insect Sci. 2011;18(1):30–9.

    Article  CAS  Google Scholar 

  25. Cahill M, Gorman K, Day S, Denholm I, Elbert A, Nauen R. Baseline determination and detection of resistance to imidacloprid in Bemisia tabaci (Homoptera: Aleyrodidae). Bull Entomol Res. 1996;86(4):343–9.

    Article  CAS  Google Scholar 

  26. Karunker I, Benting J, Lueke B, Ponge T, Nauen R, Roditakis E, et al. 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. 2008;38(6):634–44.

    Article  CAS  PubMed  Google Scholar 

  27. Basit M, Saeed S, Saleem MA, Denholm I, Shah M. Detection of resistance, cross-resistance, and stability of resistance to new chemistry insecticides in Bemisia tabaci (Homoptera: Aleyrodidae). J Econ Entomol. 2013;106(3):1414–22.

    Article  CAS  PubMed  Google Scholar 

  28. Guo L, Xie W, Wang S, Wu Q, Li R, Yang N, et al. Detoxification enzymes of Bemisia tabaci B and Q: biochemical characteristics and gene expression profiles. Pest Manag Sci. 2014;70(10):1588–94.

    Article  CAS  PubMed  Google Scholar 

  29. Zhou C-s, Cao Q, Li G-Z, Ma D-Y. Role of several cytochrome P450s in the resistance and cross-resistance against imidacloprid and acetamiprid of Bemisia tabaci (Hemiptera: Aleyrodidae) MEAM1 cryptic species in Xinjiang, China. Pestic Biochem Physiol. 2020;163:209–15.

    Article  CAS  PubMed  Google Scholar 

  30. Nauen R, Vontas J, Kaussmann M, Wolfel K. Pymetrozine is hydroxylated by CYP6CM1, a cytochrome P450 conferring neonicotinoid resistance in Bemisia tabaci. Pest Manag Sci. 2013;69(4):457–61.

    Article  CAS  PubMed  Google Scholar 

  31. Nauen R, Wolfel K, Lueke B, Myridakis A, Tsakireli D, Roditakis E, et al. Development of a lateral flow test to detect metabolic resistance in Bemisia tabaci mediated by CYP6CM1, a cytochrome P450 with broad spectrum catalytic efficiency. Pestic Biochem Physiol. 2015;121:3–11.

    Article  CAS  PubMed  Google Scholar 

  32. Yang N, Xie W, Jones CM, Bass C, Jiao X, Yang X, et al. Transcriptome profiling of the whitefly Bemisia tabaci reveals stage-specific gene expression signatures for thiamethoxam resistance. Insect Mol Biol. 2013;22(5):485–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Yang X, Deng S, Wei X, Yang J, Zhao Q, Yin C, et al. MAPK-directed activation of the whitefly transcription factor CREB leads to P450-mediated imidacloprid resistance. Proc Natl Acad Sci U S A. 2020;117(19):10246–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Yang X, Xie W, Wang S-L, Wu Q-J, Pan H-P, Li R-M, et al. Two cytochrome P450 genes are involved in imidacloprid resistance in field populations of the whitefly, Bemisia tabaci, in China. Pestic Biochem Physiol. 2013;107(3):343–50.

    Article  CAS  PubMed  Google Scholar 

  35. Feng Y, Wu Q, Wang S, Chang X, Xie W, Xu B, et al. Cross-resistance study and biochemical mechanisms of thiamethoxam resistance in B-biotype Bemisia tabaci (Hemiptera: Aleyrodidae). Pest Manag Sci. 2010;66(3):313–8.

    Article  CAS  PubMed  Google Scholar 

  36. Kang S, Lee HJ, Kim YH, Kwon DH, Oh JH, Kim BJ, et al. Proteomics-based identification and characterization of biotype-specific carboxylesterase 2 putatively associated with insecticide resistance in Bemisia tabaci. J Asia-Pacif Entomol. 2012;15(3):389–96.

    Article  CAS  Google Scholar 

  37. Xie W, Yang X, Wang SI, Wu QJ, Yang NN, Li RM, et al. Gene expression profiling in the thiamethoxam resistant and susceptible B-biotype sweetpotato whitefly, Bemisia tabaci. J Insect Sci. 2012;12(46):46–56.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Xia J, Xu H, Yang Z, Pan H, Yang X, Guo Z, et al. Genome-wide analysis of carboxylesterases (COEs) in the whitefly, Bemisia tabaci (Gennadius). Int J Mol Sci. 2019;20(20):4973–88.

    Article  CAS  PubMed Central  Google Scholar 

  39. Eakteiman G, Moses-Koch R, Moshitzky P, Mestre-Rincon N, Vassão DG, Luck K, et al. Targeting detoxification genes by phloem-mediated RNAi: a new approach for controlling phloem-feeding insect pests. Insect Biochem Mol Biol. 2018;100:10–21.

    Article  CAS  PubMed  Google Scholar 

  40. Yang N, Xie W, Yang X, Wang S, Wu Q, Li R, et al. Transcriptomic and proteomic responses of sweetpotato whitefly, Bemisia tabaci, to thiamethoxam. PLoS One. 2013;8(5):e61820.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Ilias A, Lagnel J, Kapantaidaki DE, Roditakis E, Tsigenopoulos CS, Vontas J, et al. Transcription analysis of neonicotinoid resistance in Mediterranean (MED) populations of B. tabaci reveal novel cytochrome P450s, but no nAChR mutations associated with the phenotype. BMC Genomics. 2015;16:1–23.

    Article  Google Scholar 

  42. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 2010;11(R14):1–12.

    Google Scholar 

  44. Rauch N, Nauen R. Identification of biochemical markers linked to neonicotinoid cross resistance in Bemisia tabaci (Hemiptera: Aleyrodidae). Arch Insect Biochem Physiol. 2003;54(4):165–76.

    Article  CAS  PubMed  Google Scholar 

  45. Wang Z, Yao M, Wu Y. Cross-resistance, inheritance and biochemical mechanisms of imidacloprid resistance in B-biotype Bemisia tabaci. Pest Manag Sci. 2009;65(11):1189–94.

    Article  CAS  PubMed  Google Scholar 

  46. Karunker I, Morou E, Nikou D, Nauen R, Sertchook R, Stevenson BJ, et al. Structural model and functional characterization of the Bemisia tabaci CYP6CM1vQ, a cytochrome P450 associated with high levels of imidacloprid resistance. Insect Biochem Mol Biol. 2009;39(10):697–706.

    Article  CAS  PubMed  Google Scholar 

  47. Roditakis E, Morou E, Tsagkarakou A, Riga M, Nauen R, Paine M, et al. Assessment of the Bemisia tabaci CYP6CM1vQ transcript and protein levels in laboratory and field-derived imidacloprid-resistant insects and cross-metabolism potential of the recombinant enzyme. Insect Sci. 2011;18(1):23–9.

    Article  CAS  Google Scholar 

  48. Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74.

    Article  PubMed  Google Scholar 

  49. Alin M. Puinean, Stephen P. Foster, Linda Oliphant, Ian Denholm, Linda M. Field, Neil S. Millar, Martin S. Williamson, Bass C: Amplification of a Cytochrome P450 Gene Is Associated with Resistance to Neonicotinoid Insecticides in the Aphid Myzus persicae. PLoS Genet 2010, 6(No.6):e1000999.

  50. Pan Y, Peng T, Gao X, Zhang L, Yang C, Xi J, et al. Transcriptomic comparison of thiamethoxam-resistance adaptation in resistant and susceptible strains of Aphis gossypii glover. Comp Biochem Physiol Part D Genomics Proteomics. 2015;13:10–5.

    Article  CAS  PubMed  Google Scholar 

  51. Tomizawa M, Casida JE. Neonicotinoid insecticide toxicology: mechanisms of selective action. Annu Rev Pharmacol Toxicol. 2004;45(1):247–68.

    Article  Google Scholar 

  52. Casida JE. Neonicotinoid metabolism: compounds, substituents, pathways, enzymes, organisms, and relevance. J Agric Food Chem. 2011;59(7):2923–31.

    Article  CAS  PubMed  Google Scholar 

  53. Yu X, Wang M, Kang M, Liu L, Guo X, Xu B. Molecular cloning and characterization of two nicotinic acetylcholine receptor β subunit genes from Apis cerana cerana. Arch Insect Biochem Physiol. 2011;77(4):163–78.

    Article  CAS  PubMed  Google Scholar 

  54. Markussen MD, Kristensen M. Low expression of nicotinic acetylcholine receptor subunit Mdα2 in neonicotinoid-resistant strains of Musca domestica L. Pest Manag Sci. 2010;66(11):1257–62.

    Article  CAS  PubMed  Google Scholar 

  55. Taillebois E, Beloula A, Quinchard S, Jaubert-Possamai S, Daguin A, Servent D, et al. Neonicotinoid binding, toxicity and expression of nicotinic acetylcholine receptor subunits in the aphid Acyrthosiphon pisum. PLoS One. 2014;9(5):e96669.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Honda H, Tomizawa M, Casida JE. Insect muscarinic acetylcholine receptor: pharmacological and toxicological profiles of antagonists and agonists. J Agric Food Chem. 2007;55(6):2276–81.

    Article  CAS  PubMed  Google Scholar 

  57. Collin C, Hauser F, Valdivia G, Li S, Reisenberger J, Carlsen EMM, et al. Two types of muscarinic acetylcholine receptors in Drosophila and other arthropods. Cell Mol Life Sci. 2013;70:3231–42.

    Article  CAS  PubMed  Google Scholar 

  58. Lee YS, Park YS, Nam S, Suh SJ, Lee J, Kaang BK, et al. Characterization of GAR-2, a novel G protein-linked acetylcholine receptor from Caenorhabditis elegans. J Neurochem. 2000;75:1800–9.

    Article  CAS  PubMed  Google Scholar 

  59. Lee YS, Park YS, Chang DJ, Hwang JM, Min CK, Kaang BK, Cho NJ: Cloning and expression of a G protein-linked acetylcholine receptor from Caenorhabditis elegans. J Neurochem 1999, 72(No.1):58-65.

  60. He C, Xie W, Yang X, Wang SL, Wu QJ, Zhang YJ. Identification of glutathione S-transferases in Bemisia tabaci (Hemiptera: Aleyrodidae) and evidence that GSTd7 helps explain the difference in insecticide susceptibility between B.tabaci Middle East-minor Asia 1 and Mediterranean. Insect Mol Biol. 2018;27(1):22–35.

    Article  CAS  PubMed  Google Scholar 

  61. Yang X, He C, Xie W, Liu Y, Xia J, Yang Z, et al. Glutathione S-transferases are involved in thiamethoxam resistance in the field whitefly Bemisia tabaci Q (Hemiptera: Aleyrodidae). Pestic Biochem Physiol. 2016;134:73–8.

    Article  CAS  PubMed  Google Scholar 

  62. Kaplanoglu E, Chapman P, Scott IM, Donly C. Overexpression of a cytochrome P450 and a UDP-glycosyltransferase is associated with imidacloprid resistance in the Colorado potato beetle, Leptinotarsa decemlineata. Sci Rep. 2017;7(1):1762–71.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Li X, Zhu B, Gao X, Liang P: Over-expression of UDP-glycosyltransferase gene UGT2B17 is involved in chlorantraniliprole resistance in Plutella xylostella (L.). Pest Manag Sci 2017, 73(No.7):1402–1409.

  64. Tian FJ, Wang ZB, Li CF, Liu JL, Zeng XN. UDP-glycosyltransferases are involved in imidacloprid resistance in the Asian citrus psyllid, Diaphorina citri (Hemiptera: Lividae). Pestic Biochem Physiol. 2019;154:23–31.

    Article  CAS  PubMed  Google Scholar 

  65. Pan Y, Tian F, Wei X, Wu Y, Gao X, Xi J, et al. Thiamethoxam resistance in Aphis gossypii glover relies on multiple UDP-glucuronosyltransferases. Front Physiol. 2018;9(322):1–9.

    Google Scholar 

  66. Xuewei C, Jin X, Qingli S, Dunlun S, Xiwu G. UDP-glucosyltransferases potentially contribute to imidacloprid resistance in Aphis gossypii glover based on transcriptomic and proteomic analyses. Pestic Biochem Physiol. 2019;159:98–106.

    Article  Google Scholar 

  67. Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann Entomol Soc Am. 1994;87(6):651–701.

    Article  CAS  Google Scholar 

  68. Cahill M, Byrne FJ, Gorman K, Denholm I, Devonshire AL. Pyrethroid and organophosphate resistance in the tobacco whitefly Bemisia tabaci (Homoptera: Aleyrodidae). Bull Entomol Res. 1995;85(2):181–7.

    Article  CAS  Google Scholar 

  69. Fang S-M, Hu B-L, Zhou Q-Z, Yu Q-Y, Zhang Z. Comparative analysis of the silk gland transcriptomes between the domestic and wild silkworms. BMC Genomics. 2015;16(1):60 (1–12).

    Article  PubMed  PubMed Central  Google Scholar 

  70. Kim D, Langmead B, Salzberg SL: HISAT: A fast spliced aligner with low memory requirements (article). Nat Methods 2015, 12(No.4):357-360.

  71. Liao Y, Smyth G, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    Article  CAS  PubMed  Google Scholar 

  72. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;10(R106):1–28.

    Google Scholar 

  73. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(550):1–21.

    Google Scholar 

  74. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔ C T method. Methods 2001, 25(No.4):402-408.

  77. Li R, Xie W, Wang S, Wu Q, Yang N, Yang X, et al. Reference gene selection for qRT-PCR analysis in the sweetpotato whitefly, Bemisia tabaci (Hemiptera: Aleyrodidae). PLoS One. 2013;8(1):e53006.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  78. Collins C, Patel MV, Colvin J, Bailey D, Seal S. Identification and evaluation of suitable reference genes for gene expression studies in the whitefly Bemisia tabaci (Asia I) by reverse transcription quantitative realtime PCR. J Insect Sci. 2014;14(5):1–25.

    Article  Google Scholar 

Download references

Acknowledgments

We thank an anonymous reviewer for the constructive scientific review of this manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 31760524).

Author information

Authors and Affiliations

Authors

Contributions

DYM financially supported and conceptualized the study, critically revised the manuscript and gave final approval for submission. CSZ and HHL designed the experiments. CSZ, HHL, XHG, QC and RXYZ collected samples and performed the experiments. CSZ and HHL analyzed the data. CSZ drafted and revised the manuscript. All authors reviewed and approved the final manuscript.

Corresponding author

Correspondence to De Ying Ma.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhou, C.S., Lv, H.H., Guo, X.H. et al. Transcriptional analysis of Bemisia tabaci MEAM1 cryptic species under the selection pressure of neonicotinoids imidacloprid, acetamiprid and thiamethoxam. BMC Genomics 23, 15 (2022). https://doi.org/10.1186/s12864-021-08241-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-021-08241-6

Keywords