Metabolome and Transcriptome Analysis of the Response Mechanisms of Ginseng to Rust Root Symptoms

Background: Ginseng rusty root symptoms (GRS) is one of the primary diseases of ginseng. It leads to a severe decline in the quality of ginseng. Results: Compared with Healthy ginseng (HG), 949 metabolites and 9451 genes in diseased tissues were signicantly changed at the metabolic and transcription levels. The metabolic patterns of the diseased tissues changed signicantly, and organic acids, alkaloids, alcohols, and phenols may play a vital role in the response of ginseng to this disease. There were signicant differences in the expression of plant hormone signal transduction, phenylpropanoid biosynthesis, peroxidase pathway, and multiple genes in the plant-pathogen interaction pathway. Conclusion: The current study performed a comparative metabolome and transcriptome analysis of GRS and HG. Based on the ndings at the transcriptional and metabolic levels, the mechanism model of ginseng response to rusty root symptoms was established. Our results provide new insights into ginseng's response to rusty root symptoms, which will help reveal the potential molecular mechanisms of this disease in ginseng.

metabolism between them. To obtain the critical metabolites information between the two groups, we used the supervised partial least squares discriminant analysis (PLS-DA) method to separate the samples. PLS-DA score plot showed signi cant separation between the two groups of samples ( Fig. 1C-D). Also, to evaluate the reliability of the models, we performed the permutation test by iterating 200 times. The cross-validation results showed that neither the positive ion model nor the negative ion model was over-tted R 2 =0.79 and Q 2 =-0.90, R 2 =0. 74 and Q 2 =-0.82, respectively) ( Fig. 1E-F).
To search for differential metabolites in diseased ginseng tissues, we extracted VIP values from the variable importance plots of the PLS-DA models. Differential metabolites were screened based on three criterions: VIP > 1, P < 0.05 and Fold change > 2. According to the above criteria, 595 and 344 metabolites were screened out from the positive ion mode and negative ion mode, respectively. To better understand the chemical differences between diseased and healthy tissues, we plotted the heatmaps of 12 samples in positive and negative ion modes ( Fig.   2A-B). Most of the differential metabolites screened in positive ion mode had higher accumulation in diseased tissues (419 upregulated and 176 downregulated), and the same phenomenon was observed in negative ion mode (285 upregulated and 59 downregulated). By matching with the database, we nally identi ed 669 differential accumulated metabolites, of which 508 were up-regulated and 159 were downregulated ( Table 1). The differential metabolites in the two tissues were mapped to the Kyoto Encyclopedia of Gene and Genome (KEGG) database and were mainly enriched in multiple alkaloid biosynthesis pathways, such as "Biosynthesis of alkaloids derived from histidine and purine" and "Biosynthesis of alkaloids derived from ornithine, lysine and nicotinic acid". Besides, the extremely signi cant differences in "Stilbenoid, diarylheptanoid and gingerol biosynthesis", "Phenylpropanoid biosynthesis", "biosynthesis of plant hormones" and "Purine metabolism" are also noteworthy (Fig. 2C). The abundances of differential metabolites in 9 metabolic pathways with signi cantly enriched differential accumulated metabolites in KEGG metabolic pathways and more differential metabolites (compound number≥5) in the pathways were shown in Fig. 2D.

Quality control and de novo assembly of RNA-seq reads
We mainly analyzed mRNA obtained from the cDNA libraries. The sequencing of each ginseng sample produced over 100 million reads. After ltering adaptors and low-quality reads, we got high-quality clean reads, accounting for more than 96% of the raw reads (Table S3).

DEGs of diseased tissues
Here, we analyzed GRS and HG samples, the signi cance of the expression difference was analyzed by using the Cuffdiff software and the mRNAs were analyzed according to the padj < 0.05 standard. The differential expression of mRNAs in the samples between GRS and HG was shown by volcano plot and clustering map. The results showed that there were 4570 upregulated and 4881 downregulated mRNAs in GRS groups compared with HG groups (Fig. 3A-B).
To detect the gene function of mRNAs in ginseng rusty root, we further used Gene Ontology (GO) and KEGG pathway analysis to predict the function of RNAs. GO is the international standard classi cation system of gene function. The results of GO enrichment analysis of Differentially expressed genes (DEGs) were shown by a directed acyclic graph (DAG). All of the DEGs can be divided into three categories, including biological process (BP), cellular component (CC), and molecular function (MF). We have drawn DAGs in each of these three categories. The branch represents the relationship of inclusion, which de nes the scope from increasingly small from top to bottom. The top 10 results of GO enrichment analysis were selected as the master node of DAG and are shown together with related GO term by including the relationship and systematically GO were terms shown together, where the color depth represents the enrichment degree.
DAGs of BP, CC, and MF for GO enrichment analysis of the transcriptionally up-regulated mRNAs and the histogram of GO enrichment drawn after classi cation of the enriched GO terms are shown in Fig. 4A. On the basis of GO analysis of the transcriptionally upregulated mRNAs, the most signi cantly enriched BPs were "carbohydrate metabolic process", "steroid biosynthetic process", "oxidation-reduction process", "cell wall organization or biogenesis", and "response to oxidative stress", and the most noteworthy enriched CCs were "cell wall" and "betagalactosidase complex". The most signi cantly enriched MFs were "hydrolase activity", "oxidoreductase activity", and "peroxidase activity". In addition, DAGs and the histogram drawn according to transcriptionally downregulated mRNAs are shown in Fig. 4B. The results show that the most signi cantly enriched BPs were "oxidation-reduction process" and "terpenoid biosynthetic process", and the most noteworthy enriched CCs were "apoplast" and "cell wall". The most signi cantly enriched MFs were, "iron ion binding", "hydrolase activity", and "peroxidase activity". All of the enrichment results are shown in Tables S4 and Table S5.
For KEGG analysis of DEGs, the most signi cantly involved pathways were "metabolic pathways", "biosynthesis of secondary metabolites", and "phenylpropanoid biosynthesis" (Fig. 4C). All pathways of DEGs based on KEGG enrichment analysis are given in Table S6.
Response of ginseng to disease GO and KEGG enrichment analysis showed the response of GRS. Compared with HG, numerous DEGs were involved in plant hormone signal transduction, lignin synthesis, plant-pathogen interactions, and oxidative stress.
There were 9 differentially metabolites and 124 DEGs involved in the pathway "phenylpropanoid biosynthesis" (Fig. 2D). The detailed locations of differentially metabolites and DEGs in the pathway are shown in Fig. 5. In GRS, enzymes closely related to lignin synthesis, such as phenylalanine ammonia-lyase (PAL), O-methyltransferase1 (OMT1), 4-coumarate: coenzyme A ligase1 (4CL1), and cinnamoyl coenzyme A reductase (CCR) were detected with a high expression level.
We analyzed the pathways "peroxisome", which is closely related to oxidative stress in plants. As shown in Fig. S3 and Fig. S4, we found that 59 DEGs were involved in "peroxisome". Among them, several genes of the PEX family, such as PEX14, PEX16, and PEX19, were downregulated in the GRS group. PXMP2 and MPV17 were involved in reactive oxygen species (ROS) metabolism and showed obvious downregulated expression. Additionally, catalase (CAT) and superoxide dismutase (SOD) in the antioxidant system were upregulated.

Correlation analysis of differential accumulation metabolites and DEGs
We found that many differentially expressed genes and differentially accumulated metabolites were enriched in the same KEGG metabolic pathway in disorganized tissues. These pathways include "Phenylpropanoid biosynthesis", "Phenylalanine metabolism", "Stilbenoid, diarylheptanoid and gingerol biosynthesis", and so on (Fig. 6A). Besides, considering the role of the phenanthrene biosynthesis pathway in disease resistance and stress response in plants, we constructed a correlation network graph of differentially expressed genes and differential metabolites based on Pearson correlation coe cient (PCC). As shown in Fig. 6B, a total of 36 DEGs and 8 differential metabolites in the network were highly positively correlated (PCC > 0.8). It is speculated that these genes and metabolites may play a key role in ginseng response to stress.

Real-time quantitative polymerase chain reaction (qPCR) validation of mRNAs expression
To verify the accuracy of RNA-seq results and provide the basis for further research, some genes from the four pathways of plant resistance were selected for qRT-PCR analysis and the expression level of these genes was demonstrated by qRT-PCR and transcriptome sequencing.

Discussion
GRS is one of the most common diseases at present, which has a great impact on the quality of ginseng and seriously threatens development of the ginseng industry. Previous studies have suggested that GRS is associated with the accumulation of metal ions and the metabolism of phenolic substances in the root epidermis [2]. Scholars have argued that the red substance in the epidermis of ginseng root may have been formed by a complexation reaction. Additionally, GRS involves microbial infestation [19]. The occurrence of GRS is a complex mechanism, which is not clear at present.
In this study, metabolites and total RNAs were extracted from the same variety of 5-year-old ginseng planted in the same ginseng farm (higher rusty root index) were prepared from healthy and diseased ginseng tissues. The purpose was to ensure that ginseng samples experienced the same eld management and climate change in the growth process [20].
Metabolome and transcriptome analysis methods have provided a basis for the study of many plant diseases, and a large number of previous studies have also found the response of metabolic and transcription levels to abiotic stress and biological infection in plants [21][22][23][24][25]. Therefore, we conducted UHPLC-MS/MS and RNA-Seq analysis to the diseased ginseng and HG. Our results showed that in GRS, there were 939 metabolites differentially accumulated and 9451 genes signi cantly differentially expressed.
Metabolite accumulation differences in diseased tissues compared with HG were analyzed by metabolomics. The results showed that the metabolite expression patterns of the two ginseng tissues were signi cantly different. Among the metabolites up-regulated to accumulate in diseased tissues, including organic acids, alkaloids, alcohols, and phenols, they have been closely associated with plant responses to biotic and abiotic stresses [26]. Our study revealed that many organic acids produced signi cant changes in diseased tissues of ginseng (Table 1).
Organic acids are essential intermediates in plant energy balance and ow and are involved in plant response to abiotic stresses. For example, α-linolenic acid is a jasmonic acid precursor compound, which plays a vital role in plant resistance by forming jasmonic acid and generating methyl jasmonate through enzymatic reactions [27,28]. However, considering that most organic acids act as intermediate metabolites, it is di cult to determine their speci c stress response role [29]. Alkaloids are structurally diverse and have many biological activities that affect plant growth. The accumulation of alkaloids is a vital defense strategy for plants in response to biotic stresses and is a crucial class of compounds involved in plant chemical defense [30]. In diseased tissues, we also found an upregulation of the accumulation of multiple alkaloid metabolites. The KEGG enrichment results of differential metabolites also enriched multiple alkaloid biosynthetic pathways (Fig. 2C). This result suggests that alkaloid metabolites are essential for the stress response of ginseng. Chemical defense involving alkaloids may occur in diseased tissues of ginseng against pathogenic bacteria and abiotic stresses [31]. Besides, we also noted that multiple alcohol metabolites showed higher accumulation in diseased tissues. It was shown that the accumulation of alcohols usually plays a role in the ght against oxidative stress. They can act as intermediates of redox reactions with good scavenging capacity for free radicals and superoxide [32,33]. It can be speculated that ginseng accumulates more alcohols in diseased tissues to cope with the oxidative stress in itself to improve its resistance. Besides, we found that various phenolics were also induced in diseased tissues in response to external stress. The phenolic content is often considered closely related to the total antioxidant capacity of plant tissues [34]. Studies have shown that avonoids are also one of the main chemical defense substances in plants and can reduce various forms of reactive oxygen species in plant cells, which are usually produced in response to external environmental stresses such as light intensity, temperature, and soil moisture [35]. Besides, it is also an active substance involved in plant-microbe interactions [36]. In our results, a signi cant upregulation of a large number of avonoids was found. Signi cant enrichment of the avonoid biosynthetic pathway was also found by KEGG enrichment analysis (Fig. 2C). Overall, these compounds accumulate in diseased tissues of ginseng and are jointly involved in the oxidative stress response and chemical defense of ginseng.
Terpenoids, especially saponins, are the main active substances in ginseng. We found signi cant accumulation or down-regulation of several saponin substances in diseased tissues, such as Ginsenoside Ro, 20(S)-Ginsenoside Rh2, and Pseudoginsenoside F11. In a previous study, ginsenosides had inhibitory effects on pathogens [37]. It is speculated that GRS may impact ginseng terpene metabolism, and further studies are needed to determine whether saponins play a role in disease response.
Lipids are essential membrane components, and changes in their composition may help plants maintain cellular compartmentalization and membrane integrity [38]. Besides, free radicals are metabolized and produced when cells are under stress or in the process of injury. Their scavenging capacity is still reduced, leading to an imbalance of reactive oxygen species and membrane lipid peroxidation. Therefore, reducing lipid peroxidation has an essential role in plant resilience. In diseased tissues of ginseng, there were more lipid metabolites with signi cantly altered content. For example, phosphatidylcholine (PE) and phosphatidylethanolamine (PC) are lipids that are present in large amounts in cell membranes. An increase in the PC/PE ratio has been reported to be a marker of the membrane bilayer's structural stability and prevent membrane degradation [39]. Interestingly, we detected a signi cant increase in several PC levels and a substantial decrease in PE levels in diseased tissues of ginseng. This result may imply that diseased tissues may have a stronger tendency to maintain bio lm integrity and stability than healthy tissues.
GO enrichment analysis helps to reveal the different gene functions of biological processes, molecular functions, and cellular components. GO annotations have been shown to be good predictors of gene function and trends [40,41]. Our data showed that the genes related to "cell wall organization" (GO:0071554) and "biogenesis response to oxidative stress" (GO:0006979) were abundant on BP. The CC term had the enrichment of "cell wall" (GO:0005618). Additionally, "peroxidase activity" (GO:0004601) was enriched in the MF term. This suggests that there may be a strong oxidation-reduction reaction in the diseased tissues of ginseng, as well as the cell wall organization and oxidative stress.
We found that the biosynthesis of secondary metabolites was signi cantly enriched at both the metabolic and transcriptional levels, which fully demonstrated that the secondary metabolism of diseased tissue was altered considerably in comparison with that of HG tissues. Among them, metabolic changes of alkaloids and avonoids in ginseng root tissues may be important chemical defense patterns.
Based on the results of metabolites and DEGs (GO and KEGG) enrichment analysis, we found that plant hormones, cell wall formation, oxidative stress and plant-pathogen interaction were closely related to plant response to disease. Plants produce a variety of hormones during growth. Among them, ABA, SA, and JA play an essential role in mediating plant defense responses to pathogens and abiotic stress [42,43].
Here, three plant hormones and hormone-related signaling molecules were detected by UHPLC-MS/MS, and all of them were signi cantly upregulated in diseased tissues (ABA, SA, and MeJA) (Fig S1A). At the transcriptional level, the ABA receptor PYR/PYL/RCAR, negative regulator PP2C, and positive regulator SnRK2 together comprise the regulatory system [44]. In GRS, the genes that regulate ABA signal transduction are differentially expressed and ultimately ABF transcription factor upregulation (Fig. S1B). Under JA-stimulated conditions, JA-Ile, synthesized by JAR1, binds to the receptor COI1, leading to degradation of the repressor protein JAZ, which upregulates the expression level of JA target genes (Fig. S1C) [45,46]. NPR1 is one of the key regulatory elements of SA-dependent PR gene activation, which interacts with GTA factors to regulate PR gene expression [47,48]. We found differential expression of genes of the byproducts of NPR1 (BOP2 and NPR3) and multiple genes of TGA transcription factors that culminated in upregulation of PR-1 gene expression (Fig. S1D). Obviously, ABA, SA, and MeJA were fully involved in the stress response of diseased tissues. These changes in crucial hormone levels may be an early response to stress. The affect metabolic processes and eventually lead to changes in growth patterns that are adapted to withstand environmental pressures.
Lignin is a highly branched polymer of phenylpropanoid compounds, which is the main component of the plant cell wall [49]. Because lignin is a kind of non-degradable mechanical barrier for most microorganisms, it plays an important role in plant defense. This is considered one of a series of mechanisms for plants to prevent microbial invasion. Lignin biosynthesis is a complex genetic network, in which many enzymes are involved [50,51]. Among these enzymes, several key enzymes, including PAL, 4CL, CCR, OMT, and CAD were upregulated (Fig. 5) [52,53]. Based on these analysis, we believe that lignin may accumulate in GRS, which usually occurs when plants are attacked by pathogens. In addition, we detected 9 compounds involved in phenylpropanoid metabolic pathways through UHPLC-MS/MS (Fig. 2D). The relative levels of phenolic compounds, such as eugenol, coniferin and chlorogenate, were far higher in the disease tissues than in HG.

Exposure of plants to various biotic and abiotic stress conditions triggers rapid changes in reactive oxygen species (ROS) production and
clearance [54]. ROS play an important role in signaling pathways that regulate in ammatory and defense responses in plants but its accumulation is often harmful to cells. Peroxisomes are subcellular organelles with a basic oxidative type of metabolism and may be the primary site of intracellular ROS production [55,56]. This organelle is a small, usually spherical body, encased by a lipidic double membrane. Enrichment analysis with KEGG showed that both PXMP2 and MPV17, involved in ROS metabolism, were transcriptionally downregulated ( Figure S2). PXMP2 is a perovskite membrane protein with the channel-forming activity that allows free diffusion of compounds, such as H 2 O 2 [57]. MPV17 is a protein involved in the production of ROS [56]. Several PEX family protein genes were differentially expressed, including PEX13 and PEX14, involving matrix protein import, and PEX16 and PEX19, involving membrane protein import [58,59]. Moreover, although peroxisomes regulate the production and release of ROS and CAT and SOD may also play an antioxidant role, oxidative stress may still be one of the causes of death of epidermal cells in diseased tissues [60].
Based on the available ndings, the plant immune system depends on the cell's innate immunity and systemic signals from the site of infection. Plants have two layers of defense against pathogens, PTI and ETI [61]. Here, compared to HG, there were 96 DEGs in GRS associated with the plant-pathogen interaction pathway. RBOH mediates ROS generation, which is strongly associated with both PTI and ETI [62,63]. RBOHD and RBOHF are considered to be critical components of plant defense [64]. FLS2 is a pattern recognition receptor (PRR) that can trigger innate immunity in plants and is also involved in endocytosis [65]. WRKY33 plays a vital role as a transcription factor in plant resistance to pathogens [66]. The downregulation of the above protein genes, as well as NOA1 and MEKK1, may indicate that part of PTI is repressed in GRS. Interestingly, ETI seems to be activated in GRS. RPS2 is a resistance (R) protein in plants, while RIN4 is an important negative regulator of natural plant immunity and has a regulatory effect on RPS2 [67]. HSP90.1 is required for full RPS2 resistance and is rapidly induced in plants in the face of biological stress [68]. SGT1 positively regulates disease resistance produced by many R proteins [69].
In our analysis, RIN4 was downregulated and PRS2, SGT1, and HSP90.1 were upregulated. Furthermore, NHO1, an essential universal resistance gene in plants, was upregulated in expression [70]. Overall, the immune system of diseased tissues is activated, which will indirectly lead to the hypersensitive response and cell wall reinforcement. Simultaneously, the discovery of many DEGs in the "plant-pathogen interaction" pathway also suggested that GRS is related to pathogen infection.
To further understand the critical aspects of ginseng response to GRS, we combined metabolomic data with DEGs for analysis. Many differential metabolites and DEGs were annotated in the same metabolic pathway (Fig. 7A). Among them, both parts of the data were signi cantly enriched in the phenylpropane metabolic pathway. Considering that the phenylpropanoid biosynthesis pathway is closely related to phenolic production, avonoid metabolism, and lignin formation, it is speculated that it may be a critical metabolic pathway in ginseng in response to GRS [71]. Therefore, we constructed a correlation network of differential metabolites and DEGs in this metabolic pathway based on the results of correlation analysis (Fig. 7B), which may play a direct or indirect regulatory role in altering the corresponding metabolites.

Conclusions
The current study performed a comparative metabolome and transcriptome analysis of GRS and HG. Metabolic patterns in diseased tissue have changed dramatically, especially secondary metabolism. Analysis of the pathways associated with the plant resistance response revealed that hormone signaling pathways, lignin biosynthesis, ROS regulation by peroxisomes, changes in plant-pathogen response pathways, and the secondary metabolites produced in many processes might play essential roles in GRS. Based on previous studies and the results of this metabolomics and transcriptome analysis, we propose a model to explain ginseng's response to rust root symptoms (Fig. 8).
Pathogen infestation and abiotic stress induce oxidative stress, which activates the PTI system and stimulates the synthesis of SA, MeJA, and ABA in ginseng tissues. SA and MeJA induce the accumulation of phenolic, alkaloid, and alcohol metabolites. SA and MeJA also induce lipid and fatty acid metabolism and enhance membrane stability against stress. The production of antioxidant enzymes such as CAT and SOD and the peroxisome regulation together affect ROS accumulation in tissues. Meanwhile, differential expression of genes in PTI and ETI in ginseng epidermal tissues may promote PR proteins' production and hypersensitive responses against external stresses. The activation of the phenylpropanoid biosynthesis pathway enhances the production of phenolic metabolites on the one hand. It may increase the accumulation of lignin and promote the cell wall's reinforcement on the other hand. The reinforcement of the cell wall and the enhancement of membrane stability together constitute a physical defense in ginseng tissues. Our results provide new insights into ginseng's response to rusty root symptoms, which can reveal the potential molecular mechanisms of this disease in ginseng.

Plant material and tissue collection
All ginseng samples (5 years old) were collected from the same ginseng farm in Hunchun city, Jilin Province, China (42.86′N and 130.37′E). In the previous survey, this was a ginseng farm with a high rusty root index [72].
The diseased tissues of GRS main root and the healthy tissues of HG were cut off as samples and frozen in liquid nitrogen immediately and stored at -80°C before use. The independent biological replicates were prepared and each replicate included root materials from three or more of ginseng plants.

Metabolome analysis
The rst is the extraction of metabolites. Tissues (100 mg) were individually grounded with liquid nitrogen, and the homogenate was resuspended with prechilled 80% methanol and 0.1% formic acid by a well vortex. The samples were incubated on ice for 5 min and then were centrifuged at 15,000 rpm, 4°C for 5 min. Some of the supernatants were diluted to nal concentration containing 53% methanol by UHPLC-MS/MS grade water. The samples were subsequently transferred to a fresh Eppendorf tube and then were centrifuged at 15000 g, 4°C for 10 min.
Then, the UHPLC-MS/MS analysis was performed using a Vanquish UHPLC system (Thermo Fisher, Germany) coupled with an Orbitrap Q Exactive TM HF-X mass spectrometer (Thermo Fisher Germany). Samples were injected onto a Hypesil Gold column (100×2. The raw data les generated by UHPLC-MS/MS were processed using the Compound Discoverer 3.1 (CD3.1, Thermo Fisher) to perform peak alignment, peak picking, and quantitation for each metabolite. The main parameters were set as follows: retention time tolerance, 0.2 minutes; actual mass tolerance, 5ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, 100 000. After that, peak intensities were normalized to the total spectral intensity. The normalized data were used to predict the molecular formula based on additive ions, molecular ion peaks and fragment ions. And then peaks were matched with the mzCloud, mzVault and MassList database to obtain the accurate qualitative and relative quantitative results.
PCA and PLS-DA was performed at SIMCA-P 14.1 (Umetrics, Umea, Sweden). The metabolites with VIP > 1 and P-value < 0.05 and fold change > 2 were considered to be differential metabolites. The KEGG enrichment analysis on differential metabolites using KOBAS 2.0 software [73].
RNA extraction, library preparation, clustering and sequencing A total of six complementary DNA (cDNA) libraries, three GRS and three HG, were constructed. TRIzol reagent (Invitrogen, Carlsbad, CA) was used to isolate the total RNA of each sample. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).
In this part, a total of six complementary DNA (cDNA) libraries, three GRS and three HG were constructed. A total amount of 3 µg RNA per sample was used as input material for RNA sample preparation. Sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA), following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Brie y, mRNA was puri ed from the total RNA using poly-T oligo-attached magnetic beads. Fragmentation was conducted using divalent cations under elevated temperatures in NEBNext First Strand Synthesis Reaction Buffer (5X). The rst strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). The second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments of preferentially 150-200 bp in length, the library fragments were puri ed with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95°C before PCR. PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. PCR products were puri ed (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
After cluster generation, RNA library preparations were sequenced on an Illumina Hiseq PE150 platform, and 150 bp paired-end reads were generated. The raw RNA-seq data are freely available in the NCBI database under accession no. PRJNA684799.

Transcriptome data analysis
Raw data (raw reads) of the fastq format were rst processed through in-house Perl scripts. Then, a certain range of length was chosen from clean reads to do all of the downstream analyses.
The reference genome and gene model annotation les were downloaded from the genome website directly [74]. The index of the reference genome was built using STAR, and paired-end clean reads were aligned to the reference genome using STAR (v2.5.1b). STAR used the method of the maximal mappable pre x (MMP), which can generate a precise mapping result for junction reads. HTSeq v0.6.0 was used to count the read numbers mapped to each gene. The expected number of Fragments Per Kilobase of exon model per Million mapped fragments (FPKM) of each gene was calculated based on the length of the gene and read count mapped to this gene.
Differential expression analysis of two groups was performed using the DESeq2 R package (1.10.1). The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate.
GO (http://www.geneontology.org/) enrichment analysis of DEGs was implemented by the GOseq based Wallenius non-central hypergeometric distribution, in which gene length bias was corrected [75]. Besides, we performed KEGG enrichment analysis on DEGs using KOBAS 2.0 software [73].
Quantitative real-time qPCR The quantitative PCR reaction was performed on the Real-Time PCR System (Lightcycler 96, Roche, Switzerland) using 2 × RealStar Green Fast Mixture (GenStar, China). For each reaction, 0.5 μl of the forward and reverse primers and 2 μl of cDNA template were added. All of the primers used in this study are listed in Table S7. The relative gene expression level was calculated according to the 2 -ΔΔCt method [76].
Combined analysis of DEGs and differentially expressed metabolites The obtained differential metabolites and DEGs were mapped to the KEGG pathway for integration. Besides, we performed correlation analysis on transcriptome and metabolomics data to explore correlations between differentially expressed genes and metabolites. The PCC of DEGs and differentially metabolites was calculated using the COR function in R. The genes and metabolites of PCC > 0.8 in the KEGG pathway were used to establish the related network, which was visualized by Cytoscape software.

Statistical analysis
Statistical analysis was performed with SPSS 20.0 software. All of the data were expressed as the mean ± SEM. P < 0.05 was considered statistically signi cant. RT-PCR results were analyzed using a Student's t-test. Signi cance was established at P < 0.05. Table S1. Results of UHPLC-MS/MS analysis in the positive ion mode. Table S2. Results of UHPLC-MS/MS analysis in the negative ion mode. Table S3. Summary of RNA-seq data and reads mapping for mRNA.     Figure S1. Plant hormones content and heatmaps of DEGs. Figure S2. The plant hormone signal transduction pathway. Figure S3. Heatmap of DEGs related to peroxisome synthesis. Figure S4. The peroxisome pathway. Figure S5. The plant-pathogen interaction pathway. Figure S6. Heatmap of DEGs related to plant-pathogen interaction.