- Open Access
Transcriptional analysis of Bemisia tabaci MEAM1 cryptic species under the selection pressure of neonicotinoids imidacloprid, acetamiprid and thiamethoxam
BMC Genomics volume 23, Article number: 15 (2022)
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.
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).
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.
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 . 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 . The first global invasion of whiteflies was caused by the transportation of MEAM1 cryptic species on ornamental plants through trade since the late 1980s . 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 .
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 . 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 . Resistance to imidacloprid in field populations of B. tabaci is also associated with the increased expression of CYP6CM1 and another cytochrome p450 gene, CYP4C64 . 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 . 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 .
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.
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).
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).
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) . 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)  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).
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).
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).
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.
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.
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) . 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 . 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 . 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 . 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 ; 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 . 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 , and reduced expression of nAChR subunit α2 was correlated with the neonicotinoid resistance phenotype in Musca domestica L. . 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 . 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 . 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 . 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 . 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 , and GST14 is thought to be related to the resistance of thiamethoxam of the Q-biotype . 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.
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
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) .
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 .
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) . 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 . 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  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  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 , 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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Barinaga M. Is devastating whitefly invader really a new species? Science. 1993;259(5091):30.
De Barro PJ, Liu SS, Boykin LM, Dinsdale AB. Bemisia tabaci: a statement of species status. Annu Rev Entomol. 2011;56:1–19.
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.
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.
Casida JE, Durkin KA. Neuroactive insecticides: targets, selectivity, resistance, and secondary effects. Annu Rev Entomol. 2013;58(1):99–117.
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.
Casida JE. Neonicotinoids and other insect nicotinic receptor competitive modulators: Progress and prospects. Annu Rev Entomol. 2018;63:125–44.
Bass C, Denholm I, Williamson MS, Nauen R. The global status of insect resistance to neonicotinoid insecticides. Pestic Biochem Physiol. 2015;121:78–87.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Balabanidou V, Grigoraki L, Vontas J. Insect cuticle: a critical determinant of insecticide resistance. Curr Opin Insect Sci. 2018;27:68–74.
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.
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.
Tomizawa M, Casida JE. Neonicotinoid insecticide toxicology: mechanisms of selective action. Annu Rev Pharmacol Toxicol. 2004;45(1):247–68.
Casida JE. Neonicotinoid metabolism: compounds, substituents, pathways, enzymes, organisms, and relevance. J Agric Food Chem. 2011;59(7):2923–31.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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).
Kim D, Langmead B, Salzberg SL: HISAT: A fast spliced aligner with low memory requirements (article). Nat Methods 2015, 12(No.4):357-360.
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.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;10(R106):1–28.
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.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.
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.
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.
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.
We thank an anonymous reviewer for the constructive scientific review of this manuscript.
This work was supported by the National Natural Science Foundation of China (Grant No. 31760524).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no conflicts of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
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
- Bemisia tabaci transcriptome