Transcriptome and metabolome reveal the accumulation of secondary metabolites in different varieties of Cinnamomum longepaniculatum

Background Cinnamomum longepaniculatum (Gamble) N. Chao ex H. W. Li, whose leaves produce essential oils, is a traditional Chinese medicine and economically important tree species. In our study, two C. longepaniculatum varieties that have significantly different essential oil contents and leaf phenotypes were selected as the materials to investigate secondary metabolism. Result The essential oil content and leaf phenotypes were different between the two varieties. When the results of both transcriptome and metabolomic analyses were combined, it was found that the differences were related to phenylalanine metabolic pathways, particularly the metabolism of flavonoids and terpenoids. The transcriptome results based on KEGG pathway enrichment analysis showed that pathways involving phenylpropanoids, tryptophan biosynthesis and terpenoids significantly differed between the two varieties; 11 DEGs (2 upregulated and 9 downregulated) were associated with the biosynthesis of other secondary metabolites, and 12 DEGs (2 upregulated and 10 downregulated) were related to the metabolism of terpenoids and polyketides. Through further analysis of the leaves, we detected 196 metabolites in C. longepaniculatum. The abundance of 49 (26 downregulated and 23 upregulated) metabolites differed between the two varieties, which is likely related to the differences in the accumulation of these metabolites. We identified 12 flavonoids, 8 terpenoids and 8 alkaloids and identified 4 kinds of PMFs from the leaves of C. longepaniculatum. Conclusions The combined results of transcriptome and metabolomic analyses revealed a strong correlation between metabolite contents and gene expression. We speculate that light leads to differences in the secondary metabolism and phenotypes of leaves of different varieties of C. longepaniculatum. This research provides data for secondary metabolite studies and lays a solid foundation for breeding ideal C. longepaniculatum plants. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03637-2.

Terpenoids (isoprenoids) are a class of important chemicals produced by plants [12,13]. Many terpenoids have a long history of being used as flavour-enhancing compounds, pharmaceuticals, insecticides, and industrial compounds [14]. Terpenoids in plants can be divided into monoterpenes (C10), sesquiterpenes (C15), diterpenes (C20), triterpenes (C30), tetraterpenes (C40) and polyterpenes [15]. Plants use two independent pathways to produce terpenoids: the cytosolic mevalonic acid (MVA) pathway, and the plastidial methylerythritol phosphate (MEP) pathway [16][17][18]. During the last few decades, there has been increasing interest in the terpenes of Cinnamomum plants, mainly focusing on analysis of the main compounds of essential oil extracted from leaves [12,14]. Many studies have attempted to evaluate the antibacterial activity of terpenoids that constitute essential oils [19,20]. In recent years, many studies have focused on analysing the transcriptomes related to specific terpenoids in leaf metabolic pathways and their regulatory mechanisms [21] and the differential accumulation of terpenoids in different Cinnamomum chemotypes [22]. However, combined metabolome and transcriptome profiling of terpenoids in different varieties of C. longepaniculatum has not been reported.
Cinnamomum longepaniculatum (Gamble) N. Chao ex H. W. Li [23] is an evergreen tree species [24]. C. longepaniculatum is extensively cultivated in southwestern China, especially in the Yibin region (Sichuan, China). C. longepaniculatum leaves and twigs are harvested for essential oil extraction. Research on essential oils has mainly focused on their composition [25], especially their antibacterial activities [26] and endophytic fungi [27]. Terpenoids (monoterpenes and sesquiterpenes) produced during plant secondary metabolism are the major components of essential oils [19]. However, we know little about other secondary metabolites present in the leaves of C. longepaniculatum.
In this study, two different varieties of C. longepaniculatum, namely, CLH and CLL, were selected as the materials to study the secondary metabolite differences. Previous research has shown that the production of essential oils significantly differs among varieties, so we used transcriptomic and metabolomic data from functional leaves from different varieties to investigate their secondary metabolism. This research not only provides data for C. longepaniculatum secondary metabolite studies but also lays a solid foundation for breeding ideal C. longepaniculatum plant types.

Leaf essential oil content and phenotypic changes of different varieties
We sampled fresh leaves from different varieties for three years and extracted essential oils in the same way. Under the same climatic conditions, the essential oil contents of the different varieties of C. longepaniculatum were different. The results showed that the essential oil content of CLH was significantly higher than that of CLL (Table 1). Many differences in the leaves of the two varieties were detected. The leaf area of CLH was greater than that of CLL, but the difference was not significant (Fig. 1A). We compared cross-sections of leaves from the different varieties by paraffin sectioning. The cross-section of the midrib revealed a biconvex shape, and the average crosssection area for CLH was 82.94 μm 2 , while that for CLL was 72.82 μm 2 . We also measured the leaf thickness, and the average for CLH (226.57 μm) was significantly greater than that for CLL (189.10 μm). Many ground parenchyma (gp) cells were observed on the veins of the adaxial sides of the leaves. The vascular system was composed of two vascular bundles, which were arranged in the open arc shape, and the middle was very large. Secretory idioblasts (Si) were present between the midrib and lateral vein of the CLH leaves, but these were not present in the CLL leaves (Fig. 1B, C). The mesophyll comprised one layer of the epidermis. Small collateral vascular bundles were  22:243 bordered by the endoderm and penetrate through the spongy parenchyma (Fig. 1D, E).

Analysis of transcriptomic data
To explore the molecular mechanisms underlying secondary metabolism, we compared transcriptome between the different varieties. A robust data set was generated after data processing: 48.63 million and 46.68 million high-quality reads were obtained from CLH and CLL, respectively; moreover, 7.34 and 7.05 Gb were obtained, and the GC contents were 45.59 and 45.74%, respectively. The sequencing statistics of the 6 RNA generated libraries are listed in Table 2. BUSCO analysis showed that the completeness of C. longepaniculatum leaf transcriptome was 95.2%, indicating the excellent continuity of the assembly (Table 3). To compare the gene expression in CLH with that in CLL, Venn diagram analysis showed that 34,503 genes were present in both comparison groups (Fig. S1).
The KEGG analysis showed that the phenylpropanoid, tryptophan biosynthesis and terpenoid pathways were significantly different between the two varieties (Table  S3). We were more interested in the pathways related to secondary metabolism: biosynthesis of other secondary metabolites (2 upregulated and 9 downregulated DEGs) and metabolism of terpenoids and polyketides (2 upregulated and 10 downregulated DEGs). In addition, there were 3 KEGG pathways that were related to light: photosynthesis (map00195, 35 DEGs), photosynthesisantenna proteins (map00196, 18 DEGs), and oxidative phosphorylation (map00190, 29 DEGs) (Table S2).

Metabolome profiling of different varieties
Metabolites were extracted from leaf samples (three replicates) for each experimental group and were analysed by GC-MS. Partial least squares discrimination analysis (PLS-DA) resulted in the division of the total variation into two major components (Component 1 and Component 2), which contributed 41.3 and 23.8% of the variation, respectively (Fig. 3A). In total, we detected 196 metabolites from C. longepaniculatum (Table S4): 147 were matched with KEGG database information, and 92 were characterized in the library. The analysis of the function of the metabolites revealed nine main aspects: amino acid metabolism (34), biosynthesis of other secondary metabolites (24), carbohydrate metabolism (24), energy metabolism (8), lipid metabolism (9), metabolism of cofactors and vitamins (17), metabolism of other amino acids (11), metabolism of terpenoids and polyketides (6), and nucleotide metabolism (8) (Fig. S2).

Correlation analysis between the transcriptomic and Metabolomic data
Correlation analysis between metabolites and transcripts of different varieties was performed to gain insight into the regulatory network of secondary metabolism of C. longepaniculatum. We conducted a joint KEGG pathway enrichment analysis, and identical pathways were enriched in the transcriptome and metabolome. A histogram was construct to visualize the enrichment degree of pathways with both DEGs and DAMs (Fig. 4A). The results showed that three KEGG pathways were related to phenylalanine metabolism: phenylalanine metabolism (map00360); phenylalanine, tyrosine and tryptophan biosynthesis (map00400); and phenylpropanoid biosynthesis (map00940). According to the results, 7 DEGs involved in the amino acid pathway (3) and the secondary metabolism (4) pathway were related to phenylalanine metabolism. After analysis of the secondary metabolites, 9 metabolites were found to be associated with phenylpropanoid biosynthesis. We used two-way orthogonal partial least squares (O2PLS) to evaluate the intrinsic correlation between the transcriptomic and metabolomic data, calculated the score of each sample and obtained a joint score. Then, the load value of each gene and metabolite was calculated to obtain a load map (Fig. 4B). A joint score plot was constructed to visualize the relationship between the two data matrices, and metabolites/genes with high loading values were considered necessary for the similarity of the two data sets. Finally, the absolute values of the load value of the top 15 DAMs/DEGs were selected to construct a histogram (Fig. 4C). The correlation analysis showed that transcripts and metabolites were related to phenylalanine metabolic pathways.

Confirmation of DEGs via qRT-PCR
To validate the reliability and stability of the RNA sequencing (RNA-seq) data of the DEGs, 8 DEGs were selected for qRT-PCR validation. Among them, four were upregulated in CLH (Fig. 5). The results suggested that the transcriptomic data were reliable.

Discussion
Between the two varieties, 49 (26 downregulated and 23 upregulated) DAMs, 12 flavonoids, 8 terpenoids and 8 alkaloids were identified. Combined transcriptome and metabolome analyses revealed a strong correlation between metabolite content and gene expression (Fig. 6). We identified 12 flavonoids in C. longepaniculatum leaves (Table 4). Flavonoids are low-molecular-weight secondary metabolites produced by plants [28]; are present in all plant parts; and show a variety of functions, including antioxidant activity [29,30], UV protection [31], pathogen protection and signalling [32,33]. Constituting a large group of plant secondary metabolites, flavonoids are derived from the phenylpropanoid pathway [34]. Flavonoids, which include, flavones, flavonols, flavan-3-ols, flavanones, isoflavanones, isoflavans, and pterocarpans, have diverse structures. Previous studies suggest that apigenin [35], chrysin [36], neohesperidin [37], epicatechin [38], catechin [39], and kaempferol [40,41] are bioactive substances that have antioxidant, antibacterial, antifungal, and antitumor activities and antiinflammatory properties. Therefore, we believe that the Fig. 4 Correlation analysis between metabolites and transcripts. A Histogram of the KEGG pathway enriched by the differentially expressed genes and metabolites between CLH and CLL. The ordinate represents the pathway name, and the abscissa represents the pathway to impact. B Assess the intrinsic correlation between the metabolites and transcripts by the two-way orthogonal partial least squares (O2PLS) method. Solid blue squares represent genes and solid green dots represent metabolites. The abscissa and ordinate represent the combined loading value, each gene/metabolite has a relative coordinate point in pq1 and pq2, p represents the loading value of the gene, and q represents the loading value of the metabolite. C The top15 differential expresses gene and differential metabolite loading histogram. Blue bars represent genes and green bars represent metabolites. The abscissa represents the combined loading value pq1, and the ordinate represents the differential metabolite/differential gene flavonoids in C. longepaniculatum leaves are related to the plant stress response. It may be that endophytic fungi [27] cause the accumulation of flavonoids in the leaves of C. longepaniculatum, but we believe that the accumulation of flavonoids is caused by light. There is very important evidence in that 82 DEGs were mapped to 3 KEGG pathways that were related to light. We also found that the KEGG oxidative phosphorylation pathway was enriched in DAMs. Therefore, we speculate that different varieties of C. longepaniculatum have different responses to light, which results in differences in flavonoid contents. Nonetheless, this hypothesis needs further experimental confirmation. PMFs constitute a unique class of flavonoids that are considered the biologically active constituents in citrus fruits [42]. We identified 4 kinds of PMFs from the leaves of C. longepaniculatum; this is the first report on the isolation of PMFs from C. longepaniculatum. PMFs can be used for anti-inflammatory, antiatherogenic, antiangiogenic and anticancer activities [43,44], and different PMF structures lead to differences in biological activity [45]. By comparison, the PMFs identified in the present study are different from those in citrus, and the biological activities of PMFs from C. longepaniculatum need to be studied. Some studies have pointed out that citrus PMFs can affect the modulation of gut microbiota [10], and others have pointed out that microbial community composition and diversity can affect secondary metabolite content in leaves of Cinnamomum camphora [46]. We speculate that there is a complex relationship between PMFs and endophytic fungal communities of C. longepaniculatum.
A total of 7 classes and 8 kinds of terpenoids were identified in our research (Table 4). These components are largely consistent with the findings of previous studies [13,46,47]. Owing to their pleasant fragrance, unique biological activities and favourable physicochemical properties, terpineol [48,49], eucalyptol [50], limonene [51], and squalene [52] have strong antioxidant properties and are widely applied in foods, pharmaceuticals, cosmetics, biomaterials, and biofuels [53,54]. Among the terpenoids we identified, the contents in CLL were higher than those in CLH, except for squalene and beta-ionone; squalene synthesis occurs downstream of terpenoid synthesis. This difference between varieties may be because more substances may be converted to other terpenes, such as terpineol, in CLL and because less squalene is formed downstream. We analysed 196 terpenoids from fresh leaves of C. longepaniculatum via GC-MS. However, other researchers [19] used water extracts to identify 23 compounds that were detected from C. camphora leaves. Therefore, different extraction methods could affect the identification of substances.
In recent years, much research has involved the use of transcriptomic data to investigate the key regulatory genes related to terpenoid biosynthesis in C. camphora [21,35]. Many terpenoid synthases have been isolated and characterized from various plant species [12]. In our research, our transcriptome analysis revealed that 23 DEGs were related to secondary metabolites, of which 12 DEGs were related to the metabolism of terpenoids and polyketides. The genes encoding S-(+)-linalool synthase (TRINITY_DN10559_ c0_g1) and geranylgeranyl diphosphate synthase (TRIN-ITY_DN70426_c0_g1) were differentially expressed in the different varieties (Fig. 6). The DEGs were all related to these two pathways, indicating that the essential oils in C. longepaniculatum metabolize are produced through both the MVA and the MEP pathways, which is consistent with the results of Chandran's study [55]. Interestingly, a large number of DEGs were related to photosynthesis and oxidative phosphorylation. Like for flavonoids, we speculate that differences in light responses between the different varieties led to differences in terpenoid metabolism. C. longepaniculatum leaves synthesize a large number of secondary substances, such as flavonoids and terpenoids, to cope with light stress.

Conclusions
In summary, the different varieties of C. longepaniculatum differ in essential oil content and phenotype. By performing transcriptome and metabolome analyses, we investigated the accumulation of secondary metabolites in different varieties of C. longepaniculatum. Of the 34,503 total genes discovered, we identified 1486 DEGs. KEGG analysis showed that phenylpropanoid, tryptophan biosynthesis and terpenoid pathways were significantly altered between the two varieties. In total, we detected 196 metabolites from C. longepaniculatum. We speculate that light may cause differences in secondary metabolites and phenotypes in the leaves of C. longepaniculatum.

Plant materials
The two varieties of C. longepaniculatum (Gamble) N. Chao ex H. W. Li were identified by Professor Ruizhang Feng (Yibin University, Yibin, PR China). Voucher specimens (CLH, 20171026YBU016; CLL, 20171026YBU014) were deposited in the herbarium of the Faculty of Agriculture, Forestry and Food Engineering, Yibin University, Yibin, PR China. The leaves of different C. longepaniculatum varieties (CLH, CLL) were used as experimental materials. The two varieties were propagated from one mother tree via cuttings. C. longepaniculatum was planted on Hongyan Mountain (27°50′ N, 105° 20′ E) in Yibin, Sichuan Province, PR China. No specific permits were required to collect the plant materials because Hongyan Mountain is a germplasm resource nursery of Yibin University. The leaves used for the experiments were collected in November 2019, November 2020 and November 2021. For each variety, 10 healthy plants that have been growing for 20 years were randomly selected. One hundred functional leaves were selected from the sun-facing and shaded sides of the trees. All the leaves of the same variety were mixed and then divided into three groups as 3 biological replicates. Some of the fresh materials were cut into small pieces, immediately frozen in liquid nitrogen and stored at − 80 °C until further transcriptome and metabolome analysis. The others materials were air dried at room temperature and then subjected to distillation to obtain essential oils.

Essential oil isolation
The leaves used for essential oil extraction were collected in November 2019, November 2020 and November 2021. After air drying, 100 g of leaves was put into a steam distillation instrument for 1.5 h to isolate essential oils. The isolation of essential oils from each sample was repeated three times. The essential oil yield was estimated on a dry weight basis (w/w).

Cytological observations of leaves
Leaves of different varieties were harvested in November 2020, kept in formaldehyde-acetic acid-ethanol (FAA) (70% ethanol:acetic acid:38% methanol = 90:5:5), and then embedded in paraffin after dehydration through a series of ethanol concentrations according to the manufacturer's protocol. Cross-sections of 2 to 5 mm thickness were cut with a microtome (Leica RM2235). The crosssections were double dyed with safranine and fast green and covered with a slide cover. The sections of the leaves were observed under a microscope (Motic BA200). The leaf cross-sectional area and leaf thickness were subsequently determined with ImageJ (image processing and analysis in Java).

RNA extraction, library construction and sequencing of fresh leaves
The leaves of different varieties collected in November 2020 were used for RNA-seq. Total RNA was obtained using an RNA purification reagent for plant tissue (Invitrogen, Carlsbad, CA, USA). Three biological replicates were sampled per variety. Six RNA-seq transcriptome libraries were prepared using an Illumina TruSeq ™ RNA Sample Preparation Kit (San Diego, CA). Poly(A) mRNA was purified from the total RNA using oligo-dT-attached magnetic beads and then fragmented with fragmentation buffer. Taking these short fragments as templates, double-stranded cDNA was synthesized using a SuperScript Double-Stranded cDNA Synthesis Kit (Invitrogen, CA) with random hexamer primers (Illumina). Then, the synthesized cDNA was subjected to end repair, phosphorylation and polyadenylation according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 200-300 bp on 2% low-range ultra-agarose followed by PCR amplification using Phusion DNA polymerase (New England Biolabs, Boston, MA), with 15 PCR cycles. After quantification by a TBS-380 instrument, 6 RNA-seq libraries were sequenced in a single lane on an Illumina HiSeq Xten/NovaSeq 6000 sequencer (Illumina, San Diego, CA) for 2 × 150 bp paired-end reads.

De novo assembly and annotation
The raw paired-end reads were trimmed and subjected to quality control by SeqPrep (https:// github. com/ jstjo hn/ SeqPr ep) and Sickle (https:// github. com/ najos hi/ sickle), with the default parameters. Then, clean data from the samples (C. longepaniculatum) were utilized to perform de novo assembly with Trinity (http:// trini tyrna seq. sourc eforge. net/) [56]. All the assembled transcripts were searched via BLASTX against the sequence information within the NCBI protein nonredundant (NR), COG, and KEGG databases to identify the proteins whose sequences were most similar to those of the given transcripts to retrieve their functional annotations, and a typical cut-off E-value of less than 1.0 × 10 − 5 was set. The BLAST2GO (http:// www. blast 2go. com/ b2gho me) [57] program was used to obtain GO annotations of unique assembled transcripts for describing biological processes, molecular functions and cellular components. Metabolic pathway analysis was performed using the KEGG database (http:// www. genome. jp/ kegg/) [58].

Differential expression analysis and functional enrichment
To identify DEGs between the two varieties, the expression level of each transcript was calculated according to the transcripts per million reads (TPM) method. Afterwards, RSEM (http:// dewey lab. biost at. wisc. edu/ rsem/) [59] was used to quantify gene abundance. Specifically, differential expression analysis was performed using DESeq2 [60]/EdgeR [61] with a Q value ≤0.05. Genes with |log2(FC)| > 1 and Q value <= 0.05 (DESeq2 or EdgeR)/Q value <= 0.001 (DEGseq) were considered to be significantly differentially expressed. In addition, functional enrichment analysis via the GO and KEGG databases was performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways with a Bonferroni-corrected P value of ≤0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by GOATOOLS (https:// github. com/ tangh aibao/ Goato ols) and KOBAS (http:// kobas. cbi. pku. edu. cn/ home. do), respectively [62,63].

Extraction of metabolites from fresh leaves and GC-MS analysis
The leaves of different varieties collected in November 2020 were used for metabolite extractions. Three biological replicates were sampled per variety. Fifty milligrams of leaf tissue was extracted in 500 μL of 75% methanol, 2% L-2 chlorophenylalanine and 200 μL of chloroform were added, and the mixture was homogenized at 50 Hz for 3 min at − 10 °C. The mixture was then subjected to ultrasonic treatment at 40 kHz for 10 min at 5 °C after vortexing for 30 s, the process of which was repeated 3 times. After being allowed to settle at − 20 °C for 30 min, the sample was centrifuged at 12000 rcf at 4 °C for 20 min. The supernatant was then vacuum dried. A total of 80 μL of methoxyamine hydrochloride (15 mg/mL in pyridine) was added to the sample, which was then shaken for 2 min and incubated at 37 °C for 90 min. Finally, 80 μL of bis-(trimethylsilyl) trifluoroacetamide (BSTFA) with 1% trimethylchlorosilane (TMCS) and 20 μL of n-hexane were added, after which the sample was stored at 70 °C for 60 min after shaking for 2 min. The samples were subsequently incubated at room temperature for 30 min and then analysed by GC-MS. An extract of each sample was prepared as a quality control (QC) sample. The QC samples were compiled and tested in the same manner as the analytic samples. To assess the repeatability of the analysis, the QC samples were injected at regular intervals (every 3 samples).
GC-MS analysis was conducted on an Agilent 8890B gas chromatograph coupled to an Agilent 5977B mass selective detector (inert electron impact ionization (EI) source and ionization voltage of 70 eV) (Agilent, USA). The gas chromatograph was equipped with an HP-5MS (30 m × 0.25 mm × 0.25 μm) capillary column, and the gas chromatograph column temperature was programmed to hold at 60 °C, increase to 310 °C at a rate of 8 °C/m and hold at the final temperature for 6 min. The gas chromatograph conditions included an inlet temperature of 310 °C, helium carrier gas at a constant flow rate of 1 mL/min, a resting oven temperature of 60 °C and a GC-MS transfer line temperature of 300 °C. After a sample injection of 1 μL, the sample was introduced in splitless mode with an inlet temperature of 260 °C. The ion source temperature was 230 °C, and the quadrupole temperature was 270 °C. Data acquisition was conducted in full-scan mode with a range of m/z 50-500.

Metabolome: data processing and statistics
The GC-MS data were analysed with MassHunter Workstation Quantitative Analysis (version 10.0.707.0). An internal standard was used for data QC (reproducibility). Samples with metabolite features and a relative standard deviation (RSD) of QC > 30% were discarded. A multivariate statistical analysis was performed using ropls (version 1.6.2, http:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ ropls. html). All of the metabolite variables were scaled to unit variances prior to conducting principal component analysis (PCA). VIPs were calculated via a OPLS-DA model, and P values were estimated with paired Student's t test via single-dimensional statistical analysis. The metabolites that differentially accumulated between the two groups were mapped to their biochemical pathways through metabolic enrichment and pathway analysis based on a database search (KEGG; http:// www. genome. jp/ kegg/). Scipy.stats (a Python package) (https:// docs. scipy. org/ doc/ scipy/) was exploited to identify a statistically significantly enriched pathway using Fisher's exact test.

Combined transcriptome and metabolome analyses
Pearson correlation coefficients were calculated for metabolome and transcriptomic data integration. Based on the Pearson correlation coefficient, correlations between genes and metabolites in samples can be measured, and the range of the correlation coefficient is (− 1, + 1). A correlation coefficient less than 0 indicates a negative correlation; when it is greater than 0, it indicates a positive correlation; and when it is equal to 0, it indicates no correlation. The threshold of the correlation coefficient in this study was ±0.8, and the correlation was significant (P < 0.05). The DEGs and DAMs related to essential oil biosynthesis were used to construct coexpression networks, and the candidate target genes were visualized by Cytoscape (version 3.7.2, USA).

Validation of RNA-Seq data by qRT-PCR
Total RNA was isolated from approximately 100 mg of fresh leaves of the different varieties of C. longepaniculatum, which were collected in November 2020, using an RNAprep Pure Plant Plus Kit (TIANGEN, Beijing, China). One microgram of total RNA was subjected to reverse transcription using a SYBR Green PCR Master Mix (TaKaRa) Kit with gDNA Eraser (Perfect for Real-Time). Real-time PCR was carried out by using SYBR Premix Ex Taq II (TaKaRa) on an ABI StepOne ™ Plus Real-Time PCR System (Roche, Switzerland). All the primers used for qRT-PCR are listed in Table S5. Relative gene expression was calculated using the 2 -ΔΔCt method, and owing to its stable expression via PCR amplification, the ACT (KM086738.1) gene was used as a reference control gene [22].

Statistical analysis
The control and treatment groups were analysed for statistically significant differences using one-way ANOVA followed by Duncan's multiple comparisons test. All the calculations were performed using SPSS software (version 21; IBM, Armonk, NY, USA), and all the results are presented as the mean ± SD of 3 independent biological replications. The treatment means were separated by Duncan's multiple range test at a P value less than 0.01. We then used the min-max normalization method of TBtools [63] to analyse the transcriptomic and metabolomic data representing the expression values used in our heatmap.