Identification and Functional Analysis of Differentially Expressed Genes in Myzus persicae (Hemiptera: Aphididae) in Response to Trans-anethole.

Plant essential oils, with high bioactivity and biodegradability, provide promising alternatives to synthetic pesticides for pest control. Trans-anethole is the major component of essential oil from star anise, Illicium verum Hook. The compound has a strong contact toxicity against the green peach aphid, Myzus persicae (Sulzer) (Hemiptera: Aphididae), which is a major insect pest of many vegetables and crops. However, little information is known about how M. persicae responds to trans-anethole at the molecular level. We conducted a comparative transcriptome analysis of M. persicae in response to a LD50 dose of trans-anethole. A total of 559 differentially expressed genes were detected in the treated individuals, with 318 genes up-regulated, and 241 genes down-regulated. Gene ontology (GO) analysis revealed that these genes were classified into different biological processes and pathways. We also found that genes encoding ATP-binding cassette (ABC) transporters, DnaJ, and cuticle proteins were dramatically up-regulated in response to trans-anethole. To study the function of these genes, we performed RNA interference (RNAi) analysis. Knockdown of an ABC transporter gene (ABCG4) and a DnaJ gene (DnaJC1) resulted in a significantly increased mortality rate in M. persicae following trans-anethole exposure, indicating the involvement of these two genes in the toxicity response to trans-anethole. The findings provide new insights into the mechanisms of M. persicae in coping with plant essential oils.

The green peach aphid, Myzus persicae (Sulzer) (Hemiptera: Aphididae), is a serious agricultural pest attacking many vegetables and crops (van Emden et al. 1969). M. persicae causes damage by direct sap feeding and by transmission of plant pathogenic viruses (Lai et al. 2017). Synthetic insecticides are commonly used to control M. persicae, but repeated use of pesticides has often resulted in insecticide resistance (Tang et al. 2017). Many studies have documented the resistance of M. persicae populations to organophosphate, pyrethroid, and neonicotinoid insecticides (Bass et al. 2014). Extensive pesticide use has also increased insecticide pollution and residues in the environment (Guedes et al. 2016). Environmentally-friendly insecticides are needed to manage the population of M. persicae.

Essential oils are plant secondary metabolites that are important for defense against herbivores and microbial pathogens (Isman 2020). Due to their relatively high bioactivity and biodegradability, plant essential oils have been widely used as natural pesticides to control agricultural, stored-product, and urban pests (Liao et al. 2016, Gao et al. 2019, Wang et al. 2020). Trans-anethole [trans-1-methoxy-4-(1-propenyl) benzene] is a phenylpropanoid that is an abundant component of the essential oil of star anise, Illicium verum Hook (Huang et al. 2010). Trans-anethole probably evolved to protect plants against herbivores and pathogens (Deng and Lu 2017). This compound has high bioactivity against many pest insects, such as the oriental fruit fly Bactrocera dorsalis (Hendel) (Diptera: Tephritidae), fall armyworm Spodoptera frugiperda (J.E. Smith) (Lepidoptera: Noctuidae), red flour beetle Tribolium castaneum (Herbst) (Coleoptera: Tenebrionidae), cowpea weevil Callosobruchus maculatus (Fabricius) (Coleoptera: Bruchidae), and the yellow fever mosquito Aedes aegypti L. (Diptera: Culicidae) (Chang et al. 2009, Cruz et al. 2017, Hashem et al. 2018; Pandiyan et al. 2019, Barbosa et al. 2020). Exposure to trans-anethole inhibits acetylcholinesterase (AChE) activity in the flour moth Ephestia kuehniella, and treatment with trans-anethole causes histological damage in T. castaneum (Hashem et al. 2018, Shahriari et al. 2018). Thus, trans-anethole has potential as a green insecticide active ingredient for pest control.

Exposure to plant essential oils alters gene transcription levels in insects, and differentially expressed genes (DEGs) might be involved in insect defense mechanisms against essential oils or their components (Huang et al. 2018a, Gao et al. 2020). Comparative transcriptome analysis, based on high-throughput RNA sequencing technology, can be used to determine transcriptional variation of numerous genes simultaneously (Wang et al. 2009). This method has been used to compare gene expression profiles in insects responding to essential oils, synthetic pesticides, and environmental stressors (Silva et al. 2012, Huang et al. 2018a, Xiong et al. 2019, Yang et al. 2021).

RNA interference (RNAi) is a method of post-transcriptional gene silencing (Fire et al. 1998). In insects, RNAi can be used to elucidate the physiological function of genes (Bellés 2010). For example, in T. castaneum, knockdown of two chitin synthase genes by RNAi revealed their specialized roles in synthesis of the cuticle and the midgut peritrophic matrix (Arakane et al. 2005). Silencing a gene encoding NADPH-cytochrome P450 reductase increased the susceptibility of Nilaparvata lugens (Stål) (Hemiptera: Delphacidae) to insecticides (Liu et al. 2015). In M. persicae, RNAi has been used to study the function of a number of genes including hunchback, erythropoietin-producing hepatocellular receptor, voltage-gated sodium channel, and cuticle protein genes (Mao and Zeng 2014, Mulot et al. 2016, Tariq et al. 2019, Shang et al. 2020).

We previously demonstrated that trans-anethole has strong contact toxicity against apterous adults of M. persicae. Thus, trans-anethole could be utilized to develop effective aphicide formulations (Li et al. 2017). However, the effect of trans-anethole on the gene expression of M. persicae is unknown and the function of aphid genes in defense against trans-anethole is not clear. In this study, we analyzed the transcriptome of M. persicae in response to trans-anethole exposure, and identified candidate genes involved in the detoxification and defense response using RNAi. Our findings provide insights into the possible mechanisms used by M. persicae to respond to trans-anethole exposure.

Materials and Methods

Aphid Rearing

The M. persicae individuals used in this study originated from a colony collected from a cabbage field in an experimental farm at Anhui Agricultural University, Hefei, China in 2016 (Li et al. 2017). The aphids were maintained in an insectary without contact with insecticides. Aphids were fed potted cabbage (Brassica oleracea var. capitata). Cabbages (30-day-old) were used for rearing aphids, and new plants were provided weekly. The rearing conditions were 25 ± 1°C, 60 ± 5% RH, and a photoperiod of 16:8 (L:D) h photoperiod.

Insect Treatment

Trans-anethole (CAS number: 4180-23-8; purity: 99%; purchased from Aladdin Co., Ltd., Shanghai, China) was diluted with methanol to make a stock solution. The stock solution was serially diluted with distilled water to create a working solution at a final concentration of 2.09 μmol/L (0.31 mg/L). This concentration is the 72 h LD50 dose determined previously (Li et al. 2017). Aphid treatment with trans-anethole was performed using the leaf-dip method (Kang et al. 2016). Briefly, leaf discs (3-cm diameter) were made from fresh cabbage leaves (30-day-old plants), and 10 aphids (one-day-old apterous adults) were transferred onto each leaf disc with a fine brush. Aphids and the leaf discs were dipped into the trans-anethole solution for 10 s, and residual droplets were absorbed with clean filter paper. Aphids in the control group were treated with distilled water. After exposure, aphids were reared under normal rearing conditions [25 ± 1°C, 60 ± 5% RH, and a photoperiod of 16:8 (L:D) h]. Surviving aphids at 72 h after treatment were frozen in liquid nitrogen and stored at −80°C until RNA extraction was conducted. The experiment was biologically repeated three times, and each replicate included at least 100 individuals.

RNA Preparation and Transcriptome Sequencing

Total RNA was extracted from the whole body of M. persicae using RNAiso Plus reagent (Takara, Dalian, China) following manufacturer’s protocol. The total RNA was then treated with RNase-free DNase I (Takara) to remove genomic DNA contamination. The quality and concentration of each total RNA sample were determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). The Hieff NGS MaxUp Dual-mode mRNA Library Prep Kit for Illumina (Yeasen, Shanghai, China) was used for cDNA library construction. The mRNA was purified from 10 μg total RNA using oligo (dT) magnetic beads and fragmented into short sequences. Then, first-strand cDNA was synthesized using a random hexamer primer, followed by second-strand cDNA synthesis. After end repairing and adaptor ligation, the second-strand cDNA was amplified by PCR and concentrated by Hieff NGS DNA selection beads (Yeasen, Shanghai, China) to create a cDNA library. The prepared cDNA libraries were sequenced on the Illumina HiSeq2000 platform to generate 2 × 100 bp paired-end reads. Each library was sequenced in a single-flow cell lane.

Data Filtering and Reads Mapping

To obtain the trimmed reads, raw data generated from the Illumina platform were filtered by removing adaptors, low-quality reads (reads with >20% of the bases with quality scores <20), and reads with unknown nucleotides (N bases > 5%) using Trimmomatic v0.35 software (https://github.com/usadellab/Trimmomatic) (Bolger et al. 2014). The reference genome (assembly name: M. persicae O genome v2.0; assembly date: May 5, 2020) and gene model annotation file (name: M. persicae O OGS2.0) of M. persicae clone O were downloaded from AphidBase (https://bipaa.genouest.org/sp/myzus_persicae/), and HISAT2 v2.0.5 software (https://github.com/DaehwanKimLab/hisat2) was used as the mapping tool to build the index of the reference genome and align the trimmed reads with the genome sequences (Kim et al. 2019). The default parameters of HISAT2, with a maximum of two mismatches were used. In addition, gene splicing events were taken into account by HISAT2 and the maximum allowed intron length during read mapping was set to 1,000,000.

Analysis of DEGs

DEGs were analyzed using the DEGseq v3.0 software (http://bioconductor.org/packages/release/bioc/html/DEGseq.html) (Wang et al. 2010). The gene expression within each group was normalized by the transcript per million (TPM) method using the formula: Normalized expression = mapped readcount/Total reads × 1,000,000. The significantly up- and down-regulated genes were selected according to the threshold: log2 fold-change (treatment/control) > 2.0 and P-value < 0.05. The P-value was calculated using Fisher’s exact test and likelihood ratio test (Wang et al. 2010). Pearson’s correlation analysis and principal component analysis (PCA) were performed by Corrplot and Vegan packages, respectively, in R software (version 3.1.3) using the log-transformed TPM values of all DEGs. To better understand the function of DEGs in M. persicae, gene ontology (GO) enrichment analysis was performed to investigate the biological processes and pathways that the DEGs were involved in. ClusterProfiler v3.0.2 software (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) was used to analyze GO terms, with default parameters and corrected P-value < 0.05 (Yu et al. 2012). The ClusterProfiler package provides groupGO, a gene classification method to classify genes based on their projection at a specific level of the GO corpus, and enrichGO, a tool to calculate an enrichment test for GO terms based on hypergeometric distribution. The compareCluster tool in the ClusterProfiler package was used for automatically calculating enriched functional categories of each gene cluster and for visualization (Yu et al. 2012).

Reverse Transcription Quantitative PCR (RT-qPCR) Validation

Another cohort of M. persicae (different from those used in transcriptome sequencing) treated with 2.09 μmol/L trans-anethole solution and distilled water (control; see ‘Insect treatment’ section) was used for RT-qPCR analysis. Total RNA was extracted as described above, and first-strand cDNA was synthesized from 1 μg of total RNA using a ReverTra Ace qPCR RT kit (Toyobo, Osaka, Japan) following manufacturer’s instructions. Primers for RT-qPCR are listed in Supp Table S1 (online only). The housekeeping genes β-actin and 18S rRNA were used as internal references to normalize target gene expression (Kang et al. 2017). In a pretest, the PCR products amplified by these primers were sequenced to verify the amplification of the correct targets. In addition, the amplification efficiency (Supp Table S1 [online only]) of each pair of primers was calculated from the slope of the standard curve generated from a serial dilution of cDNA samples using the formula: E = 10[–1/slope] – 1 (Bustin et al. 2009).

RT-qPCR was performed using SYBR Green Real-Time PCR Master Mix (Toyobo). Each reaction (20 μL volume) contained 10 μL SYBR Green Mix (Toyobo, Osaka, Japan), 0.4 μL (0.2 μM) sense primer, 0.4 μL (0.2 μM) antisense primer, 1 μL (10 ng) cDNA template, and 8.2 μL nuclease-free water. Reactions were run on a CFX96 Real-Time System (Bio-Rad, Hercules, CA) with the following parameters: one cycle at 95°C for 30 s, and 40 cycles at 95°C for 5 s and 60°C for 25 s. A heat-dissociation protocol was added at the end of the thermal cycle to confirm that only a single gene was detected by the fluorescent dye. A no-template control and a no-reverse transcriptase control were also included to detect potential contamination (Bustin et al. 2009). The RT-qPCR experiments were biologically repeated three times, and relative expression levels for target genes were normalized to the geometric mean of two reference genes (Vandesompele et al. 2002). The formula used is as follows:

where ‘E’ refers to the primer efficiency.

DsRNA Synthesis and RNAi

Gene-specific primers containing T7 RNA polymerase promoter sequences (Supp Table S1 [online only]) were designed and used to amplify a short fragment (~400 bp) from each of the selected genes. A 427-bp fragment of green fluorescent protein (GFP; control) gene was amplified from the pEGFP-N1 vector (Clontech, Mountain View, CA; Supp Table S1 [online only]). The PCR products were run on an agarose gel, and DNA bands of the expected size were purified and used as templates to synthesize the dsRNA. The synthesis was carried out using the MEGAscript T7 Transcription Kit (Thermo Fisher Scientific, Wilmington, DE) following manufacturer’s instructions. The quality and concentration of the dsRNAs were verified by agarose gel electrophoresis and a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific). The dsRNAs were diluted to a final concentration of 4 μg/μL with nuclease-free water and stored at –80°C. One-day-old apterous aphids were anesthetized on ice, and 13.8 nL (55.2 ng) of the dsRNA was injected into the intersegmental membrane between the prothorax and the mesothorax, using a Nanoinject II automatic injector (Drummond Scientific, Broomall, PA) with glass capillary needles. Aphids injected with GFP-dsRNA (dsGFP) were used as controls. Afterwards, aphids were transferred to fresh cabbage leaves and reared at 25 ± 1°C, 60 ± 5% RH, and a photoperiod of 16:8 (L:D) h. Variations in gene expression were determined by RT-qPCR at 24 h, 48 h, and 72 h following the dsRNA injection. Individuals at 48 h following the injection were exposed to a LD50 dose of trans-anethole (as described in ‘Insect treatment’ section), and mortality was recorded 72 h post trans-anethole exposure. Aphids were considered dead if they did not move when probed with a fine brush. The experiments were repeated biologically three times, and 50 individuals were used in each replicate.

Statistics

The RT-qPCR results and mortality data were analyzed using Data Processing System (DPS) software v9.5 (Tang and Zhang 2013). Student’s t-test was used to analyze the difference between two samples. The level of significance was set at P < 0.05.

Results

Transcriptome Sequencing and Assembly

The cDNA libraries constructed from trans-anethole-treated and control individuals were sequenced on an Illumina platform, which generated 41.4–48.7 million raw reads for treated aphids and 47.5–51.6 million raw reads for untreated aphids (Table 1). After data filtration, 38.6–45.7 million and 45.3–48.5 million trimmed reads were produced from the treated individuals and controls, respectively (trimmed reads ratio > 77%; Table 1). About 93% of the trimmed reads of trans-anethole-challenged and control individuals were mapped directly to the reference genome data of M. persicae (Table 1). The data presented here have been deposited into NCBI’s Sequence Read Archive (SRA) under accession numbers SRX10441482–SRX10441487.

Table 1.

Summary of statistical data for the transcriptomes of M. persicae

	Treatment 1	Treatment 2	Treatment 3	Control 1	Control 2	Control 3
Number of raw reads	41.4 × 106	48.7 × 106	46.9 × 106	49.3 × 106	47.5 × 106	51.6 × 106
Length of raw reads (Mb)	6.20	7.31	7.03	7.40	7.12	7.74
Number of trimmed reads	38.6 × 106	45.7 × 106	44.0 × 106	45.9 × 106	45.3 × 106	48.5 × 106
Length of trimmed reads (Mb)	5.18	6.10	5.80	6.28	5.48	6.48
Trimmed reads ratio (%)	83.5	83.5	82.5	84.8	77.0	83.7
Q20 (%)	97.3	97.5	95.7	97.5	97.8	97.5
GC (%)	48.1	46.7	44.6	45.3	50.9	42.3
Total map	36.3 × 106 (93.8%)	42.4 × 106 (92.7%)	41.3 × 106 (93.7%)	42.8 × 106 (93.2%)	42.1 × 106 (92.8%)	45.3 × 106 (93.2%)

Treatment 1–3 and control 1–3 represent different biological repeats, respectively. Raw reads, the original sequence data; trimmed reads, the filtered sequencing data; trimmed reads ratio, the proportion of trimmed reads to the total raw reads; Q20 (%), the percentage of bases with a Phred value of >20; GC (%), the percentage of the number of guanine and cytosine in the total bases; total map, the number of trimmed reads that can be located on the genome.

In total, 559 DEGs were detected in the trans-anethole-treated aphids, compared with the numbers in the control groups (Supp Fig. S1 [online only]; Supp Table S2 [online only]). Pearson correlation analysis indicated that there was a positive correlation between the biological replicates (Supp Table S3 [online only]). The PCA results showed that 81.9% of the transcriptional changes were related to the trans-anethole treatment (Supp Fig. S2 [online only]). Of these, the first (PC1), second (PC2), and third (PC3) principal components explained 42.76%, 26.93%, and 12.23% of the variance, respectively (Supp Fig. S2 [online only]). Among the DEGs, 318 genes were up-regulated and 241 genes were down-regulated (Supp Fig. S1 [online only]; Supp Table S2 [online only]). The most up-regulated gene was acyl-CoA synthetase, followed by serine/threonine-protein kinase (Table 2). Most of the top 10 down-regulated genes were related to rRNA-processing and enzymatic catalysis (Table 2). In addition, we found that genes annotated as ATP-binding cassette (ABC) transporters (ABCG1 and ABCG4), DnaJ (DnaJC1), and cuticle proteins (CP19, CP7, CP65, and CP12.5) were dramatically up-regulated in response to trans-anethole (Table 3). The full repertoire of DEGs is listed in Supp Table S2 (online only).

Table 2.

Top 10 up- and down-regulated genes in M. persicae responding to trans-anethole

Gene ID	Log2 fold-change	P-value	Result	Gene description
XM_022320847	15.94	1.56E-17	Up	Acyl-CoA synthetase family member 3, mitochondrial, transcript variant X1
XM_022307950	15.66	1.87E-04	Up	Serine/threonine-protein kinase STE20-like, transcript variant X3
XM_022318480	15.60	4.22E-08	Up	Flightin, transcript variant X1
XM_022308595	14.98	2.09E-04	Up	Sodium-independent sulfate anion transporter-like, transcript variant X1
XM_022326720	14.82	8.89E-09	Up	Double-stranded RNA-binding protein Staufen homolog 2, transcript variant X1
XM_022320225	14.75	3.51E-04	Up	Tetratricopeptide repeat protein 27, transcript variant X2
XM_022308145	14.74	7.45E-07	Up	Gamma-aminobutyric acid type B receptor subunit 2, transcript variant X1
XM_022311808	14.72	2.85E-04	Up	Tudor domain-containing protein 7-like, transcript variant X2
XM_022312345	14.72	2.67E-03	Up	Uncharacterized LOC111032120, transcript variant X1
XM_022311917	14.59	1.29E-05	Up	Signal recognition particle subunit SRP68, transcript variant X2
XM_022320783	−14.68	2.57E-08	Down	rRNA-processing protein UTP23 homolog
XM_022310374	−14.49	1.95E-02	Down	Gamma-soluble NSF attachment protein
XM_022320935	−14.46	5.49E-05	Down	Chromodomain-helicase-DNA-binding protein Mi-2 homolog
XM_022305482	−14.34	3.00E-02	Down	CDGSH iron-sulfur domain-containing protein 2 homolog A
XM_022321141	−14.23	3.10E-03	Down	Peptidylglycine alpha-hydroxylating monooxygenase,
XM_022305972	−14.20	4.21E-03	Down	Esterase FE4-like
XR_002604214	−14.20	2.23E-02	Down	Uncharacterized LOC111031524
XM_022321589	−14.10	4.92E-03	Down	Adenylyltransferase and sulfurtransferase MOCS3
XM_022304755	−14.05	4.80E-07	Down	AP-3 complex subunit sigma-1
XM_022316265	−14.01	1.26E-02	Down	Uncharacterized LOC111034863

Table 3.

DEGs encoding ABC transporters, DnaJ, and cuticle proteins from M. persicae responding to trans-anethole

Gene ID	Log2 fold-change	P-value	Gene description
XM_022310265	14.06	1.56E-03	ATP-binding cassette sub-family G member 1 (ABCG1)
XM_022323461	12.37	2.54E-02	ATP-binding cassette sub-family G member 4 (ABCG4)
XM_022306112	7.30	2.85E-04	DnaJ homolog subfamily C member 1 (DnaJC1)
XM_022315866	5.68	3.19E-03	Cuticle protein 19(CP19)
XM_022322620	5.42	1.08E-06	Cuticle protein 7(CP7)
XM_022316814	4.70	1.77E-06	Cuticle protein 65(CP65)
XM_022326838	2.69	6.63E-03	Cuticle protein 12.5(CP12.5)

GO Analysis

GO enrichment analysis was performed to elucidate the function of the DEGs. As shown in Supp Fig. S3 (online only), the 559 DEGs were classified in three categories, namely biological process, cellular component, and molecular function. Within the biological process category, the DEGs were mostly enriched in the term ‘cellular process’, followed with the terms ‘biological regulation’ and ‘metabolic process’. While of the cellular component category, the terms ‘cell’ and ‘cell part’ were the most abundant. For the molecular function category, most DEGs were related to the term ‘binding’ (Supp Fig. S3 [online only]).

RT-qPCR Validation of DEGs

To validate the DEG data, seven genes were selected and subjected to RT-qPCR experiments. These genes included two ABC transporter genes (ABCG1 and ABCG4), one DnaJ gene (DnaJC1), and four cuticle protein genes (CP19, CP7, CP65, and CP12.5). All of these genes displayed up-regulation after trans-anethole treatment based on DEG data (Table 3). As expected, the RT-qPCR results were consistent with the transcriptomic data (Fig. 1).

Fig. 1.

Verification of DEGs in trans-anethole-treated and control individuals using RT-qPCR. Black-filled columns represent the fold-change of gene expression based on comparative transcriptome analysis, whereas grey-filled columns represent data from RT-qPCR analysis. Error bars represent standard error (SE). ‘*’ denotes a significant difference in expression levels between treated and control insects (Student’s t-test, P < 0.05).

Effect of Gene Knockdown on Trans-anethole Susceptibility

We conducted RNAi analysis to study the function of seven DEGs (ABCG1, ABCG4, DnaJC1, CP19, CP7, CP65, and CP12.5) in response to trans-anethole. To verify the specificity of RNAi, dsRNA sequences corresponding to ABC transporter and cuticle protein genes were aligned, and the results showed no consecutive identical nucleotides between pairs of the dsRNA fragments (Supp Fig. S4 [online only]). After injection of dsRNAs, the expressions of all tested genes were successfully suppressed from 24 h to 72 h compared with the mRNA levels in the control group (Fig. 2). Of these, the transcription levels of ABCG1, DnaJC1, CP7, and CP12.5 decreased rapidly at 24 h, whereas RNAi of ABCG4, CP19, and CP65 decreased slowly at 48 h (Fig. 2).

Fig. 2.

Relative expression levels of selected genes (ABCG1, ABCG4, DnaJC1, CP19, CP7, CP65, and CP12.5) in M. persicae after injection of corresponding dsRNAs. Individuals injected with the green fluorescent protein (GFP) dsRNA were used as controls. Data are presented as mean ± SE. ‘*’ represents a significant difference between two samples (Student’s t-test, P < 0.05); ‘ns’ represents ‘not significant’ (Student’s t-test, P > 0.05).

Among the tested genes, silencing of ABCG4 and DnaJC1 resulted in a significantly increased mortality rate in individuals after exposure to trans-anethole, compared with the mortality in controls which were injected with dsRNA of GFP (P < 0.05; Fig. 3). However, knockdown of the remaining five genes (ABCG1, CP19, CP7, CP65, and CP12.5) had negligible effects on mortality under trans-anethole treatment (P > 0.05; Fig. 3).

Fig. 3.

Effects of RNAi targeting selected genes on the mortality in M. persicae exposed to trans-anethole. Individuals injected with GFP dsRNA were the controls. Data were presented as mean ± SE. ‘*’ denotes a significant difference (Student’s t-test, P < 0.05).

Discussion

We constructed a transcriptome dataset from M. persicae exposed to trans-anethole, the most abundant bioactive compound in the essential oil of I. verum. Prior to this study, transcriptomic analyses for M. persicae have been used to investigate the genetic response to insecticides and ultraviolet-B (UV-B) irradiation, and the genes regulating development (Silva et al. 2012, Ji et al. 2016, Meng et al. 2019, Yang et al. 2021). However, none of these studies provided information on genes responding to trans-anethole. Hence, the dataset reported here expands understanding of the molecular mechanisms underlying trans-anethole-regulated gene expression in M. persicae.

Exposure to trans-anethole led to a total of 318 up-regulated and 241 down-regulated genes in M. persicae. Previously, a number of DEGs in M. persicae responding to pirimicarb, imidacloprid, and UV-B irradiation were studied (Silva et al. 2012, Meng et al. 2019, Yang et al. 2021). The number of DEGs varies dramatically in the aphids under different stressors. For example, 559 DEGs were identified in the aphids treated with trans-anethole (this study), whereas exposure to imidacloprid and UV-B resulted in 252 and 758 DEGs, respectively (Meng et al. 2019, Yang et al. 2021). Fewer DEGs were found in pirimicarb-treated M. persicae; there are 7–183 up-regulated genes and 17–78 down-regulated genes in different aphid genotypes (Silva et al. 2012). Among the trans-anethole-induced DEGs, the most up-regulated gene was acyl-CoA synthetase, followed with serine/threonine-protein kinase. In insects, proteins belonging to the acyl-CoA synthetase family can activate fatty acids to acyl-CoA and thus play a role in energy metabolism (Alves-Bezerra et al. 2016). Serine/threonine protein kinase is a master regulator of cellular energy metabolism due to its ability to regulate glucose, lipid, and protein metabolism (Witczak et al. 2008). Our findings indicated that trans-anethole activates energy metabolism pathways in M. persicae. This result is consistent with previously reported data. Many genes related to energy metabolism are up-regulated by pirimicarb in M. persicae (Silva et al. 2012). Although trans-anethole and pirimicarb have different structures, both compounds are toxic to insects and would be expected to activate the detoxificaton pathways. It is possible that detoxification of xenobiotic compounds requires substantial energy. Therefore, energy metabolism-related genes are activated thereafter and play an essential role in maintaining a balance between energy production and consumption (Even et al. 2012). In addition, trans-anethole can be hydroxylated in the larvae of Spodoptera litura (Fabricius) (Lepidoptera: Noctuidae) and Trichoplusia ni (Hübner) (Lepidoptera: Noctuidae) (Passreiter et al. 2004). It is possible that similar transformation mechanisms may exist in M. persicae. The hydroxylation of trans-anethole requires enzymes with oxygenase and/or hydroxylase activity. However, none of the oxygenase and hydroxylase genes were found in the DEGs (Supp Table S2 [online only]). We hypothesize that these genes are constitutively expressed in M. persicae.

Trans-anethole is a plant-derived compound with high toxicity against M. persicae (Li et al. 2017). Exposure to trans-anethole may activate the detoxification/defense pathways in M. persicae, and these pathways may participate in the detoxification of trans-anethole. Therefore, although many DEGs were identified in M. persicae, seven up-regulated genes (two ABC transporters, one DnaJC1, and four cuticle proteins) with potential roles in detoxification and defense were selected and their expression patterns were validated by RT-qPCR. The functions of these genes were also investigated using RNAi.

Two significantly up-regulated genes annotated as ABC transporters were identified in the DEGs. ABC transporters are a family of transmembrane proteins critical for the detoxification of xenobiotic compounds in insects (Dermauw and van Leeuwen 2014). Members from the family can translocate a large range of substrates, particularly insecticides and plant secondary metabolites, across cellular membranes in an ATP-dependent manner (Wu et al. 2019, Rösner and Merzendorfer 2020). In Drosophila melanogaster Meigen (Diptera: Drosophilidae) and T. castaneum, knockdown or knockout of ABC transporter genes resulted in increased sensitivity of insects to different pesticides (Denecke et al. 2017, Rösner and Merzendorfer 2020). Several reports showed that essential oil treatments can enhance the expression levels of multiple ABC transporter genes in insects (Liao et al. 2016, Huang et al. 2018a). Chen et al. (2018) identified 67 ABC transporter genes in the M. persicae genome. In this study, we found that silencing one (ABCG4) of the ABC transporters significantly increased the mortality rate in M. persicae following trans-anethole treatment. This indicated that the gene is involved in the detoxification of the compound. The role of ABCG subfamily genes in insecticide detoxification has been reported in other insect species. For example, knockdown of the ABCG15 gene significantly increased thiamethoxam sensitivity in the cotton aphid Aphis gossypii (Glover) (Hemiptera: Aphididae) (Pan et al. 2020), and knockdown of ABCG3 significantly reduced the survival rate of imidacloprid-treated Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae) (He et al. 2019). In addition, silencing of ABCG6, ABCG9, and ABCG14 significantly increased the chlorantraniliprole susceptibility in Plutella xylostella (L.) (Lepidoptera: Plutellidae) (Shan et al. 2021).

A significant induction of the expression level of a DnaJ homolog subfamily C member 1 gene (DnaJC1) was detected after trans-anethole exposure (Table 3). DnaJ, also called heat shock protein 40, is a member of a family of chaperone proteins (Qiu et al. 2006). The most important function of DnaJ is to promote recovery of misfolded and damaged proteins in cells subjected to environmental stresses and insecticides (Chen et al. 2018). In the honeybee Apis cerana Fabricius (Hymenoptera, Apidae), two DnaJC subfamily genes (DnaJC3 and DnaJC8) are required for the tolerance to stress conditions such as exposure to heat, abamectin, and lambda-cyhalothrin (Li et al. 2018, Zhang et al. 2019, Li et al. 2020). Previous studies demonstrated that trans-anethole can cause abnormally high-level reactive oxygen species (ROS), leading to oxidative damage in insect cells (Shahriari et al. 2018). Here, we observed that knockdown of the DnaJC1 gene increased mortality in M. persicae after trans-anethole exposure. This suggests that up-regulation of DnaJC1 could be an important mechanism to protect against oxidative damage and maintain the cell homeostasis in M. persicae under trans-anethole stress.

In addition to ABC transporter and DnaJ genes, up-regulation of four cuticle protein genes was detected in M. persicae. Cuticle proteins are main components of the insect cuticle, which is the primary protective barrier against xenobiotic compounds and pathogenic microbes (Charles 2010). Enhanced transcription of cuticle protein genes is often observed in insects exposed to insecticides, and this up-regulation suggests the contribution of cuticle proteins to tolerance or resistance to pesticides (David et al. 2014, Gao et al. 2018, Zhou et al. 2018). In the mosquito Culex pipiens L. (Diptera: Culicidae), a cuticle protein enhances pyrethroid resistance by increasing the thickness of the cuticle (Huang et al. 2018b). A silkmoth Bombyx mori (L.) (Lepidoptera: Bombycidae) strain that carries a mutation in a cuticle protein gene is more sensitive to deltamethrin (Xiong et al. 2018). In this study, knockdown of the four cuticle protein genes had negligible effect on the mortality of M. persicae exposed to trans-anethole. These unexpected results suggest that these genes may not be required for protecting against trans-anethole in M. persicae. Rather, these genes may be involved in important pathways such as defense against other insecticides, pathogenic microbes, and transmission of plant viruses (Ortiz-Urquiza and Keyhani 2013, Balabanidou et al. 2018, Deshoux et al. 2018). Further analyses are needed to evaluate the underlying mechanisms.

Cytochrome P450s (CYPs), glutathione S-transferases (GSTs), and esterases are key enzymes in insect metabolism of xenobiotics (Li et al. 2007, Liu et al. 2017, Shi et al. 2018). Unexpectedly, we did not identify up-regulated genes encoding CYP, GST, and esterase, but found that a GST gene (GenBank accession no. XM_022315636) and an esterase gene (XM_022305972) were down-regulated following trans-anethole treatment (Supp Table S2 [online only]). In a number of insect species such as D. melanogaster, S. exigua, S. litura, and T. castaneum, the cap ‘n’ collar isoform-C/muscle aponeurosis fibromatosis (CncC/Maf) and aryl hydrocarbon receptor/aryl hydrocarbon receptor nuclear translocator (AhR/ARNT) signal pathways are responsible for the regulation of detoxification genes such as CYPs, GSTs, and esterases (Palli 2020, Amezian et al. 2021). It is possible that the transcription of the GST and esterase genes in M. persicae are under the influence of CncC/Maf and/or AhR/ARNT pathways and the downregulation is a result of a cascade response caused by disturbing the pathways.

In conclusion, this study provides a comprehensive transcriptomic resource and expression profiles of genes in M. persicae responding to trans-anethole. We discovered numerous DEGs that belong to different biological pathways. We identified a number of genes such as ABC transporters, DnaJ, and cuticle proteins and revealed the involvement of ABCG4 and DnaJC1 in the defensive response to trans-anethole. The findings enhance understanding of the manner in which M. persicae responds to plant essential oils.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.