Transcriptome and co-expression network analysis uncover the key genes mediated by endogenous defense hormones in tea plant in response to the infestation of Empoasca onukii Matsuda

The tea plant [Camellia sinensis (L.) O. Kuntze] is an evergreen woody cash crop native to southwest China. Its leaves, especially tender leaves, are the main raw materials for tea, the world's most popular beverage (with the exception of water)^[1]. As a thermophilic, moisture-loving and shade-tolerant plant, the warm and humid growing environment makes it vulnerable to pathogens and herbivorous pests, which seriously threaten the yield and quality of tea^[2]. Effective prevention of losses caused by biotic stress in tea production is a significant problem that needs to be solved.

The tea green leafhopper [Empoasca (Matsumurasca) onukii Matsuda] is one of the most dominant pests in major tea production regions of China, which usually causes 10% to 30% of tea loss in summer and autumn, and more than 50% in some years^[3,4]. As a type of piercing-sucking pest, the tea green leafhopper mainly damages tea plants by sucking the sap from the back of new shoots and tender stems^[5]. Both adults and larvae can cause harm, and the symptoms of damage include yellowing of the leaf margins, reddening of the leaf veins and curling of the leaves. In severely infected plants, the shoot tips turn reddish-brown and scorched, the buds and leaves exhibit stagnant growth and abscission, the whole leaf easily falls off, the leaf texture becomes brittle, and the stem internodes become shorter^[6]. In addition, female adults lay eggs in the tender shoots of tea plants, resulting in damage to the tender shoot conduction tissue and insufficient water and nutrient supply^[5]. Because of the overlapping generations of leafhoppers in tea plantations and the eggs laid on the tender stems of tea plants, the control effect of chemical pesticides is not good^[6,7], and the excessive use of chemical pesticides can cause potential environmental pollution and pesticide residue problems. Given the differences in the resistance of different tea plant cultivars to tea green leafhoppers^[8,9], it is possible to breed resistant cultivars to control tea green leafhoppers. However, crop breeding, particularly for woody cash crops, is a prolonged process. The lack of understanding of the molecular regulatory mechanism of the tea plant defense response to tea green leafhoppers infection also limits the application of bioengineering technology to accelerate the genetic improvement of tea plants.

During the study, it was found that indirect resistance mediated by the special metabolites of tea plants itself plays a key role in the response of tea plants to herbivores^[10]. For instance, in tea plants, the infestation of tea green leafhoppers increases the release of a large number of volatiles, such as (Z)-3-hexenyl acetate, (E)-β-ocimene, (E)-4,8-dimethyl-1,3,7-nonatriene (DMNT), (E,E)-α-farnesene, (E)-nerolidol, methyl salicylate (MeSA), and indole^[4,11]. The volatile composition, defense-related enzyme activity, and gene expression of MeSA-treated tea shoots were similar to those induced by tea green leafhoppers feeding^[11,12]. Exogenous application of MeSA in tea plantations obviously reduced the populations of tea green leafhopper, while the populations of its natural enemies increased significantly^[13]. Therefore, an increasing number of studies have focused on its potential defense mechanisms. Currently, the roles of cis-3-hexenol^[14,15], α-farnesene^[16], β-ocimene^[16−18], nerolidol^[19], linalool^[20], DMNT^[21], and indole^[22] in tea plant defense against insect pests have been elucidated, and their synthase genes were significantly induced to be expressed in insect pests, and most of them could mediate the defense response of neighboring undamaged plants by activating the jasmonic acid (JA) pathway. In turn, when tea plant leaves were treated with JA, the expression of the synthetic genes for the above volatiles were significantly induced, which promoted their accumulation^[16,19]. In addition, feeding by different types of pests can also directly trigger different hormone-mediated signaling pathways, such as JA, salicylic acid (SA), abscisic acid (ABA), ethylene or auxin in damaged tea plants^{[10, 23−25]}. Exogenous application of JA, ethylene or auxin induces the biosynthesis of (+)-catechin (C), epicatechin (EC) and epigallocatechin gallate (EGCG), which in turn enhances the resistance of tea plants to Ectropis grisescens^[25]. Exogenous gallic acid (GA) improves the resistance of tea plants to E. obliqua larvae similarly through the activation of JA signaling^[26]. Volatiles emitted from plants infested with E. obliqua could induce the accumulation of JA, jasmonoyl-isoleucine (JA-Ile) and SA in the neighboring tea plants, and the expression of key genes in the JA synthesis pathway and the PI gene related to direct defense were also significantly up-regulated^[21]. These results suggest that hormone-mediated signaling pathways may be centrally located in the secondary metabolite-mediated defense responses of tea plants against herbivorous pests. However, although this phenomenon has been found in different studies^[6,7,27], the mechanism of the defense hormones-mediated tea plants response to E. onukii damage is unclear. As an analysis of this mechanism is crucial for formulating effective prevention and control strategies, more focused exploration is needed.

Here, based on transcriptome sequencing, we systematically analyzed the differentially expressed genes at different time points (6, 24 and 48 h) after feeding the E. onukii with the same treatment of non-feeding leaves at the same time point as a control, and then focused on the expression patterns of key genes in the main secondary metabolic pathways of tea plants. After that, we clarified the content changes of the main defense related phytohormones (JA, JA-Ile, ABA and SA) in tea plants in this process, and constructed two gene regulatory networks closely related to the changes in hormone content by weighted gene co-expression network analysis (WGCNA), and identified several key hub genes. The results obtained herein provide new insights into the defense response of tea plants to E. onukii infection, and also provide useful genetic resources for the subsequent elucidation of defense mechanisms and breeding resistant cultivars.

In this study, two-year-old potted seedlings of Camellia sinensis cv. 'Qiancha 1' (three plants per pot) grown in the tea germplasm resource nursery of the Tea Research Institute of the Guizhou Academy of Agricultural Sciences (26°30′ N, 106°39′ E) were used as materials. The tea green leafhoppers (Empoasca onukii Matsuda) used in this study were also captured from the tea germplasm resource nursery. For the E. onukii feeding experiments, healthy potted seedlings were moved to the greenhouse, and 24 pots were placed in a 2 m × 1 m × 1 m dense net shed. The temperature was set at 25 ± 2 °C and the photoperiod was set at 16 h of light with 8 h of darkness. After a week of adaptive growth of the tea seedlings in dense nets, approximately 200 adults of E. onukii were placed into it for treatment (M). In another room, under the same conditions as the tea seedlings, E. onukii were not used as a control (CK). When the leaves showed signs (conspicuous brown spots appeared near the middle vein) of E. onukii infection, dense nets and E. onukii were removed and the timing began (0 h). Samples were taken from the damaged young leaves at 6, 24 and 48 h (M6, M24, and M48), and the control samples were taken simultaneously (CK6, CK24, and CK48). Each sample was obtained from at least three tea seedlings as a repetition, and three repetitions were performed at each time point. The samples were frozen in liquid nitrogen and stored at −80 °C for further use.

RNA sequencing and transcriptome data analysis

Transcriptome sequencing and data analysis were performed as described in our previous study^[28]. Briefly, total RNA was extracted using the TaKaRa MiniBEST Plant RNA Extraction Kit (TaKaRa) according to the manufacturer's instructions. Then the genomic DNA was removed using DNase I (TaKaRa). The RNA concentration and purity were measured using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). The total RNA qualified by quality inspection was sent to Biomarker Technology Co. Ltd (Beijing, China) for cDNA library construction and transcriptome sequencing. A total of 18 libraries were constructed and sequenced using the Illumina NovaSeq 6000 platform with 150 bp paired-end reads. The raw sequence data have been deposited in the Genome Sequence Archive^[29] at the National Genomics Data Center, China National Center for Bioinformation, under accession number CRA008188, which is publicly accessible at https://bigd.big.ac.cn/gsa.

The clean reads were mapped to the 'Shuchazao' tea plant reference genome^[30] using HISAT2 tools with default parameters^[31]. Gene expression levels were estimated as fragments per kilobase of transcript per million fragments mapped (FPKM) using StringTie^[31]. Differential expression analysis was performed using DESeq2^[32], and the genes with |log2 (fold change)| ≥1 and adjusted P-value < 0.01 were assigned as differentially expressed genes (DEGs). The gene function was annotated based on public databases, including NR, Pfam, KOG, COG, Swiss-Prot, eggNOG, KEGG and GO. The GOseq R packages and KOBAS software were used for GO and KEGG pathways enrichment analysis of the DEGs, respectively. The PlantTFDB database (http://planttfdb.gao-lab.org) was used to predict transcription factor (TF) genes^[33].

Phytohormones extraction and quantitative analysis

The extraction and quantitation of endogenous phytohormones ABA, JA, JA-Ile and SA in tea leaves were performed according to the method described by Yang et al.^[23] with minor modifications. Firstly, 50−100 mg (fresh weight, FW) of tea leaves samples fully ground with liquid nitrogen were weighed into a 2 mL pre-cooled centrifuge tube, and then 1 mL of ethyl acetate extract containing 5 ng of ²H-SA, ²H-ABA, D6-JA-Ile and 20 ng of D6-JA internal standard were added and vortexed for 15 min. After centrifugation at 12,000 r/min for 15 min at 4 °C, the supernatant was transferred to a new 2 mL centrifuge tube and placed in a rotary evaporator to evaporate the extract. The dried residues were diluted with 200 μL of 50% methanol and vortexed again for 15 min then centrifuged at 13,000 r/min, 4 °C for 20 min. The supernatant was transferred to a 1.5 mL centrifuge tube, and then centrifuged at 13,500 r/min for 30 min at 4 °C. Finally, 100 μL of supernatant was pipetted for sample detection. Phytohormones were measured using an HPLC-MS/MS system equipped with a C₁₈ analytical column (ProntoSIL, 50 mm × 2 mm, 5 μm) and quantified based on the ratio between their ion intensity and their respective internal standards^[23].

Co-expression regulatory network construction by WGCNA

The weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA R package^[34]. The genes with average FPKM > 1 and coefficient of variation (cv) > 0.5 in 18 samples were used to construct the co-expression network. A power of β = 14 was selected as the best soft threshold with 'unsigned' topographical overlap matrix (TOM) types. The minModuleSize was set to 30, mergeCutHeight was set to 0.25, and the other parameters were set to the default settings for module detection. After calculating the eigengene value of each module, the phytohormones phenotype data were introduced for module-trait associations. The gene co-expression networks were visualized using Cytoscape 3.8.2 software.

Quantitative real-time PCR (qRT-PCR) analysis

For qRT-PCR assay, total RNA was used to synthesize the first strand cDNA using PrimeScriptTM RT reagent Kit with gDNA Eraser Perfect Real Time (TaKaRa) according to the manufacturer's instructions. The Primer Premier 5 software was used to design the gene specifc primers, and the TB Green Premix Ex Taq™ II (TaKaRa) was used for qRT-PCR reactions on a BioRad CFX96 Real-Time PCR System (BioRad). PCR amplification was performed in 10 μL reaction mixtures, which contained 5 μL of TB Green Premix Ex Taq, 1 μL of template cDNA, 0.3 μL of each primer (10 μM) and 3.4 μL of sterile double distilled water. The reaction conditions were as follows: 95 °C for 3 min; 40 cycles of 95 °C for 30 s, 57 °C for 30 s, and 72 °C for 20 s. Each sample was tested with three biological replicates. The CsGAPDH gene (CSS0050274) was used as an internal reference and gene expression levels were calculated using the 2^−ΔΔCᴛ method^[35]. The primers used in this study are listed in Supplemental Table S1.

A total of 18 libraries were constructed for transcriptome sequencing, and 126.40 Gb of clean data were obtained, with an average of 5.83 Gb per sample. The Q30 value ranged from 91.00% to 93.97%, averaging 93.16%. Approximately 84.22%−87.40% of the clean reads were successfully mapped to the tea plant reference genome (Supplemental Table S2). In addition, a total of 13,087 new genes were predicted, and 7,789 were functionally annotated using the public database (NR, Pfam, KOG, COG, Swiss-Prot, eggNOG, KEGG and GO). Finally, a total of 56,141 annotated genes were used for the downstream analysis.

As shown in Fig. 1a, 4,533, 5,499 and 8,643 DEGs were identified at 6 h (2,768 up-regulated and 1,765 down-regulated), 24 h (2,266 up-regulated and 3,233 down-regulated) and 48 h (5,775 up-regulated and 2,868 down-regulated) after feeding, respectively. Except for 24 h, the number of up-regulated DEGs at 6 h and 48 h was greater than that of down-regulated DEGs, especially at 48 h, and the number of up-regulated DEGs was more than twice that of down-regulated. We focused on the up-regulated DEGs and found that 624, 843 and 3,477 were specifically expressed at the three time points. At the same time, we found that 695 up-regulated DEGs overlapped at all the three time points (Fig. 1b).

Figure 1. 

Differentially expressed genes at different time points after E. onukii Matsuda feeding. (a) Differentially expressed gene numbers between different comparison groups (CK6-vs-M6, CK24-vs-M24 and CK48-vs-M48). (b) Venn diagram of the number of up-regulated genes in the three comparison groups. Top 15 KEGG enrichment pathways of up-regulated genes in comparison group CK6-vs-M6 (c), CK24-vs-M24 (d) and CK48-vs-M48 (e). Top 15 (f) KEGG and (g) GO enrichment pathways of co-upregulated genes in the three comparison groups.

The enrichment of KEGG pathways of the up-regulated DEGs showed that 'biosynthesis of secondary metabolites' was most significantly enriched at the three time points (Fig. 1c−e). In addition, 'flavonoid biosynthesis', 'α-linolenic acid metabolism', 'phenylpropanoid biosynthesis', 'monoterpenoid biosynthesis', 'plant hormone signal transduction' and different amino acid biosynthesis pathways were also significantly enriched, but the degree of enrichment was different at different time points. For the up-regulated DEGs overlapped at the three time points, they were also mainly enriched in the above KEGG pathways, while GO functional enrichment showed that they were mainly (top five) enriched in the 'regulation of jasmonic acid mediated signaling pathway' (GO:2000022), 'response to wounding' (GO:0009611), 'regulation of defense response' (GO:0031347), 'response to hormone' (GO:0009725) and 'response to endogenous stimulus' (GO:0009719) biological processes (Fig. 1f, g).

Expression patterns of key structural genes in catechins, caffeine and amino acid synthesis pathways

In tea plants, the synthesis of catechins mainly involves the phenylpropanoid and flavonoid pathways. A total of 46 unigenes encoding 14 enzymes in the phenylpropanoid and flavonoid biosynthesis pathways were identified. As shown in Fig. 2a, the transcriptional expression of most members was significantly induced in the leaves fed with E. onukii, especially at 6 h and 48 h after feeding. Many genes were up-regulated more than 10-fold, such as phenylalanine ammonialyase (PAL:CSS0041448, CSS0001442, CSS0018870 and CSS0021474), cinnamate 4-hydroxylase (C4H:CSS0002737), p-coumaroyl-CoA ligase (4CL:CSS0011741 and CSS0006766), flavanone 3-hydroxylase (F3H:CSS0016177 and CSS0019002), flavonoid 3'-hydroxylase (F3'H:CSS0048905 and CSS0030176) and flavonoid 3',5'-hydroxylase (F3'5'H: CSS0018419), while chalcone synthase (CHS), favonol synthase (FLS), leucoanthocyanidin reductase (LAR) and anthocyanidin synthase (ANS) were up-regulated more than 20 times. We identified five N-methyltransferase (NMT) genes for the caffeine synthesis pathway. Their expression was also induced by E. onukii feeding, especially CSS0035688 and CSS0032602 (tea caffeine synthase, TCS), which were up-regulated by 46-fold and 35-fold at 6 h after feeding, respectively (Fig. 2b).

Figure 2. 

Expression characteristics of key genes in (a) catechins and (b) caffeine synthesis pathways at different time points after E. onukii Matsuda feeding. Expression level is based on the log2 value of FPKM. Red and blue indicate higher or lower expression levels of genes, respectively. The heatmaps are drawn using TBtools software, and the original data of the drawing is shown in Supplemental Table S3. PAL, phenylalanine ammonialyase; C4H, cinnamate 4-hydroxylase; 4CL, p-coumaroyl-CoA ligase; CHS, chalcone synthase; CHI, chalcone isomerase; F3H, flavanone 3-hydroxylase; F3'H, flavonoid 3'-hydroxylase; F3'5'H, flavonoid 3',5'-hydroxylase; DFR, dihydroflavonol 4-reductase; FLS, flavonol synthase; ANS, anthocyanidin synthase; LAR, leucoanthocyanidin reductase; ANR, anthocyanidin reductase; SCPL1A, type 1A serine carboxypeptidase-like acyltransferases. NMT, N-methyltransferase; 7-NMT, xanthosine methyltransferase; MXMT, 7-methylxanthine methyltransferase/theobromine synthase; TCS, tea caffeine synthase/3,7-dimethylxanthine methyltransferase; SAMS, S-adenosylmethionine synthetase; SAH, S-adenosylhomocysteine.

For the synthesis pathways of free amino acids, we divided them into three categories according to the properties of amino acids: theanine synthesis-related pathway (Fig. 3a), aromatic amino acid synthesis pathway (Fig. 3b), and protein-coding amino acid synthesis pathway (Fig. 3c). We found that the glutamate decarboxylase gene (GAD) involved in GABA synthesis and most of the key genes such as 3-deoxy-D-arabinoheptulosonate 7-phosphate synthase (DAHPS), chorismate mutase (CM), prephenate dehydratase (PDT), prephenate aminotransferase (PPAT), and arogenate dehydrogenase (ADH) in the aromatic amino acid synthesis pathway were obviously up-regulated after feeding (Fig. 3a, b). The theanine synthase gene (CsTSI, CSS0049154) was also up-regulated by 1.7 times at 6 h after feeding, and then returned to normal. In addition, the key genes involved in Asp, Lys, Thr, Met and Ile synthesis and metabolism also showed different degrees of up-regulation (Fig. 3c). These results indicate that GABA, Phe, Tyr, Asp, Lys, Thr, Met and Ile also play special roles in the defense response of tea plants to E. onukii.

Figure 3. 

Expression characteristics of key genes in amino acid metabolic pathways at different time points after E. onukii Matsuda feeding. (a) Expression patterns of key genes for theanine, proline, GABA, glutamic acid and arginine synthesis. (b) Expression patterns of genes related to the synthesis of aromatic amino acids. (c) Expression patterns of proteinaceous amino acid synthesis-related genes. Expression level is based on the log2 value of FPKM. Red and blue indicate higher or lower expression levels of genes, respectively. The heatmaps are drawn using TBtools software, and the original data of the drawing is shown in Supplemental Table S3.

Expression patterns of key structural genes in terpene synthesis and α-linolenic acid metabolism pathways

Terpene volatiles (monoterpenes and sesquiterpenes) and green leaf volatiles (GLVs) are mainly synthesized through the terpene synthesis pathway and α-linolenic acid metabolism pathway, respectively. As shown in Fig. 4, most of the key genes of the terpene synthesis pathway were significantly up-regulated after feeding with the E. onukii, including acetoacetyl-CoA thiolase (AACT), 3-hydroxy-3-methylglutaryl-CoA synthase (HMGS), 3-hydroxy-3-methylglutaryl-CoA reductase (HMGR), mevalonate kinase (MK), phospho-mevalonate kinase (PMK), and mevalonate-5-diphosphate decarboxylase (MPDC) involved in sesquiterpenes synthesis and 1-deoxy-D-xylulose 5-phosphate synthase (DXS), 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase (CMK), (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate synthase (HDS), and geranyl diphosphate synthase (GPPS), which are involved in monoterpene synthesis. Although the degree of induced expression of most members was most dramatic at 6 h after feeding, the induction trend continued until 48 h after feeding. In addition, for terpene synthase (TPS), the direct gene of terpene synthesis, we identified 15 differentially expressed members with full-length sequences. Among them, except for CsTPS07, CsTPS15, CsTPS42 and CsTPS43 being significantly down-regulated in CK24-vs-M24, the other members were significantly up-regulated at all the three time points. In particular, CsTPS10, CsTPS08 and CsTPS30 were upregulated more than 100-fold compared to the control at 6, 24, and 48 h after feeding, respectively.

Figure 4. 

Expression patterns of key genes in the terpene synthesis pathway at different time points after E. onukii Matsuda feeding. Expression level is based on the log2 value of FPKM. Red and blue indicate higher or lower expression levels of genes, respectively. The heatmaps are drawn using TBtools software, and the original data of the drawing is shown in Supplemental Table S3. AACT, acetoacetyl-CoA thiolase; HMGS, 3-hydroxy-3-methylglutaryl-CoA synthase; HMGR, 3-hydroxy-3-methylglutaryl-CoA reductase; MK, mevalonate kinase; PMK, phospho-mevalonate kinase; MPDC, mevalonate-5-diphosphate decarboxylase; IDI, isopentenyl diphosphate isomerase; FPPS, farnesyl diphosphate synthase; GPPS, geranyl diphosphate synthase; TPS, terpene synthase; DXS, 1-deoxy-D-xylulose 5-phosphate synthase; DXR, 1-deoxy-D-xylulose 5-phosphate reductase; MCT, 2-C-methyl-D-erythritol 4-phosphate (MEP) cytidyltransferase; CMK, 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase; MDS, 2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase; HDS, (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate synthase; HDR, (E)-4-hydroxy-3-methyl-but-2-enyl diphosphate reductase.

For the α-linolenic acid metabolism pathway, as shown in Fig. 5, the expression of all key genes was dramatically induced by feeding, except for the ADH gene, whose expression was down-regulated at 24 h after feeding compared to the control. In particular, the expression of the acetyl CoA:(Z)-3-hexen-1-ol acetyltransferase gene (CsCHAT1) was up-regulated by more than 400 times at 6 h after feeding. We demonstrated its activity in converting cis-3-hexenol into cis-3-hexenyl acetate (unpublished data), indicating its specific role in mediating tea plant defense against the E. onukii. In addition, the α-linolenic acid metabolism pathway is also the key pathway for plants to synthesize the defense hormone JA. The up-regulated expression of its upstream synthetic genes allene oxide synthase (AOS), allene oxide cyclase (AOC), OPDA reductase 3 (OPR3), OPC-8:0 CoA ligase 1 (OPCL1), acyl-CoA oxidase (ACX), and multifunctional protein 2 (MFP2) indicated that the feeding of the E. onukii activated JA synthesis in tea plants.

Figure 5. 

Expression characteristics of key genes in α-linolenic acid metabolism pathway at different time points after E. onukii Matsuda feeding. (a) The α-linolenic acid metabolism pathway and its key genes in tea plant. (b) Heatmap of key gene expression patterns in α-linolenic acid metabolism pathway. The number on the heat map represents the average FPKM value. The heatmaps are drawn using TBtools software. 13-LOX, 13S-lipoxygenase; HPL, hydroperoxide lyase; ADH, alcohol dehydrogenase; CHAT, acetyl CoA:(Z)-3-hexen-1-ol acetyltransferase; UGT, uridine-diphosphate-dependent glycosyltransferase; AOS, allene oxide synthase; AOC, allene oxide cyclase; OPR3, OPDA reductase 3; OPCL1, OPC-8:0 CoA ligase 1; ACX, acyl-CoA oxidase; MFP2, multifunctional protein 2; JMT, jasmonic acid carboxyl methyltransferase.

Identification of TFs in response to tea green leafhopper feeding

In addition to structural genes, TFs play a key role in plant stress responses. Here, a total of 2,377 TFs were identified and further classified into 53 families. The expression pattern based on the sum of the FPKM values of all members of the same family showed that the 53 families could be divided into four groups (Supplemental Fig. S1a). The TF family of the fourth group showed an overall up-regulated expression pattern, particularly the AP2, bHLH, C2H2, ERF, HD-ZIP, MYB, NAC, TCP and WRKY families. We then selected representative members from these families for separate expression pattern analysis. As shown in Supplemental Fig. S1b, these members were significantly induced by E. onukii feeding at the transcriptional level, suggesting that they could serve as key regulatory genes for enhancing tea plant defenses against E. onukii.

Variation characteristics of JA, JA-Ile, ABA and SA content

Defense responses mediated by plant hormone signal transduction pathways play a critical role in activating plant defense against herbivorous pests. Therefore, we further investigated the effects of feeding by the E. onukii on the changes in defense-mediated hormones JA, JA-Ile, SA, and ABA content in tea plant leaves. As shown in Fig. 6, the levels of JA, JA-Ile and ABA increased sharply after the feeding of the E. onukii, and although they decreased with time, the content of JA-Ile and ABA were still significantly higher than those of the control 48 h after ingestion. For SA, although its content continued to increase after feeding, compared with the control at the same period, it was only significantly higher than the control at 24 h after feeding, and significantly lower than the control at 6 h and 48 h.

Figure 6. 

Variation characteristics of main hormone (JA, JA-Ile, ABA and SA) content in tea plant leaves at different time points after E. onukii Matsuda feeding. Data are expressed as mean+standard deviation (SD) (n = 3). ** Indicates significant difference at the p < 0.01 level.

Identification of gene co-expression modules based on WGCNA

To comprehensively understand the gene expression patterns of tea plants responding to the feeding of E. onukii, the genes with average FPKM > 1 and coefficient of variation (cv) > 0.5 in 18 samples were used for WGCNA. Finally, 7,978 genes were divided into nine modules (Fig. 7a, b). The JA, JA-Ile, ABA and SA content in the corresponding samples were used as phenotypic data to analyze gene module-trait correlations. As shown in Fig. 7b, the black module, containing 82 genes, was most correlated with SA content (r = 0.88), whereas the green module, containing 666 genes, was highly associated with JA (r = 0.92), ABA (r = 0.97) and JA-Ile (r = 0.98) content at the same time. KEGG enrichment analysis of the genes in these two modules showed that the former was mainly concentrated in the 'biosynthesis of secondary metabolites' and 'phenylpropanoid biosynthesis' pathways, while the latter was mainly enriched in the 'plant hormone signal transduction', 'MAPK signaling pathway-plant' and 'α-linolenic acid metabolism' pathways (Supplemental Fig. S2), which was basically consistent with the enrichment pathways of the up-regulated DEGs. This indicates that these hormones may participate in the defense response of tea plants to E. onukii by regulating the expression of genes involved in the above-mentioned pathways.

Figure 7. 

WGCNA of genes responding to E. onukii Matsuda feeding in tea plant leaves. (a) Clustering dendrograms of genes. (b) Heatmap of correlation between module and trait. The vertical axis represents different modules, the left column of numbers represents the number of genes in the module and each set of data on the right represents the correlation coefficient. Gene co-expression networks of hub genes of (c) black module and (d) green module.

Co-expression networks were constructed based on the weighting values of genes in the two modules, and hub genes were screened based on their connectivity. In the black module, as shown in Fig. 7c, the amino acid permease (CsAAP, CSS0000586), glutathione S-transferase (CsGST, CSS0010310), mitogen-activated protein kinase (CsMPK3, CSS0008271), farnesoic acid carboxyl-O-methyltransferase (CsSABATH, CSS0048882), 7-deoxyloganetin glucosyltransferase (CsUGT, CSS0001530), β-glucosidase (CsGH1BG, CSS0014242) and two TFs CsWRKY (CSS0023344) and CsERF (CSS0007955) had a high degree of connectivity with other genes in the network. Therefore, they were identified as hub genes. The expression pattern showed that most of them were significantly up-regulated at 24 h after feeding (Supplemental Table S4), which was consistent with the significant accumulation of SA content at that time. For the green module, as shown in Fig. 7d, many key structural genes including those in the α-linolenic acid metabolism, terpene synthesis, and ABA synthesis and metabolism pathways, as well as some TFs (such as AP2, ERF, NAC, MYB, bHLH and WRKY) showed high connectivity in the network. We then focus on the top 50 connectivity members. Expression patterns showed that they were all remarkably up-regulated after feeding with E. onukii, and almost all had the largest up-regulation folds at 6 h after feeding (Supplemental Table S4). Among them, two aquaporin genes (CSS0042658 and CSS0044037), three lipoxygenase genes (CSS0003467, CSS0021769 and CSS0033612), one jasmonate-zim-domain protein (JAZ) gene (CSS0010510), one cationic peroxidase gene (CSS0005600), one methanol O-anthraniloyltransferase gene (CSS0032557) and one WRKY transcription factor gene (CSS0021684) were most dramatically up-regulated. Their functions should be studied further in the future.

Validation of the DEGs by qRT-PCR

Twenty differentially expressed genes including those from the α-linolenic acid metabolic pathway, the terpene synthesis pathway, and hub genes identified by WGCNA were used for qRT-PCR analysis to validate the RNA-seq data. As shown in Fig. 8, the relative expression of the selected genes was highly consistent with the RNA-seq results, indicating the reliability of the RNA-seq data. In addition, according to the expression pattern of these genes, it is suggested that they can also be used as key candidate genes for tea plants to respond to the feeding of E. onukii to conduct subsequent functional analysis to explore the mechanisms.

Figure 8. 

qRT-PCR validation of the expression levels of the genes in α-linolenic acid metabolism pathway, terpene synthesis pathway and the two co-expression networks in Fig. 7.

The defense response of tea plants to the feeding of E. onukii is a complex process, involving a large number of DEGs, DEPs and differentially accumulated metabolites^{[2,6,7,24, 27, 36, 37]}. It is worth noting that the accumulating evidence has shown that the crosstalk among various defense hormones to coordinate the expression of genes and their regulatory factors is the key to triggering the defense response of tea plants to phytophagous pests^[7,10]. In this study, based on transcriptomic and phytochemical analysis, we explored the response of tea plants at different time points after feeding with E. onukii, focusing on the changes of major defense hormones in this process and investigating their mediated gene regulatory networks. Similar to previous studies^[6,7,37], our results also showed that the E. onukii feeding elevated the expression of key genes involved in phenylpropanoid biosynthesis, flavonoid biosynthesis, α-linolenic acid metabolism and terpenoids synthesis pathways, especially in the early stage after feeding. In contrast to previous studies, we found that the degree of up-regulation of these genes at different times after feeding was significantly different. The expression levels of many genes did not reach significant levels compared with the control at 24 h after feeding, and were even lower than those in control. However, 48 h after feeding, their expression levels recovered and were significantly higher than those of the control. This phenomenon was also observed in previous studies. For example, Jin et al. showed that after E. onukii were fed for 0, 12, and 24 h, many differentially expressed genes in both resistant tea plant cultivar and sensitive cultivar showed a trend of increasing first and then decreasing^[6]. Another study showed that the number of DEGs at 24 h after feeding was higher than that at 12 h, but the number of up-regulated DEGs was lower than that of down-regulated^[37]. This is consistent with the results obtained 24 h after feeding in the present study. Therefore, we speculate that there may be a process of 'stress response – adaptation - defense response' in tea plants in response to E. onukii feeding. The early stage after feeding (such as 6 h, 12 h) is a 'stress response process' of tea plant to pest stimulation, then gradually adapt, and then start the 'defense response'. Similar expression patterns were also found at different time points after tea geometrid (a chewing insect) feeding^[38], indicating that this process is a general process by which tea plants respond to pests. However, in that study, many defense-related genes were still significantly up-regulated at 24 h after E. oblique feeding, while some genes (such as PAL and TPS) that were significantly up-regulated at 3 h after feeding were not significantly different from the control at 6 h after feeding^[38]. These results suggest that the timing difference in the emergence of this induction pattern may be related to the different feeding modes of the different pests. In addition, the accumulation of secondary metabolites in tea leaves is also regulated by circadian rhythm^[39]. Studies have found that the release of HIPVs and the expression of related synthetic genes are interactively regulated by the plant circadian clock and herbivorous pests^[40]. A recent study also showed that light could promote the accumulation of defense metabolites induced by the infestation of E. onukii in tea plants^[41]. Therefore, whether the circadian rhythm affects the defensive response of tea plants to the feeding of E. onukii needs to be further explored.

Previous studies have shown that plant hormone signals, especially JA signaling plays a key role in the response of tea plants to pests, including the feeding of tea green leafhoppers^[7,19,24]. Here, we found that almost all key genes in the JA synthesis pathway were significantly up-regulated, which persisted until 48 h after feeding. The JA content also increased dramatically. However, the accumulation of JA decreased over time after feeding. This finding was also observed in a previous study by Liao et al. Their results showed that although the JA content was significantly higher than that of the control at 96 h after feeding with E. onukii, its accumulation was lower than that at 48 h after feeding^[24]. In another study, the content of JA in the feeding leaves was only slightly increased after three weeks of continuous feeding^[7]. In addition to JA, JA-Ile and ABA showed similar trends. These results indicate that plants can induce the production of defense hormones to regulate the plant defense response rapidly after the recognition of insect pests, but do not induce their production constantly. In addition, the levels of these defense hormones declined inconsistently over time, suggesting a possible synergistic effect in terms of defense response activation. JA and JA-Ile are the major regulators of plant defense responses^[42], which will not be repeated here. It should be noted that ABA signaling is generally believed to play a key role in plant abiotic stress responses^[43]. However, Liao et al. found that the content of ABA increased significantly after the attack of E. onukii, while the tea geometrid attack did not^[24]. In this study, we also confirmed this phenomenon. In addition, the exposure of tea leaves to (E)-nerolidol increased the production of JA and ABA^[19], whereas (E)-nerolidol was one of the main volatiles induced by tea green leafhopper feeding^[11,19,44], indicating the special role of ABA signaling in tea plant response to E. onukii feeding. The E. onukii attack induced the formation of theaflavins from catechins, whereas exogenous ABA induced the accumulation of theaflavins^[24], suggesting that ABA signaling also plays a defensive role by inducing the formation of defensive metabolites. Furthermore, ABA is also involved in the synthesis of callose, which can not only resist the penetration of pathogens but also increase the obstacles to the attack of piercing-sucking pests^[45,46]. Therefore, the role of ABA signaling in the response of tea plants to the feeding of E. onukii deserves further investigation. Generally, pests with chewing mouthparts mainly trigger the JA signaling pathway, whereas pests with piercing-sucking mouthparts mainly activate SA signaling^[47,48]. In tea plants, it was found that feeding by tea geometrid and tea green leafhopper could significantly increase the content of JA and SA^[24], and the volatiles released from E. obliqua feeding on tea plants could induce significant accumulation of JA and SA in the adjacent uninjured tea plants^[21]. Another study found that feeding by tea green leafhopper only induced SA signals but not JA signals^[7]. In this study, SA was the only hormone that continued to accumulate after the attack of E. onukii, although it was significantly higher than that of the control only at 24 h after feeding. This phenomenon was also observed in the leaves of tea plants after feeding on tea geometrids. In that study, although the SA content in the leaves of tea geometrid at different time points after feeding was always lower than that in the corresponding control, the absolute SA content also showed an increasing trend with the increase in time after feeding^[23]. In this study, the trend for SA was opposite to that of the other three hormones. Antagonism between JA and SA has been widely confirmed in plants^[42], and the findings of this study indicate that there may also be antagonism between SA, ABA and JA-Ile. However, this antagonism was also synergistic. For example, after three weeks of continuous infestation by tea green leafhoppers, the content of JA was not significantly different from that of the control, whereas the content of SA was significantly increased^[7]. This indicates that activation of the SA signal may lag behind the JA signal, and may participate in the defense response of tea plants in the form of relays with JA signals. A recent study showed that the excitation timing and exciter concentration of JA and SA pathways decide the nature of JA-SA interactions in tea plant^[49]. In addition, the resistance of the cultivar itself also affects the accumulation of JA and SA after insect infestation^[50]. These results also suggest that there is a complex relationship between different hormones, which may be involved in plant defense responses through synergistic or antagonistic actions.

Given that the plant hormone signal transduction pathway is a common enrichment pathway of DEGs in tea plants in response to different treatments of tea green leafhoppers^[4,6,7], we further focused on the genes significantly related to the changes in the above different defense hormone contents. Not surprisingly, our WGCNA results showed that the four hormones were involved in two main co-expressed gene networks. In the SA-related co-expression network, two TFs (WRKY and ERF) and six functional genes (CsAAP, CsGST, CsMAPK, CsSABATH, CsUGT and CsGH1BG) were identified as hub genes. A previous study showed that the disease resistance elicitor laminarin could trigger early defense responses of tea plants to tea green leafhoppers by activating MAPK and WRKY expression and promoting SA accumulation^[51]. In addition, CsAAP and CsGST can participate in the transport of amino acids^[52] and flavonoids^[53], respectively, while CsSABATH, CsUGT and CsGH1BG can participate in the formation and metabolism of different volatiles^[54−56]. Therefore, we speculated that SA signaling might also mediate the defense response of tea plants to E. onukii by regulating the formation of related metabolites. JA, JA-Ile and ABA were in the same co-expression network due to their similar change trends at different time points after feeding with E. onukii. In addition to the upstream and downstream genes in the metabolic pathways of these hormones, many direct genes involved in HIPVs synthesis as well as different TFs were also included as hub genes in this network. Although the role of CsOPR3 in regulating tea plant defense against tea green leafhopper damage has been elucidated^[19], the roles of other genes in this process remain unclear. However, similar to CsOPR3, the expression of other hub genes were also significantly induced by E. onukii feeding, so it is clear that they also play a key role in regulating the response of tea plants to tea green leafhopper infestation, and their molecular functions in this process deserve further investigation.

In conclusion, this study investigated the defense response of tea plants to E. onukii from the perspective of changes in transcriptional dynamics at different time points after the tea plants were attacked. The simultaneous involvement of multiple defense hormones in the above defense response was clarified, and the specific roles of SA and ABA signaling in this process were further emphasized. The co-expression network construction and hub gene screening based on the dynamic changes in defense hormones provided important candidate genes and pathways to explore the resistance mechanisms of JA, JA-Ile, SA and ABA-mediated defense response of tea plants to E. onukii. In the future, we will focus on the verification of the interaction between hub genes to deeply elucidate the defense mechanism of tea plants against E. onukii regulated by defense hormone signals, thus laying a foundation for the breeding of resistant cultivars and the ecological control of E. onukii.