De novo transcriptome analysis of gene responses to pest feeding in leaves of Panax ginseng C. A. Meyer.

The aim of the present study was to investigate the transcriptomic differences between Panax ginseng [Renshen (RS)] plants bitten by pests (n=3, test group; samples defined as RS11‑13) or not (n=3, control group; samples defined as RS1‑3) using de novo RNA sequencing on an Illumina HiSeq™ 2000 platform. A total of 51,097,386 (99.6%), 49,310,564 (99.5%), 59,192,372 (99.6%), 60,338,540 (99.5%), 56,976,410 (99.6%) and 54,226,588 (99.6%) clean reads were obtained for RS11, RS12, RS13, RS1, RS2 and RS3, respectively. De novo assembly generated 370,267 unigenes, 927 of which were differentially expressed genes (DEGs), including 782 significantly upregulated and 145 significantly downregulated genes. Function enrichment analysis revealed that these DEGs were located in 28 significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways, including phenylpropanoid biosynthesis (for example, TRINITY_DN30766_c0_g2_i1, encoding peroxidase 20) and mitogen‑activated protein kinase (MAPK) signaling (TRINITY_DN85589_c0_g1_i1, encoding WRKY transcription factor 75). Weighted gene co‑expression network analysis identified modules including TRINITY_DN85589_c0_g1_i1, TRINITY_DN58279_c0_g1_i1 [encoding aspartyl protease (AP)] and TRINITY_DN74866_c0_g2_i1 [encoding 12‑oxophytodienoate reductase (OPR)] that may be the most significantly associated with pest responses. In this module, TRINITY_DN85589_c0_g1_i1 may co‑express with TRINITY_DN58279_c0_g1_i1 or TRINITY_DN74866_c0_g2_i1. WRYK and AP have been suggested to promote the activity of antioxidant peroxidase. Collectively, the findings from the present study suggested that a MAPK‑WRKY‑OPR/AP‑peroxidase signaling pathway may be a potentially important mechanism underlying defense responses against pests in ginseng plants.

Introduction

Panax ginseng C. A. Meyer is a popular medicinal plant species grown in northeast China. Previous studies have reported that ginseng exhibits a wide range of pharmacological effects (1), including antifatigue (2), antitumor (3), antioxidant (4), antidiabetic (5), anti-obesity (6) and immunomodulatory (3) effects. Thus, there is notable demand for ginseng products on the market; however, in the wild, ginseng plants are susceptible to attack from a range of native and invasive pests (7), including Locusta migratoria L., Loxostege sticticalis and Xestia c-nigrum, which lead to substantial losses in production and quality. Thus, it is necessary to understand the molecular mechanisms underlying plant-pest interaction, particularly resistance and defense against pest feeding, to optimize the environmental conditions and develop resistant ginseng varieties.

Previous studies have investigated the activity of molecular response mechanisms to pest herbivory in various plants, including plant hormone signal transduction [involving jasmonic acid (JA), ethylene, abscisic acid (AA) and salicylic acid] and transcriptional activation of defense-associated genes [superoxide dismutase (SOD), peroxidase, ascorbate peroxidase (APX), polyphenol oxidase, phenylalanine ammonia lyase, catalase (CAT) and glutathione-S-transferase (GST)] (8,9); however, there is limited information regarding the defense responses of ginseng against pests. In the present study, RNA sequencing (RNA-Seq) was conducted to analyze transcriptomic responses to pest attacks in ginseng plants.

Materials and methods

Plant materials

The 4-year-old Panax ginseng C. A. Meyer was cultivated in the experimental fields of Jilin Agricultural Science and Technology College (44°02‘33.34’ N, 126°06‘22.64’ W). Jilin City is located in a temperate continental monsoon climate, with an annual mean temperature of 5.6°C (high, 22.1°C) and a mean annual rainfall of 679 mm. A total of 6 ginseng plants were included in the study; 3 (RS11, RS12 and RS13; test group) were exposed to feeding by pests (mainly Locusta migratoria L.; Fig. 1), whereas 3 (RS1, RS2 and RS3; control group) were not. Leaves were harvested following exposure to pests for 3–4 days, and three replicates were conducted for each plant to pool the samples. Following cleaning, the leaves were immediately frozen in liquid nitrogen, and stored at −80°C until further use.

Figure 1.

Feeding of pests on Panax ginseng C. A. Meyer. Representative image of Panax ginseng C. A. Meyer following bites from pests, mainly Locusta migratoria L.

RNA isolation and sequencing

Total RNA was extracted from the samples using TRIzol® reagent (Thermo Fisher Scientific, Inc.). The integrity of the total RNA was determined via 1% agarose gel electrophoresis and its concentration was quantified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc.). Total RNA (~1 µg) with RNA Integrity Number ≥8 was used for library construction using a NEBNext® Ultra™ RNA Library Prep kit for Illumina® (New England BioLabs, Inc.) following the manufacturer's protocols: The NEBNext Poly(A) mRNA Magnetic Isolation Module was used for isolation of poly(A) mRNA. mRNA fragmentation and priming was conducted using NEBNext First Strand Synthesis Reaction Buffer and NEBNext Random Primers. Fragmented RNA was reverse transcribed into first-strand cDNA using the ProtoScript II Reverse Transcriptase, and second-strand cDNA was synthesized using the Second Strand Synthesis Enzyme Mix (all New England BioLabs, Inc.). Double-stranded cDNA was purified via AxyPrep Mag PCR Clean-up (Axygen; Corning, Inc.) and then treated with the End Prep Enzyme Mix (New England BioLabs, Inc.) to repair the ends, and attach a dA-tail to one end and adaptors to the two ends. Size selection of adaptor-ligated DNA was also performed using AxyPrep Mag PCR Clean-up, and fragments of ~360 bp were recovered. Then, 11 cycles of PCR amplification were performed using P5 (5′-AGATCGGAAGAGCGTCGTGTAGGGAAAGA-3′) and P7 (5′-GATCGGAAGAGCACACGTCTGAACTCCAGTCACAAGACGGAATCTCGTATGCCGTCTTCTGCTTG-3′) primers with Phusion® Hot Start Flex 2X Master Mix (New England Biolabs, Inc.) under the following thermocycling conditions: 98°C for 10 sec, 60°C for 30 sec, and 72°C for 15 sec, and 72°C for 10 min, to enrich the purified cDNA. The PCR products were cleaned up using AxyPrep Mag PCR Clean-up, validated using an Agilent 2100 Bioanalyzer and quantified using a Qubit 2.0 Fluorometer (Invitrogen; Thermo Fisher Scientific, Inc.). The cDNA library was sequenced by Genewiz, Inc. using an Illumina HiSeq 2000 sequencer (Illumina, Inc.) in 2×150 bp paired-end (PE) mode.

RNA-Seq data analysis

Raw Illumina data were demultiplexed using BCL2FASTQ software (version 2.20; Illumina, Inc.). Raw read quality was determined using FastQC (version 0.10.1; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The reads were pre-processed using Cutadapt (version 1.9.1; http://cutadapt.readthedocs.io/en/stable/) (10) to remove residual adaptor sequences, and reads with low-quality bases (<20 nt in length), N content >10% and length <75 bp following trimming. High-quality clean data in fastq format were assembled de novo to generate the unigene sequence file using the Trinity program (version 2.2.0) with the default parameters (11). Unigenes were annotated via Basic Local Alignment Search Tool against public databases, including non-redundant protein database (nr; http://blast.ncbi.nlm.nih.gov/Blast.cgi), Clusters of Orthologous Groups (COG; http://www.ncbi.nlm.nih.gov/COG/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp).

Bowtie2 version 2.1.0 (12) with the default parameters was used to map the clean reads to the unigenes. RNA-Seq by Expectation Maximization version 1.2.6 (13) was used to estimate the expression levels (fragments per kilobase per million mapped reads) of genes and isoforms from the PE clean data. The differentially expressed genes (DEGs) between the test and control groups were identified using the Bioconductor package DESeq2 (version 1.6.3; http://bioconductor.org/packages/release/bioc/html/DESeq2.html) (14), a model based on a negative binomial distribution. The P-value was adjusted by Benjamini and Hochberg's method (15) to control for false discovery rate (FDR). FDR <0.05 and |log2 fold change (FC)| >1 (FC >2) were set as the threshold value. The underlying functions of DEGs were predicted via KEGG pathway enrichment analyses with a hypergeometric test. An adjusted P-value (Q-value) <0.05 was considered to be statistically significant. In addition, weighted gene co-expression network analysis (16) was performed to identify significant modules of highly associated genes related to pest responses from the DEGs and all unigenes. Highly connected genes may be regarded as hub genes. The top co-expression pairs (weight >0.6) were used to construct the co-expression network using Cytoscape software (version 2.8; www.cytoscape.org/) (17).

Results

Illumina sequence analysis

To determine the global transcriptome profile of ginseng in response to pests, three RNA libraries were constructed and deep RNA-Seq was performed on the leaves of ginseng plants from the test and control groups. A total of 51,323,700, 49,534,966, 59,458,254, 60,617,462, 57,222,860 and 54,469,800 raw reads were generated for the RS11, RS12, RS13, RS1, RS2 and RS3 samples, respectively (Table I). Following quality control to remove the low-quality reads and adaptor sequences, 51,097,386 (99.6%), 49,310,564 (99.5%), 59,192,372 (99.6%), 60,338,540(99.5%), 56,976,410 (99.6%) and 54,226,588 (99.6%) clean reads were retained for further analysis (Table I). Additionally, >98% of the reads exhibited an average quality score of >20 (Q20) and the GC content was consistently ~43% for all samples, suggesting that the sequencing was highly accurate.

Table I.

Quality control results.

A, Raw reads
Samples	Length	Reads	Bases	Q20, %	Q30, %	GC, %	N, ppm
RS11	150	51,323,700	7,698,555,000	97.86	94.98	43.60	413.32
RS12	150	49,534,966	7,430,244,900	97.83	94.93	43.38	410.77
RS13	150	59,458,254	8,918,738,100	97.95	95.14	43.77	406.53
RS1	150	60,617,462	9,092,619,300	97.80	94.85	43.56	406.93
RS2	150	57,222,860	8,583,429,000	97.82	94.89	44.07	406.57
RS3	150	54,469,800	8,170,470,000	97.85	94.97	43.88	410.32
B, Clean reads
Samples	Length	Reads	Bases	Q20, %	Q30, %	GC, %	N, ppm
RS11	148.87	51,097,386	7,606,807,843	98.13	95.32	43.65	8.09
RS12	148.85	49,310,564	7,339,855,636	98.11	95.27	43.43	8.01
RS13	148.80	59,192,372	8,807,982,127	98.22	95.47	43.81	8.09
RS1	148.88	60,338,540	8,983,240,212	98.07	95.20	43.62	8.03
RS2	148.89	56,976,410	8,482,933,397	98.09	95.22	44.12	8.16
RS3	148.87	54,226,588	8,072,497,768	98.12	95.31	43.93	8.08

Q20 and Q30, the percentage of bases with Phred values >20 and >30, respectively; GC content, the GC ratio of the total base number.

De novo assembly of the clean reads produced 11,548,589 contigs of 678,729,555 nucleotides (nt); the average length of these contigs was 58.77 nt, with an N50 of 48 nt. Further assembly of these contigs generated 370,267 unigenes with a mean length and N50 of 626.17 and 839 nt, respectively. A total of 230,086 unigenes (62.14%) were 200–500 nt in length; 83,195 unigenes (22.47%) were 500–1,000 nt; and 56,985 unigenes (15.39%) were >1,000 nt.

Functional annotation results revealed that 200,394 unigenes (54.1%) were annotated to at least one public database. In total, 191,132 unigenes were annotated to the Nr database, among which 62,196 unigenes were identified in Daucus carota subsp. sativus [including TRINITY_DN85589_c0_g1_i1, which may encode WRKY transcription factor 75 (gi|1040876417|ref|XP_017248451.1|); TRINITY_DN74866_c0_g2_i1, which may encode 12-oxophytodienoate reductase (OPR)2-like isoform X1 (gi|1040859474|ref|XP_017240506.1|); and TRINITY_DN30766_c0_g2_i1, which may encode peroxidase 20 (gi|1040813078|ref|XP_017228400.1|); Fig. 2]. Of the 97,892 unigenes that were assigned to the COG database, 11,938 unigenes belonged to the cluster ‘Post-translational modification, protein turnover, chaperones’ [including TRINITY_DN58279_c0_g1_i1, aspartyl protease (AP; KOG1339)], followed by ‘general function prediction only’ [11,648 unigenes, including TRINITY_DN74866_c0_g2_i1, NADH:flavin oxidoreductase/ODR (KOG0134); Fig. 3]. A total of 53,451 unigenes were mapped to 128 KEGG pathways, among which ‘signal transduction’ was the most enriched, featuring 15,940 unigenes, including TRINITY_DN85589_c0_g1_i1 (encoding WRKY transcription factor 33; Fig. 4). WRKY transcription factor 33 was enriched in the ‘MAPK signaling pathway-plant’ (ko04016) and in ‘plant-pathogen interactions’ (ko04626) (data not shown). ‘Carbohydrate metabolism’ was the second most enriched pathway, featuring 8,759 unigenes, including those involves in phenylpropanoid biosynthesis, such as TRINITY_DN30766_c0_g2_i1.

Figure 2.

Functional classification of Panax ginseng unigenes based on NR. NR, non-redundant protein database.

Figure 3.

Functional classification of Panax ginseng unigenes based on COG. COG, Clusters of Orthologous Groups.

Figure 4.

Functional classification of Panax ginseng unigenes based on Kyoto Encyclopedia of Genes and Genomes annotation.

DEG analysis

Comparative transcriptome profiling yielded 927 DEGs, including 782 significantly upregulated and 145 significantly downregulated genes (Table II). The heat map indicated that these DEGs (Fig. 5) could clearly distinguish between the test and control groups. These DEGs were significantly enriched into 28 KEGG pathways (Fig. 6; Table III), including ‘phenylpropanoid biosynthesis’ (TRINITY_DN30766_c0_g2_i1) and ‘MAPK signaling pathway-plant’ (TRINITY_DN85589_c0_g1_i1).

Table II.

Top 20 upregulated and downregulated differentially expressed genes.

Gene ID	log2FC	FDR
TRINITY_DN178238_c0_g1_i1	5.80	1.25×10−103
TRINITY_DN16862_c0_g2_i1	5.12	1.73×10−62
TRINITY_DN91119_c0_g4_i2	4.86	4.09×10−55
TRINITY_DN16862_c0_g1_i1	4.69	1.36×10−50
TRINITY_DN71994_c0_g2_i1	4.66	1.30×10−49
TRINITY_DN101501_c0_g1_i2	4.61	4.50×10−51
TRINITY_DN86652_c0_g1_i1	4.60	3.58×10−51
TRINITY_DN109207_c1_g3_i1	4.44	2.91×10−46
TRINITY_DN87037_c0_g2_i2	4.22	3.23×10−41
TRINITY_DN106559_c0_g1_i1	4.13	5.08×10−38
TRINITY_DN109207_c1_g7_i1	4.04	1.77×10−35
TRINITY_DN101501_c0_g1_i1	3.98	1.90×10−34
TRINITY_DN109207_c1_g4_i1	3.68	5.59×10−28
TRINITY_DN38537_c0_g1_i1	3.54	8.21×10−38
TRINITY_DN109207_c1_g4_i2	3.51	4.92×10−25
TRINITY_DN105289_c0_g1_i2	3.48	1.19×10−24
TRINITY_DN58279_c0_g1_i1	1.59	1.06×10−06
TRINITY_DN30766_c0_g2_i1	1.49	6.01×10−04
TRINITY_DN85589_c0_g1_i1	1.35	5.32×10−03
TRINITY_DN111795_c3_g1_i7	1.045	3.00×10−02
TRINITY_DN118817_c1_g2_i6	−1.66	8.48×10−05
TRINITY_DN119005_c2_g8_i2	−1.68	6.85×10−05
TRINITY_DN117632_c1_g3_i4	−1.68	1.35×10−05
TRINITY_DN94357_c3_g1_i6	−1.69	5.64×10−05
TRINITY_DN119337_c3_g11_i6	−1.71	3.85×10−05
TRINITY_DN108111_c1_g12_i4	−1.76	2.38×10−05
TRINITY_DN113154_c3_g3_i11	−1.85	6.97×10−06
TRINITY_DN117824_c1_g15_i8	−1.87	4.71×10−06
TRINITY_DN100396_c0_g3_i2	−1.90	3.17×10−06
TRINITY_DN111108_c0_g1_i9	−1.91	2.95×10−06
TRINITY_DN116269_c4_g2_i9	−1.95	1.49×10−06
TRINITY_DN119178_c2_g5_i2	−1.95	6.83×10−07
TRINITY_DN114880_c1_g1_i1	−1.99	7.82×10−07
TRINITY_DN113891_c0_g1_i2	−2.12	6.78×10−08
TRINITY_DN118920_c1_g3_i21	−2.20	1.98×10−08
TRINITY_DN114499_c0_g1_i5	−2.22	8.01×10−09
TRINITY_DN116506_c3_g5_i18	−2.30	2.10×10−09
TRINITY_DN117970_c1_g2_i1	−2.40	1.79×10−10
TRINITY_DN117566_c0_g1_i1	−2.82	3.79×10−15
TRINITY_DN113064_c0_g3_i4	−3.55	8.26×10−26

FC, fold change; FDR, false discovery rate.

Figure 5.

Cluster analysis of DEGs in Panax ginseng following pest biting. Cluster analysis of 927 DEGs between test and control Panax ginseng C. A. Meyer leaves following attack by pests (mainly Locusta migratoria L.) based on hierarchical clustering. RS11-13 are test samples (exposure to pest feeding); RS1-3 are control samples. DEG, differentially expressed gene; RS, Renshen.

Figure 6.

KEGG pathway enrichment analysis of gene expression in Panax ginseng following pest feeding. KEGG pathway enrichment analysis of 927 differentially expressed genes between test (exposure to pest feeding) and control Panax ginseng C. A. Meyer leaves following attack by pests (mainly Locusta migratoria L.). KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table III.

KEGG pathways for the differentially expressed genes.

Pathway ID	Pathway	Gene list	Q-value
ko00940	Phenylpropanoid	TRINITY_DN38537_c0_g1_i1, TRINITY_DN119027_c3_g1_i5,	1.40×10−21
	biosynthesis	TRINITY_DN90698_c1_g1_i1, TRINITY_DN93128_c0_g1_i1,
		TRINITY_DN104257_c1_g1_i1, TRINITY_DN104257_c1_g1_i5,
		TRINITY_DN104257_c1_g1_i2, TRINITY_DN117744_c5_g27_i4,
		TRINITY_DN119027_c3_g1_i1, TRINITY_DN111296_c1_g1_i1,
		TRINITY_DN119027_c3_g1_i2, TRINITY_DN119027_c3_g1_i7,
		TRINITY_DN56484_c1_g1_i1, TRINITY_DN110990_c3_g4_i4,
		TRINITY_DN2548_c0_g1_i1, TRINITY_DN110990_c3_g4_i3,
		TRINITY_DN30766_c0_g2_i1, TRINITY_DN102469_c0_g1_i1,
		TRINITY_DN77906_c0_g1_i1, TRINITY_DN111296_c1_g4_i3,
		TRINITY_DN102469_c0_g2_i1, TRINITY_DN119565_c6_g21_i1,
		TRINITY_DN114431_c3_g1_i1, TRINITY_DN116123_c1_g9_i3
ko01040	Biosynthesis of	TRINITY_DN115911_c2_g17_i2, TRINITY_DN70767_c1_g1_i1,	9.31×10−16
	unsaturated fatty acids	TRINITY_DN115911_c2_g17_i5, TRINITY_DN38619_c0_g1_i1,
		TRINITY_DN62336_c0_g1_i1, TRINITY_DN108852_c3_g3_i1,
		TRINITY_DN35750_c0_g1_i1, TRINITY_DN77948_c0_g1_i1,
		TRINITY_DN115911_c2_g17_i7, TRINITY_DN105028_c0_g5_i1,
		TRINITY_DN105028_c0_g7_i1, TRINITY_DN75773_c0_g1_i1,
		TRINITY_DN117863_c6_g4_i1, TRINITY_DN115911_c2_g22_i2,
		TRINITY_DN115911_c2_g34_i1, TRINITY_DN2487_c0_g1_i1,
		TRINITY_DN105028_c0_g1_i1
ko00073	Cutin, suberine and wax	TRINITY_DN106559_c0_g1_i1, TRINITY_DN105289_c0_g1_i2,	1.20×10−12
	biosynthesis	TRINITY_DN110555_c1_g3_i1, TRINITY_DN103944_c0_g2_i1,
		TRINITY_DN108493_c1_g1_i3, TRINITY_DN112895_c0_g1_i2,
		TRINITY_DN112895_c0_g1_i1, TRINITY_DN110677_c1_g1_i2,
		TRINITY_DN110677_c1_g1_i1
ko01212	Fatty acid metabolism	TRINITY_DN112064_c0_g1_i3, TRINITY_DN115911_c2_g17_i2,	1.33×10−10
		TRINITY_DN70767_c1_g1_i1, TRINITY_DN115911_c2_g17_i5,
		TRINITY_DN38619_c0_g1_i1, TRINITY_DN62336_c0_g1_i1,
		TRINITY_DN108852_c3_g3_i1, TRINITY_DN35750_c0_g1_i1,
		TRINITY_DN77948_c0_g1_i1, TRINITY_DN115911_c2_g17_i7,
		TRINITY_DN105028_c0_g5_i1, TRINITY_DN88167_c0_g1_i2,
		TRINITY_DN105028_c0_g7_i1, TRINITY_DN75773_c0_g1_i1,
		TRINITY_DN117863_c6_g4_i1, TRINITY_DN115911_c2_g22_i2,
		TRINITY_DN115911_c2_g34_i1, TRINITY_DN2487_c0_g1_i1,
		TRINITY_DN105028_c0_g1_i1
ko00360	Phenylalanine	TRINITY_DN119027_c3_g1_i5, TRINITY_DN117744_c5_g27_i4,	3.63×10−08
	metabolism	TRINITY_DN119027_c3_g1_i1, TRINITY_DN111296_c1_g1_i1,
		TRINITY_DN119027_c3_g1_i2, TRINITY_DN119027_c3_g1_i7,
		TRINITY_DN59880_c0_g2_i1, TRINITY_DN77906_c0_g1_i1,
		TRINITY_DN111296_c1_g4_i3, TRINITY_DN119565_c6_g21_i1,
		TRINITY_DN114431_c3_g1_i1
ko00941	Flavonoid biosynthesis	TRINITY_DN117744_c5_g27_i4, TRINITY_DN103733_c0_g4_i1,	5.32×10−06
		TRINITY_DN56484_c1_g1_i1, TRINITY_DN25782_c0_g1_i1,
		TRINITY_DN119565_c6_g21_i1
ko00040	Pentose and glucuronate	TRINITY_DN102935_c1_g1_i1, TRINITY_DN105582_c5_g2_i3,	1.54×10−04
	interconversions	TRINITY_DN107139_c0_g1_i2, TRINITY_DN61993_c1_g1_i1,
		TRINITY_DN65703_c1_g1_i1, TRINITY_DN108358_c0_g2_i5,
		TRINITY_DN108358_c0_g2_i1, TRINITY_DN96538_c0_g1_i2,
		TRINITY_DN113543_c1_g1_i1
ko00591	Linoleic acid metabolism	TRINITY_DN164530_c0_g1_i1, TRINITY_DN38928_c0_g1_i1,	1.83×10−04
		TRINITY_DN118754_c5_g7_i4, TRINITY_DN118754_c5_g7_i1
ko04745	Phototransduction-fly	TRINITY_DN116200_c1_g13_i2, TRINITY_DN103329_c0_g1_i3,	3.82×10−04
		TRINITY_DN116190_c4_g2_i4, TRINITY_DN113735_c2_g1_i1,
		TRINITY_DN118958_c1_g1_i1, TRINITY_DN118958_c1_g3_i1
ko00945	Stilbenoid, diarylheptanoid	TRINITY_DN117744_c5_g27_i4, TRINITY_DN56484_c1_g1_i1,	1.26×10−03
	and gingerol biosynthesis	TRINITY_DN119565_c6_g21_i1
ko05412	ARVC	TRINITY_DN116200_c1_g13_i2, TRINITY_DN116190_c4_g2_i4,	1.82×10−03
		TRINITY_DN113735_c2_g1_i1, TRINITY_DN118958_c1_g1_i1,
		TRINITY_DN118958_c1_g3_i1
ko05414	Dilated cardiomyopathy	TRINITY_DN116200_c1_g13_i2, TRINITY_DN116190_c4_g2_i4,	1.82×10−03
		TRINITY_DN113735_c2_g1_i1, TRINITY_DN118958_c1_g1_i1,
		TRINITY_DN118958_c1_g3_i1
ko04210	Apoptosis	TRINITY_DN113765_c2_g1_i4, TRINITY_DN80010_c0_g2_i1,	8.98×10−03
		TRINITY_DN109158_c4_g3_i4, TRINITY_DN116200_c1_g13_i2,
		TRINITY_DN116190_c4_g2_i4, TRINITY_DN113735_c2_g1_i1,
		TRINITY_DN118958_c1_g1_i1, TRINITY_DN111795_c3_g1_i7,
		TRINITY_DN107791_c3_g1_i5, TRINITY_DN118958_c1_g3_i1
ko00592	α-linolenic acid	TRINITY_DN164530_c0_g1_i1, TRINITY_DN38928_c0_g1_i1,	9.96×10−03
	metabolism	TRINITY_DN118754_c5_g7_i4, TRINITY_DN118754_c5_g7_i1,
		TRINITY_DN61851_c1_g2_i1
ko00908	Zeatin biosynthesis	TRINITY_DN116264_c1_g3_i5, TRINITY_DN95422_c1_g1_i2	1.32×10−02
ko05130	Pathogenic	TRINITY_DN113765_c2_g1_i4, TRINITY_DN80010_c0_g2_i1,	1.52×10−02
	Escherichia coli infection	TRINITY_DN109158_c4_g3_i4, TRINITY_DN116200_c1_g13_i2,
		TRINITY_DN116190_c4_g2_i4, TRINITY_DN113735_c2_g1_i1,
		TRINITY_DN118958_c1_g1_i1, TRINITY_DN118958_c1_g3_i1
ko05410	HCM	TRINITY_DN116200_c1_g13_i2,TRINITY_DN116190_c4_g2_i4,	1.52×10−02
		TRINITY_DN113735_c2_g1_i1,TRINITY_DN118958_c1_g1_i1,
		TRINITY_DN118958_c1_g3_i1,
ko02024	Quorum sensing	TRINITY_DN112064_c0_g1_i3, TRINITY_DN102935_c1_g1_i1,	1.52×10−02
		TRINITY_DN61993_c1_g1_i1, TRINITY_DN108358_c0_g2_i5,
		TRINITY_DN108358_c0_g2_i1, TRINITY_DN97205_c0_g2_i1,
		TRINITY_DN96538_c0_g1_i2
ko05416	Viral myocarditis	TRINITY_DN116200_c1_g13_i2, TRINITY_DN116190_c4_g2_i4,	1.78×10−02
		TRINITY_DN113735_c2_g1_i1, TRINITY_DN118958_c1_g1_i1,
		TRINITY_DN107791_c3_g1_i5, TRINITY_DN118958_c1_g3_i1
ko05133	Pertussis	TRINITY_DN112623_c3_g1_i4, TRINITY_DN106201_c1_g2_i7,	2.30×10−02
		TRINITY_DN103329_c0_g1_i3, TRINITY_DN75555_c1_g1_i1,
		TRINITY_DN31358_c0_g2_i1, TRINITY_DN109488_c1_g1_i11,
		TRINITY_DN111795_c3_g1_i7, TRINITY_DN109488_c1_g1_i6,
		TRINITY_DN88066_c0_g1_i1
ko00904	Diterpenoid biosynthesis	TRINITY_DN105174_c1_g1_i1, TRINITY_DN108094_c1_g2_i1	2.47×10−02
ko04611	Platelet activation	TRINITY_DN116200_c1_g13_i2, TRINITY_DN116190_c4_g2_i4,	2.47×10−02
		TRINITY_DN113735_c2_g1_i1, TRINITY_DN118958_c1_g1_i1,
		TRINITY_DN111795_c3_g1_i7, TRINITY_DN118958_c1_g3_i1
ko05140	Leishmaniasis	TRINITY_DN112623_c3_g1_i4, TRINITY_DN106201_c1_g2_i7,	2.77×10−02
		TRINITY_DN75555_c1_g1_i1, TRINITY_DN31358_c0_g2_i1,
		TRINITY_DN109488_c1_g1_i11, TRINITY_DN111795_c3_g1_i7,
		TRINITY_DN109488_c1_g1_i6, TRINITY_DN88066_c0_g1_i1
ko00052	Galactose metabolism	TRINITY_DN114499_c0_g1_i5, TRINITY_DN112361_c0_g1_i1,	3.33×10−02
		TRINITY_DN105582_c5_g2_i3, TRINITY_DN114499_c0_g1_i7,
		TRINITY_DN114401_c0_g1_i7, TRINITY_DN116821_c1_g3_i6
ko04614	Renin-angiotensin system	TRINITY_DN113178_c1_g1_i7, TRINITY_DN118282_c3_g5_i8	3.33×10−02
ko04016	MAPK signaling	TRINITY_DN101657_c1_g1_i1, TRINITY_DN115200_c1_g1_i3,	3.33×10−02
	pathway-plant	TRINITY_DN103013_c0_g1_i1, TRINITY_DN103329_c0_g1_i3,
		TRINITY_DN113050_c3_g1_i3, TRINITY_DN85589_c0_g1_i1,
		TRINITY_DN87820_c0_g1_i2
ko04670	Leukocyte	TRINITY_DN116200_c1_g13_i2, TRINITY_DN116190_c4_g2_i4,	3.71×10−02
	transendothelial	TRINITY_DN113735_c2_g1_i1, TRINITY_DN118958_c1_g1_i1,
	migration	TRINITY_DN118958_c1_g3_i1
ko01524	Platinum drug resistance	TRINITY_DN104989_c0_g2_i1, TRINITY_DN80125_c0_g2_i1,	4.25×10−02
		TRINITY_DN111271_c1_g1_i3, TRINITY_DN111271_c1_g1_i3,
		TRINITY_DN111795_c3_g1_i7, TRINITY_DN107791_c3_g1_i5

ARVC, arrhythmogenic right ventricular cardiomyopathy; HCM, hypertrophic cardiomyopathy; MAPK, mitogen-activated protein kinase; Q-value, adjusted P-value.

To determine the important genes involved in pest responses, a co-expression network was constructed. DEGs with similar patterns of expression were grouped into four modules via hierarchical average linkage clustering (Fig. 7A). The turquoise module may be the most strongly associated with pest responses, as the majority of the top 50 genes (including TRINITY_DN30766_c0_g2_i1 and TRINITY_DN85589_c0_g1_i1) in the co-expression network (weight >0.74; Fig. 7B) were located in the turquoise module. TRINITY_DN85589_c0_g1_i1 co-expressed with TRINITY_DN58279_c0_g1_i1 or TRINITY_DN74866_c0_g2_i1 (data not shown) in co-expression networks for DEGs or all unigenes.

Figure 7.

Weighted gene co-expression network analysis of gene expression in Panax ginseng following pest feeding. (A) Eigengene dendrogram of meta-module clusters. The multi-colored bar below the dendrogram indicates the significant modules identified. (B) Gene co-expression networks composed of the top 50 genes (weight >0.74) in the turquoise and blue modules.

Discussion

In the present study, the genetic response to pest bites in the leaves of Jilin ginseng plants was sequenced. Analysis of gene sequencing and expression revealed that activation of the MAPK pathway via the upregulation of WRKY transcription factors, and the co-expression of AP or ODR, may be an important mechanism in response to pest stress in ginseng plants.

Increasing evidence has indicated that WRKY transcription factors are expressed in response to various types of stress, including salt (18), drought (19), heat (20), abscisic acid (19), salicylic acid (21), pathogens (21–23) and herbivores (24,25). Overexpression of WRKY is reported to increase the transcription of antioxidant enzyme genes, including APX, CAT, GST and SOD to reduce reactive oxygen species content, positively regulating plant responses to stress (26,27) and suppressing leaf senescence (28). Furthermore, a series of other stress-associated genes, including cold-regulated 15a (COR15A), COR15B, COR413, COR6.6 (29), OsFRDL4 (30), DgNCED3A, DgNCED3B, DgP5CS, DgCSD1 and DgCSD2 (31) were also reported to be upregulated in WRKY transgenic plants compared with in wild-type plants. Consistent with these studies, it was demonstrated in the present study that WRKY75 and WRKY33 may be significantly upregulated in the leaves of ginseng plants in response to herbivore bites.

The roles of WRKY75 or 33 remain unclear; however, the present study predicted that WRKY may be regulated by upstream MAPKs and interact downstream with ODR or AP. In a previous study, chromatin immunoprecipitation assays revealed that WRKY33 is a substrate of MAPK3/MAPK6, two pathogen-responsive MAPKs, involved in the induction of phytoalexin camalexin production in Arabidopsis thaliana (32). Additionally, mutations of the MAPK3/MAPK6 phosphorylation sites in WRKY33 reduces its ability to promote camalexin induction (32). Similarly, the levels of WRKY53 transcription were reported to be positively regulated by MAPK3/MAPK6 in rice (33). Adachi et al (34) revealed that WRKY transcription factors functioned as substrates of the MAPK kinase 2/salicylic acid-induced protein kinase/wound-induced protein kinase signaling cascade. In the present study, MAPK10 (TRINITY_DN111795_c3_g1_i7) was demonstrated to be significantly upregulated, indicating that a MAPK10/WRKY33 signaling cascade may be a mechanism underlying pest responses in ginseng plants.

OPRs belong to a family of flavin-dependent oxidoreductases. OPRs are reported to convert 12-oxophytodienoate into 12-oxophytoenoic acid and participate in the biosynthesis of JA from linolenic acid via the Vick-Zimmerman pathway (35–37). JA is considered to be a signaling molecule involved in stress responses to wounds and herbivore infestation; Xin et al (38) reported that OPR3 was highly expressed in the leaves of Camellia sinensis (L.) exposed to Ectropis obliqua Prout, accompanied by increased JA levels. AP was also demonstrated to regulate fungal and osmotic stress responses in plants (39–42). Overexpression of AP may lead to increased AA levels and promote the activities of various antioxidants, inducing protective autophagy and conferring resistance to stress (40,43). Consistent with these studies, AP was also found to be upregulated in the present study. Thus, it is hypothesized that OPRs and AP may be important downstream targets of WRKY33 in ginseng plants during pest responses.

In addition to the MAPK pathway, the phenylpropanoid biosynthesis pathway was also identified to be significantly enriched with DEGs, including peroxidase 20. Peroxidase is an important antioxidant for pest resistance (44). WRKY (45) and AP (40) promote the transcription of peroxidase; the MAPK and phenylpropanoid biosynthesis pathways may be associated regulatory mechanisms underlying pest resistance.

There are certain limitations to the present study. The sample size was small, which may explain the lack of statistical significance observed for the expression of OPRs. Furthermore, additional experiments (i.e., silencing, PCR, determination of enzyme activity and hormone level detection) (8) are required to validate the importance of the hypothesized MAPK10-WRKY33-OPR/AP-peroxidase 20 pathway or pathways, and their dependence on JA or AA in the pest resistance of ginseng plants.

To the best of our knowledge, the present study reports the first investigation of the transcriptional responses of Panax ginseng C. A. Meyer to pest feeding. The findings suggested that MAPK10-WRKY33-OPR/AP-peroxidase 20 signaling may be a important mechanism underlying defense responses against pests. Further experiments should be conducted to support these conclusions and increase our understanding of plant resistance to pest feeding.