Abstract

   Myzus persicae is a serious and widespread agricultural pest, against
   which, imidacloprid remains an effective control measure. However,
   recent reports indicate that this aphid has evolved and developed
   resistance to imidacloprid. This study aimed to elucidate the
   underlying mechanisms and genetic basis of this resistance by
   conducting comparative transcriptomics studies on both
   imidacloprid-resistant (IR) and imidacloprid-susceptible (IS) M.
   persicae. The comparative analysis identified 252 differentially
   expressed genes (DEGs) among the IR and IS M. persicae transcriptomes.
   These candidate genes included 160 and 92 genes that were down- and
   up-regulated, respectively, in the imidacloprid-resistant strain. Using
   functional classification in the GO and KEGG databases, 187 DEGs were
   assigned to 303 functional subcategories and 100 DEGs were classified
   into 45 pathway groups. Moreover, several genes were associated with
   known insecticide targets, cuticle, metabolic processes, and oxidative
   phosphorylation. Quantitative real-time PCR of 10 DEGs confirmed the
   trends observed in the RNA sequencing expression profiles. These
   findings provide a valuable basis for further investigation into the
   complicated mechanisms of imidacloprid resistance in M. persicae.

   Subject terms: Transcriptomics, Entomology

Introduction

   The green peach aphid, Myzus persicae (Sulzer), is a widely distributed
   agricultural pest that has been reported to inflict serious damage on
   more than 400 plant species, both directly, through phloem feeding, and
   indirectly, through the transmission of viruses^[28]1. Over the past
   several decades, the pest has typically been controlled using synthetic
   insecticides. However, the excessive use of such insecticides has
   promoted the development of resistance of M. persicae to many chemical
   products, including organophosphates, carbamates and pyrethroids^[29]2.
   In contrast, neonicotinoid insecticides, which exhibit high binding
   affinity to insect nicotinic acetylcholine receptors (nAChRs), remain
   an effective measure for controlling M. persicae, owing to their
   efficacy, long-lasting effects, and harmlessness to mammals^[30]3. In
   fact, the first neonicotinoid insecticide, imidacloprid, is still the
   main product used to control both sucking and biting insect pests and
   is the world’s most popular insecticide^[31]4,[32]5. However, as in
   other types of insecticides, the failure to incorporate insecticide
   resistance management strategies can increase resistance levels in
   target pest populations; imidacloprid resistance in M. persicae has now
   been reported in the USA, Europe, China, and Japan^[33]6–[34]8.

   The insecticides used in pest management can be considered as
   environmental stress factors to insect populations. Therefore, it is a
   common adaptive strategy for pests to develop insecticide
   resistance^[35]9. The main mechanisms of neonicotinoid resistance
   include reduced target-site sensitivity and enhanced metabolic
   detoxification^[36]10. A field-evolved instance of imidacloprid
   resistance in M. persicae was associated with a single mutation (R81T)
   in the loop D region of the nAChR β1 subunit, which reduced the binding
   affinity of nAChR for imidacloprid. This is the first example of
   field-evolved resistance to imidacloprid that is mediated via a
   target-site mechanism^[37]4. However, nAChR mutations have also been
   reported to play dominant roles in imidacloprid resistance in other
   insects. In Nilaparvata lugens, for example, a target site mutation
   (Y151S) in the α1 and α3 subunits of nAChRs appears to be responsible
   for high-level imidacloprid resistance^[38]11, and in Musca domestica,
   the reduced expression of the α2 subunit of nAChRs has been associated
   with imidacloprid resistance^[39]12.

   In addition, a number of studies have associated high expression levels
   of cytochrome P450 genes with neonicotinoids resistance in
   insects^[40]4. For example, the overexpression of CYP6CY3, CYP6G1,
   CYP6CM1 and CYP4C64 has been associated with imidacloprid resistance in
   M. persicae, Drosophila melanogaster and Bemisia
   tabaci^[41]10,[42]13,[43]14. An increase in the activities of
   detoxification enzymes, such as glutathione S-transferases (GSTs) and
   carboxylesterases, is also known to be associated with imidacloprid
   resistance in aphids^[44]15. Although reduced target-site sensitivity
   and enhanced metabolic detoxification are known to contribute to
   imidacloprid resistance, it is also possible that both resistance
   mechanisms and adaptation strategies are complex processes that involve
   an array of metabolic and genetic factors and that such complex
   processes are responsible for the development of imidacloprid
   resistance in M. persicae.

   Global surveys of transcriptional changes in insecticide-treated
   insects could help elucidate the metabolic and regulatory mechanisms
   that underlie insecticide resistance. Currently, powerful
   next-generation RNA sequencing (RNA-Seq) technology can be used to
   determine the gene expression profiles and, thereby, helping to
   elucidate the development of insecticide resistance^[45]16–[46]18. In
   this study, high-throughput RNA-Seq was used to determine the
   transcriptome profiles of imidacloprid-resistant (IR) and
   imidacloprid-susceptible (IS) M. persicae, with a focus on genes that
   could provide insight into the mechanisms of physiological adaptation
   of insects to imidacloprid stress.

Results

RNA sequencing

   Samples of IS and IR M. persicae were subjected to Illumina-based
   RNA-Seq, with three replicates for each strain (IS1, IS2, IS3, IR1,
   IR2, and IR3). The Illumina sequencing data are shown in Table [47]1.
   After filtering, 24,205,298, 24,416,009, and 26,269,490 clean reads
   were obtained from the IS strain and 28,656,579, 23,803,901, and
   23,359,309 clean reads were obtained from the IR strain. The Q20
   percentage and GC content of clean reads in the digital gene expression
   (DGE) libraries ranged from 96.78% to 98.25% and from 39.71% to 40.99%,
   respectively, and a mean of 92.19% of clean reads was mapped to the M.
   persicae genome database.

Table 1.

   Summary of Illumina RNA-sequencing data.
   Samples Raw reads Clean reads Q20 GC content Total mapped Mapped ratio
   Notes
   IS1 24,688,350 24,205,298 96.78% 40.55% 44,350,857 91.61% Replicate 1
   IS2 24,914,645 24,416,009 97.67% 39.80% 43,792,242 89.68% Replicate 2
   IS3 26,805,241 26,269,490 98.25% 40.03% 47,999,205 91.36% Replicate 3
   IR1 29,072,436 28,656,579 98.23% 39.55% 53,969,698 94.17% Replicate 1
   IR2 24,481,332 23,803,901 97.66% 39.71% 44,727,119 93.95% Replicate 2
   IR3 24,443,951 23,359,309 97.62% 40.96% 43,139,976 92.34% Replicate 3
   [48]Open in a new tab

Differentially expressed genes

   Using the DEG-Seq. 2R package, a total of 252 imidacloprid-responsive
   DEGs were identified in the IS and IR M. persicae transcriptomes; these
   included 160 and 92 genes that were down- and up-regulated,
   respectively, in the IR M. persicae (Fig. [49]1, Supplementary
   Table [50]S1).

Figure 1.

   [51]Figure 1
   [52]Open in a new tab

   Differentially expressed gene (DEG) analysis of imidacloprid-resistant
   and imidacloprid-susceptible Myzus Persicae. (a) Volcano plot of DEGs.
   Dots represent individual genes. Red dots represent up-regulated genes,
   and green dots down-regulated genes. Blue dots indicate genes that are
   not differentially expressed. (b) Heatmap analysis of hierarchical
   clustering of DEGs. Red and blue indicate high and low expression in
   the imidacloprid-resistant strain, respectively.

Functional annotation and classification

   Gene Ontology (GO) indicated that 186 (73.8%) DEGs were assigned to 303
   subcategories, including 45 (24.2%) biological process (BP) terms (of
   which, 88 were significant; corrected P-values < 0.05), 36 (19.4%)
   cellular component (CC) terms, and 105 (56.5%) molecular function (MF)
   terms (Fig. [53]2, Supplementary Table [54]S2). Most DEGs in the BP
   category were putatively attributed to “signal transduction”, “single
   organism signaling”, “cell communication”, and “signaling”, those in
   the CC category were attributed to “extracellular region” and
   “cytoskeleton” and those in the MF category were attributed to
   “structural molecule activity”, “metal ion binding”, “cation binding”,
   and “structural constituent of cuticle”.

Figure 2.

   [55]Figure 2
   [56]Open in a new tab

   Gene ontology classification of differentially expressed genes.
   Category of biological process (BP); Category of cellular component
   (CC); Category of molecular function (MF).

   Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation
   revealed that 100 DEGs were annotated for M. persicae. These annotated
   genes were classified into 45 groups based on the secondary pathway
   hierarchy (Supplementary Table [57]S3). The major pathways included
   metabolic pathways (21 proteins), oxidative phosphorylation (9
   proteins), Hippo signaling pathway – fly (5 proteins), and
   phenylalanine metabolism (4 proteins; Fig. [58]3). In the database of
   the present study, phenylalanine metabolism (api00360), oxidative
   phosphorylation (api00190), Hippo signaling pathway – fly (api04391),
   tyrosine metabolism (api00350), and ECM-receptor interaction (api04512)
   pathways were found to be with higher corrected P-values < 0.05. These
   annotations establish an invaluable basis for elucidating the specific
   processes, functions, and pathways involved in the imidacloprid
   resistance of M. persicae.

Figure 3.

   [59]Figure 3
   [60]Open in a new tab

   Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment
   analysis.

qRT- PCR validation of DEGs

   qRT-PCR analysis of 10 randomly selected genes confirmed the expression
   trends observed in the RNA-Seq results (Fig. [61]4, Supplementary
   Table [62]S4), thereby, suggesting that the DEG profiles were reliable.
   Increases in the expression of the selected genes were small, and the
   biological relevance of these changes is likely minimal.

Figure 4.

   [63]Figure 4
   [64]Open in a new tab

   Quantitative real-time PCR (qRT-PCR) validation of the expression of
   differentially expressed genes identified using RNA-sequencing. The
   expression levels were normalised to GAPDH, EF1α, and β-actin genes.

Discussion

   The green peach aphid, M. persicae, is an economically important pest
   that is typically controlled using insecticides. However, the
   irrational use of imidacloprid has promoted the rapid development of
   insecticide resistance in many M. persicae populations, resulting in
   the failure in controlling the pest^[65]6–[66]8. Undoubtedly, the
   development of insecticide resistance is complex and is a common
   adaptation to insecticide exposure. This study aimed to identify
   imidacloprid-responsive genes in M. persicae, to establish a basis for
   investigating the responses and adaptive physiological changes that
   contribute to insect resistance.

   For insect pests, insecticide exposure can be considered a common
   environmental stress and the cuticle is the first and main barrier for
   insects against environmental stresses. As such, cuticular proteins
   have been reported to play crucial roles in the insecticide resistance
   and tolerance of a variety of insect species, including M. persicae,
   Lymantria dispar, Aphis gossypii and Plutella
   xylostella^[67]13,[68]16–[69]18. Indeed, cuticle protein genes (larval
   cuticle protein A2B/A3A, 111039230/111039229; cuticle protein
   111027628, 111034584, 111038352; cuticle protein 12.5/65/38/21/7/3,
   111042794/111035272/111032404/111039228, 111039055/111039056/111032178;
   endocuticle structural glycoprotein SgAbd-2/4/8/9, 111031116/111028315,
   111031121/111031114/111041497) were down-regulated in the IR strain of
   M. persicae, when compared to the expression levels in the IS strain.
   In this study, the down-regulation of cuticle protein genes in the IR
   strain was consistent with previous studies that reported the
   down-regulation of eight cuticle protein genes in
   thiamethoxam-resistant A. gossypii^[70]17. Such changes in expression
   suggest that cuticle protein genes contribute to the protection of M.
   persicae from the mechanical damage caused by imidacloprid exposure.

   The ABC transporters are responsible for the translocation of a variety
   of substrates (e.g. nutrients, lipids, inorganic ions and xenobiotics)
   and can be categorised into eight subfamilies from ABCA to ABCH^[71]19.
   The ABC transporters present in the blood-brain barrier of insects can
   protect the nervous system from insecticides^[72]20 and are reportedly
   involved in insecticide resistance^[73]21,[74]22. Indeed, the
   contribution of ABC transporters to insecticide resistance has been
   reported for 27 different insecticides including, imidacloprid,
   pyrethroids, and avermectin^[75]21,[76]23,[77]24. With regard to
   imidacloprid, ABC transporters have been reported to enhance its
   exclusion from the brain and to hinder its access to target
   sites^[78]24. Recently, the down-regulation of these genes was reported
   to be linked to insecticide resistance^[79]22,[80]25. Similarly, in our
   study, it was found that ABCG23 (ABCG23, 111040927) was down-regulated
   in IR M. persicae, which suggests that the ABC transporter is involved
   in imidacloprid metabolism and transport.

   Glutathione S-transferases (GSTs) are a widespread superfamily of genes
   that occur in almost all living organisms and participate in a variety
   of cellular physiological processes, including the detoxification of
   harmful endobiotic and xenobiotic compounds. Insect GSTs are generally
   categorised into two main groups: cytosolic and microsomal, based on
   their cellular location. In insects, most GSTs are cytosolic proteins
   and are classified as delta, epsilon, omega, sigma, theta, and
   zeta^[81]26. In recent years, a number of studies have investigated the
   correlation between insect GST genes and insecticide resistance. These
   studies have demonstrated that insect GSTs play important roles in
   insecticide detoxification and eliminate the oxidative stress caused by
   insecticide exposure^[82]27. The up-regulation of GST genes has been
   associated with insecticide detoxification^[83]28. However, the
   down-regulation of GST genes following insecticide exposure has also
   been reported in several insect species. In Leptinotarsa decemlineata,
   for example, the expression of LdGSTe4 and LdGSTe6 was significantly
   down-regulated after cyhalothrin, fipronil, and endosulfan exposure,
   and that of PxGSTd2, PxGSTe2, PxGSTe5, PxGSTo1, PxGSTs1, PxGSTs2, and
   PxGSTt1 was down-regulated by β-cypermethrin exposure^[84]29.
   Similarly, our results indicated that GST (111036096, 111036826) was
   down-regulated in IR M. persicae. Bautista et al. suggested that this
   phenomenon may be an adaptive mechanism to insecticide pressure and an
   energy trade-off strategy to ensure that energy is allocated to the
   most effective genes responsible for detoxification when stimulated by
   insecticide exposure^[85]30. This suggests that GST plays a relatively
   minor role in the imidacloprid metabolism of M. persicae. However,
   because insect GST genes exhibit a variety of expression responses to
   insecticide exposure^[86]29, functional studies are needed to examine
   this hypothesis and identify the specific GST members involved in
   insecticide detoxification.

   Trypsin is a serine protease responsible for digestion; it also
   contributes to insecticide detoxification^[87]31. Recently, Zhu et al.
   (2015) reported that the midgut trypsin activity of Bt-resistant
   Spodoptera frugiperda was relatively lower than that of a susceptible
   strain^[88]32. Indeed, in our study, trypsin gene (111033016) was also
   down-regulated in IR M. persicae. These findings suggest that the
   expression and function of trypsin are associated with insecticide
   resistance. However, little is known about the exact role of trypsin in
   imidacloprid resistance. It has been demonstrated that the lack of
   midgut trypsin made some insects adapt to insecticide toxins by a
   mechanism where incomplete or non-activation of the protoxin occurs,
   and finally induced resistance development^[89]33.

   Mitochondria play critical roles in a variety of cellular processes,
   among which energy generation is the most critical and which is
   primarily achieved through coupled oxidative phosphorylation. The
   electron transport chain (ETC) consists of four macromolecular protein
   complexes (complex I–IV), that coordinate to maintain mitochondrial
   inner membrane potential^[90]34. In this study, the mitochondrial
   complex-related genes NADH dehydrogenase (ETC I, 111029398), succinate
   dehydrogenase (ETC II, 111027244), cytochrome b-c1 complex (ETC III,
   111038823, 111027246), and cytochrome c oxidase (ETC IV, 111034852)
   were all down-regulated in IR M. persicae, which indicated that
   imidacloprid exposure reduced the expression of ETC I, II, III, and IV
   component genes. It is possible that the respiration and energy
   production of IR aphids may have been weakened by imidacloprid
   exposure. There is no evidence that the ETC contributed to the enhanced
   imidacloprid tolerance, but the apparent alteration in the expression
   of complex I–IV in IR M. persicae clearly associates the mitochondrial
   complex-related genes with imidacloprid resistance. Previous studies
   have suggested that the overexpression of cytochrome P450s is the
   primary reason for neonicotinoid resistance^[91]35,[92]36. Indeed,
   CYP6CY3, which is a cytochrome P450 gene in M. persicae, has been
   suggested to play a primary role in the development of insecticide
   resistance^[93]4,[94]13. However, the association between mitochondrial
   complexes and cytochrome P450s with the resistance of M. persicae to
   imidacloprid requires further investigation.

   In conclusion, this study provides, to our knowledge, the first
   description of genes related to imidacloprid resistance in M. persicae
   using an RNA-Seq approach. The results indicate that the response
   patterns of aphids are complex during the development of imidacloprid
   resistance, as demonstrated by changes in the expression of genes
   involved cuticle structure, binding, metabolic processes, and oxidative
   phosphorylation. Further investigations are needed to assess the
   specific roles of these genes in the response of M. persicae to the
   stress of insecticide exposure. These findings will provide a basis for
   investigating the development and mechanisms of insecticide resistance.

Materials and Methods

Aphid strains

   One IS strain and one IR strain of M. persicae were used. The IS strain
   was obtained from tobacco in Guizhou province, China, in 2009, and was
   subsequently reared in the absence of insecticides. The IR strain was
   generated from the IS population under continuous imidacloprid
   selection pressure^[95]37. In this study, the imidacloprid resistance
   of the IR strain was ~45-fold greater than that of the IS strain. Both
   the IS and IR strains were maintained on tobacco plants at 23–25 °C,
   with a 16-h photoperiod (16 h light, 8 h dark) and relative humidity of
   60%.

RNA isolation

   Total RNA was isolated using TRIzol reagent (Invitrogen, CA, USA). The
   RNA concentration and purity were measured using a NanoPhotometer^®
   spectrophotometer (IMPLEN, CA, USA). The integrity of RNA was confirmed
   using the Bioanalyzer 2100 system (Agilent, CA, USA).

Library preparation and sequencing

   DGE-Seq was performed using the mRNA isolated from the IS and IR
   strains, with three biological replicates per strain. Sequencing
   libraries were constructed using NEBNext^® Ultra^TM RNA Library Prep
   Kit for Illumina^® (NEB, USA). Briefly, Poly-T oligo-attached magnetic
   beads were used to purify mRNA, which was then fragmented using
   divalent cations, and first-strand cDNA was synthesised using random
   primers and reverse transcriptase, whereas second-strand was
   synthesised using DNA Polymerase I and RNase H. After 3′-end
   adenylation, the cDNA fragments were ligated to adapters and then
   selectively enriched by PCR. The libraries were purified using the
   AMPure XP system. The quality of the sample libraries was assessed
   using the Agilent Bioanalyzer 2100 system. Finally, DGE sequencing was
   implemented using an Illumina HiSeq 2000 instrument.

Read mapping and expression quantification

   Raw reads in fastq format were processed using in-house Perl scripts.
   The adapter, ploy-N, and low-quality sequences from raw reads were
   eliminated to obtain clean reads. Q score, GC content, and sequence
   duplication level were calculated to obtain high-quality clean reads,
   which were subsequently used for all the downstream analyses. The clean
   reads were mapped to the M. persicae genome (GenBank no.
   GCA_001856785.1), and the expression levels of the genes were
   calculated using fragments per kilobase per million reads (FPKM)
   values^[96]38.

Differentially expressed gene (DEG) analysis and annotation

   Differential expression analysis was performed using DESeq^[97]39, and
   the Benjamini and Hochberg’s method was used to adjust the resulting
   P-values, to minimise the false detection rate (FDR). Genes were
   identified as differentially expressed (i.e. DEGs) if the adjusted
   P-value was <0.05^[98]40. The functional annotation and classification
   of the genes were performed using the GO database, and biological
   pathway annotations were obtained using the KEGG database.

Quantitative real-time PCR analysis

   Quantitative real-time PCR (qRT-PCR) analysis was performed to validate
   the expression profiles of 10 randomly selected DEGs and three internal
   controls, namely the elongation factor 1α (EF1α), β-actin, and
   glyceraldehyde-3-phosphate dehydrogenase (GAPDH) genes. The primers
   used in qRT-PCR are summarized in Table [99]2. qRT-PCR reactions were
   performed using a SYBR Premix DimerEraser Kit (Takara, Dalian, China)
   on an ABI 7500 system (ABI, CA, USA). All qRT-PCR experiments were
   performed in triplicate using independent samples. The expression
   levels were determined by the 2^−ΔΔCt method^[100]41, using the
   geometric mean of three selected internal control genes for
   normalisation.

Table 2.

   Primers used for qRT-PCR validation of differentially expressed genes.
   Gene name Forward primer (5′-3′) Reverse primer (5′-3′) Gene
   description
   LOC111034675 5’TGCGGGAGGTGTGAGAGCTG TCGCCGTTTTTCAATGTATCG maltase
   A3-like
   LOC111036239 CGCGGTACATGAATTGCACAACTG ACGCAATGTCGAAGAACGGTATC aspartate
   aminotransferase, cytoplasmic-like
   LOC111032872 CCGCGTGAGGATATGTGTTGAC ACGGCCAGAGGACACACGATG protein
   O-linked-mannose beta-1, 2-N-acetylglucosaminyltransferase 1-like
   LOC111040357 GAGCCAAGAAAATGCAGATGAATAC TCCGCATGAATGAGACCCAAATC
   homeodomain-interacting protein kinase 2
   LOC111028270 TCCCGGGTTTATCGTGGCAAG CCCAACAACATGAGCAACAAATAAC
   serine/threonine-protein kinase tousled-like 2
   LOC111029398 CGCCCGATGCCATTAGTTCAAC TGGCATTCAAGTCATCTGTCTCATC NADH
   dehydrogenase [ubiquinone] 1 beta subcomplex subunit 11,mitochondrial
   LOC111036096 GGCAGCATACAAACTCACTTACTTC AGGGCATTTTGGGCTTGATTG
   glutathione S-transferase-like
   LOC111042314 GCCGCCGAACAGTCTGCAAAAC AGCGGCCAGGTAGGGTGAAG
   serine/arginine-rich splicing factor 8-like
   LOC111037432 CGCCAAAACATCAACAATCAACAAG GGTTGGCGTGGTTGTTAAGATTTG GATA
   zinc finger domain-containing protein 14-like
   LOC111039230 GCGACGACGTGACCGGTTACTAC TGGCGGCCTTAGCGACGATC larval
   cuticle protein A2B-like
   EF1α CCGATGTCTATGTCTGCTAAGG CATGATTTGAGCCTCGCCAA
   β-actin CGGTTCAAAAACCCAAACCAG TGGTGATGATTCCCGTGTTC
   GAPDH GCGGTTTCGACGTGTCAGTTTG CCGGAGCCCACAATGCACAC
   [101]Open in a new tab

Supplementary information

   [102]Supplementary Table S1^ (79KB, xls)
   [103]Supplementary Table S2^ (92KB, xls)
   [104]Supplementary Table S3^ (28.5KB, xls)
   [105]Supplementary Table S4^ (29.5KB, xls)

Acknowledgements