Transcriptional responses for biosynthesis of flavor volatiles in methyl jasmonate-treated Chrysanthemum indicum var. aromaticum leaves

Industrial Crops & Products 147 (2020) 112254

Contents lists available at ScienceDirect

Industrial Crops & Products journal homepage: www.elsevier.com/locate/indcrop

Transcriptional responses for biosynthesis of flavor volatiles in methyl jasmonate-treated Chrysanthemum indicum var. aromaticum leaves

T

Wenjie Gaoa,e, Qingran Mengb, Hong Luoc, Feng Chend, Yunwei Zhoue,*, Miao Hee,* a

Department of Ecological Technology and Engineering, Shanghai Institute of Technology, Shanghai, 201418, PR China School of Perfume and Aroma Technology, Shanghai Institute of Technology, Shanghai, 201418, PR China c Department of Genetics and Biochemistry, Clemson University, 110 Biosystems Research Complex, Clemson, SC, 29634, USA d Department of Food, Nutrition and Packaging Sciences, Clemson University, Clemson, SC 29634, USA e Department of Landscape Architecture of Northeast Forestry University, Harbin, Heilongjiang, 150040, PR China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Chrysanthemum indicum var. aromaticum Methyl jasmonate Transcriptome Differentially expressed genes Secondary metabolism

Chrysanthemum indicum var. aromaticum (Ci. aromaticum) is an important aromatic plant in genus Chrysanthemum for its extraordinary flavors and unique aroma. Methyl jasmonate (MeJA) has been successfully used as an effective elicitor to improve the quality of aroma by promoting flavor volatiles in Ci. aromaticum leaves. Previous studies focused mainly on the extraction of essential oils as well as their aromatic components from gelatinous secretion of glandular hair of the Ci. aromaticum. However, there has been no study on their aroma formation at the molecular level. Therefore, the information about the transcriptome of Ci. aromaticum leaves, and specifically, the description of changes in gene expression in response to MeJA, is needed for a better understanding of the biological mechanisms of the MeJA-mediated biosynthesis of volatiles in Ci. aromaticum. Transcriptome sequencing was performed on leaves of Ci. aromaticum under different induction treatments of MeJA. A total of 176,091 unigenes were obtained, and 69.94 % of them could be annotated in at least one database. In the comparison groups of 0.25 % MeJA-treated for 4 h vs control (MJ4vsCtrl), 0.25 % MeJA-treated for 24 h vs control (MJ24vsCtrl), and 0.25 % MeJA-treated for 24 h vs 0.25 % MeJA-treated for 4 h (MJ24vsMJ4), 750, 650, and 313 DEGs were obtained, respectively. These DEGs were enriched in 33 secondary metabolic pathways. After analyses of the expression heatmaps of DEGs involved in jasmonic acid biosynthetic and signal transduction pathway, terpenoid biosynthetic pathway, phenylpropanoid metabolism pathway, fatty acid metabolism pathway, as well as differentially expressed post-plant modification enzymes, and transcription factors related to secondary metabolism, we found that most of the DEGs were up-regulated by the MeJA. The results of qRT-PCR were basically consistent with those of transcriptome sequencing, indicating that the gene expression profile data of transcriptome sequencing were reliable. This is the first de novo transcriptome sequencing reported in Ci. aromaticum. The transcriptome analysis showed that exogenous application of MeJA could induce expression of genes in signal transduction and secondary metabolism, especially for the biosynthetic pathways of flavor volatiles. The present results have provided some insights of the mechanisms how MeJA activates the biosynthetic pathways of volatile metabolites in the Ci. aromaticum leaves.

Abbreviations: Ci. Aromaticum, Chrysanthemum indicum var. aromaticum; MeJA, Methyl jasmonate; DEGs, Differentially expressed genes; SDB, Scopadulcic acid B; GGDP, Geranylgeranyl diphosphate; RNA-Seq, RNA sequencing; Nr, Non-redundant; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene ontology; TF, Transcription factor; qRT-PCR, Quantitative real-time PCR; UDPG, UDP-glycosyl- transferase; JA, Jasmonic acid; AOC, Allene oxide cyclase; ACX, Acyl-CoA oxidase; ACAA1, Acetyl-CoA acyltransferase 1; PMVK, Phosphomevalonate kinase; ISPD, 2-C-Methyl-d-erythritol-4-phos-phatecyti-dyltransferase; ISPF, 2-C-Methyl-d-erythritol 2,4-cyclodiphosphate synthase; ISPG, 1-hydro- xy-2-methyl-2-(E)-butenyl-4-diphosphate synthase; MEP, 2-C-Methyl-D -Erythritol-4-Phosphate pathway; IPP, Isopentenyldiphosphate; DMAPP, Dimethylallyldiphosphate; GPS, Geranyl diphosphate synthase; FPS, Pyrophosphate synthase; GPP, Geranyl pyrophosphate; FPP, Farnesyl pyrophosphate; R-LIS, R-linalool synthase; LIS, Linalool synthase; HPT, Phytylprenyltransferase; TAT, Tyrosine aminotransferase; C4H, Cinnamic acid-4hydroxylase; HPL, Hydroperoxidelyase; ADH, Alcohol dehydrogenase; FAD2, Fatty acid desaturase; CYP450s, Cytochrome P450s; GSS, Genome survey sequences; TPS, Terpene synthase ⁎ Corresponding authors at: Department of Landscape Architecture of Northeast Forestry University, No. 26 Hexing Road, Harbin, Heilongjiang, 150040, PR China. E-mail addresses: [email protected] (Y. Zhou), [email protected] (M. He). https://doi.org/10.1016/j.indcrop.2020.112254 Received 6 November 2019; Received in revised form 22 January 2020; Accepted 19 February 2020 0926-6690/ © 2020 Elsevier B.V. All rights reserved.

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

1. Introduction

2. Materials and methods

Plants used primarily for their aromatic properties in perfumery are defined as aromatic plants (Ramya et al., 2013). Aromatic plants are widely used in landscaping and molecular breeding because of their unique fragrance. Besides the aromatic properties, many aromatic plants are also used for medicinal and cosmetic purposes. In terms of production of industrial products and fine chemicals from plants, some aromatic plants can also be used to produce colorants, dyes, and crop protection products (Lubbe and Verpoorte, 2011). Chrysanthemum, one of the top four plant species in the world for cut flowers and the top ten most famous traditional Chinese flowers, plays an important role in the modern flower production with both economic and ornamental values (Lu et al., 2018). The chrysanthemum cultivars exist in various colors and flower types, though the flowers generally have subtle scent or even scentless, rather than the strong fragrance. In addition, most chrysanthemum varieties are suitable for planting in pots but not in the open field. Hence, the development of chrysanthemum with special characteristics is a top priority research in Chrysanthemum breeding (He et al., 2016). Chrysanthemum indicum var. aromaticum (Ci. aromaticum) is an herbaceous perennial plant of the Compositae family that is commonly found in Shennongjia, Hu Bei province, P.R. China (Qihong and Shufan, 1983). The plant is an important source of aromas in the genus Chrysanthemum due to its special scent in its flowers, leaves, and roots (Gao et al., 2018). The leaves and petals of Ci. aromaticum are densely covered with glandular trichomes, and the secretions (especially the composition and the content) from glandular trichomes are closely related to the aroma components of Ci. aromaticum (Jin et al., 2017). Methyl jasmonate (MeJA) belongs to a kind of aromatic compound. It was originally isolated from the essential oil of Jasminum grandiflorum and is a class of plant stress hormones similar to the salicylates (Cohen and Flescher, 2009). As a plant stress hormone, MeJA is critical to plant growth and development, stress response, and secondary metabolic processes (Zang et al., 2015). Cho et al. (2008) found that exogenous MeJA could act as a signal transduction molecule: it mediated the production of secondary metabolites such as phenols and terpenoids, thereby boosting contents of natural plant chemicals (Cho et al., 2008). It has been reported that MeJA could selectively stimulate the biosynthesis-related enzymes for scopadulcic acid B (SDB) (such as geranylgeranyl diphosphate (GGDP) synthetase and GGDP cyclase) in wild licorice culture tissues, which in turn could promote the accumulation of SDB (Nkembo et al., 2005). Ge and Wu (2005) observed that, compared with the controls, the treatment of Salvia miltiorrhiza hairy roots with 100 μmol/L MeJA led to a significant increase in tanshinone content, an elevated expression of HMGR and DXR, two key enzyme genes crucial to the biosynthesis of steroids (Ge and Wu, 2005). Hence, exogenous application of MeJA provides an effective way to explore the key enzyme genes in the process of plant secondary metabolism and clarify the related pathways, and the production of plant secondary metabolites could also be further regulated. Our previous study showed that MeJA treatment could induce changes in the number of epidermal hairs and the content of volatile terpenoids in Ci. aromaticum leaves (He et al., 2016). In present study, RNA-Seq technique was adopted in detection of the differences in the regulation of terpenoids metabolism related genes in Ci. aromaticum leaves under the MeJA treatment at different times, which was important for us to understand the secondary metabolism-related genes, especially genes involved in the synthesis of terpenoids. The information obtained from the present study is expected to provide a basis for the development of genetic engineering technology to generate terpenoids-rich aromatic materials, and to produce new variety of chrysanthemum with predominance aroma possible.

2.1. Plant material and MeJA treatment Ci. aromaticum was introduced from the Shennongjia area of Hubei Province (110°23′57″ E, 31°28′7″ N), and was successfully cultivated in the nursery garden of Northeast Forestry University. Robust foot buds (3–5 cm) of Ci. aromaticum were cut and inserted into nutrient medium (nutrient soil: vermiculite = 1 : 1) in the greenhouse for conventional management (16 h light and 8 h dark, the light intensity was 300 μE m−2 s-1, 23 ± 2 °C). After about 2.5 months, 45 pots of Ci. aromaticum plants with good condition and consistency (plant height was around 25 cm) were selected. The Ci. aromaticum plants were divided into 3 groups, with 15 plants on each group. The three groups were placed in a closed chamber (1.0 m × 1.0 m × 0.8 m) separated by plastic film in the plant culture room. The distribution of the three groups was as follows: Group 1: plants were sprayed with 1.4 % ethyl alcohol solution, and the total spraying volume was 100 mL; Group 2 and Group 3: plants were sprayed with 0.25 % MeJA solution dissolved in 1.4 % ethyl alcohol solution, and the total spraying volume was 100 ml. After the treatments, the three groups were placed in the plant culture room for routine management. Plant leaves (the third to fifth fully expanded leaves) were collected after being treated for 4 h and 24 h, respectively, then quick-frozen in liquid nitrogen and stored at −80 °C for further use. 2.2. Total RNA extraction, cDNA library construction and sequencing Leaf total RNA samples from the Ctrl (treatment with 1.4 % ethyl alcohol solution as control), the MJ4 (treatment with 0.25 % MeJA solution dissolved in 1.4 % ethyl alcohol solution for 4 h) and the MJ24 plants (treatment with 0.25 % MeJA solution dissolved in 1.4 % ethyl alcohol solution for 24 h) were extracted using Trizol reagent according to the manufacturer’s manual. The quality and integrity of the total RNA were evaluated using 0.8 % agarose gel electrophoresis and a 2100 Bioanalyzer RNA Nano chip device (Agilent, Santa Clara, CA, USA), respectively. The concentrations were measured using a ND-1000 spectrophotometer (NanoDrop, Thermo Scientific, DE, USA). RNA samples in line with the quality standard (OD260/OD280 ≥ 2.0) were treated with DNase I and then sent to the Beijing Genomics InstituteShenzhen (BGI, Shenzhen, China) for mRNA purification, cDNA library construction and sequencing using the Illumina technology. 2.3. De novo transcriptome assembly and annotation In order to ensure the quality of information analysis, raw reads were processed and clean reads without adaptor sequences, repeat, or low-quality ones (Q ≤ 5) were obtained. Trinity (http://trinityrnaseq. sourceforge.net/) was used to assemble the short reads into contigs. The homologous transcript clustering was performed using Corset software (Davidson and Oshlack, 2014), and the longest cluster obtained after clustering was used for subsequent analysis. The cluster sequences obtained after assembly were functionally annotated using the BLASTx program (http://www.ncbi.nlm.nih.gov/BLAST/) with the E-value < 1 × 10−5. Databases for annotating Ci. aromaticum transcripts included the NCBI non-redundant (Nr) databases (http://www.ncbi.nlm.nih. gov/COG), Swiss-prot databases (http://www.expasy.ch/sprot), Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (http://www. genome.jp/kegg), and gene ontology (GO) databases (http://www. geneontology.org/). 2.4. DEGs functional annotation and enrichment analysis Differential gene expression analysis between the sample groups was performed using DESeq2 (Love et al., 2014), and the gene with |log2FoldChange| > 1 and P-adj < 0.05 (P-adj was the corrected p 2

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

splice out 176,091 unigenes with an average read length and N50 lengths of 1196 and 1722 bp, respectively. The length distribution of unigenes was shown in Fig. 1. The results indicated that the number of the unigene fragments with the length of less than 300, 300–500, 500–1000, 1000–2000, and greater than 2000 bp were 14,133 (8.02 %), 33,446 (19 %), 47,393 (26.91 %), 51,293 (29.13 %), and 29,826 (16.94 %), respectively.

value) was identified as a significant DEG. The number of common and unique DEGs between the sample groups were presented by means of Venn diagram. K-means hierarchical cluster analysis was used to analyze the DEGs set. GO enrichment analysis on the differential genes was performed using GOseq (Young et al., 2010) method, and KEGG pathway enrichment analysis was performed by comparing the number of genes in each pathway and their background of the differential genes.

3.2. Functional annotation of unigenes 2.5. Transcription factors analysis For functional annotation analysis, all the assembled unigenes were searched against the public databases including Nr, Nt, KO, SwissProt, PFAM, GO, and KOG using the BLASTx program. The results in Fig. 2 showed that 108,988 (61.89 %) unigenes could be annotated in the Nr database, 55,932 (31.76 %) in the Nt database, 41,618 (23.63 %) in the KO database, 84,341 (47.89 %) in the SwissProt database, 84,680 (48.08 %) in the PFAM database, 85,585 (48.6 %) in the GO database, and 27,593 (15.66 %) in the KOG database. Besides, there were 13,480 (7.65 %) unigenes annotated in each database, and 123,175 (69.94 %) in at least one database. The species statistics of the unigenes annotated by the Nr database were conducted, and the top five annotated species were Vitis vinifera, Sesamum indicum, Coffea canephora, Nicotiana sylvestris, and Nicotilana tomentosiformis, accounting for 11.7 %, 7.2 %, 6.9 %, 4.5 %, and 4.4 % of the total number of Nr annotations, respectively (Fig. 3).

Plant transcription factors were predicted using iTAK software, and the transcription factor (TF) was identified by hmmscan according to the TF family and rules defined in the database. 2.6. Quantitative real-time PCR (qRT-PCR) analysis To validate the results of RNA-Seq, 16 DEGs which were related to the terpenoid biosynthetic pathway genes, terpenoid synthase genes, and transcription factor that regulated the synthesis of terpenoids were chosen for qRT-PCR. The reverse transcription of RNA-seq RNA sample with same batch were performed using a reverse transcription kit (Toyobo, Japan). qRT-PCR was run using SYBR premix Ex Taq Kit (Takara, Japan) and a Bio-Rad MiniOpticon Real-Time PCR machine (Bio-Rad, USA) using the same cDNA samples as that in the RNA-seq experiment. The qRT-PCR reaction system with the total volume of 20 μl contained 10 μl of SYBR MIX, 0.4 μl of upstream primers, 0.4 μl of downstream primers, 0.5 ng of cDNA template, and the amplification program was as follows: pre-denaturation at 94 °C for 30 s, followed by 40 cycles of denaturation at 94 °C for 5 s, renaturation at 60 °C for 15 s, and extension at 72 °C for 10 s. The Chrysanthemum CmEF1α (KF305681) was used as the internal reference gene, and the relative expression levels were calculated using the 2−ΔΔCT method (Livak and Schmittgen, 2001). Three independent biological replicates of each sample and three technical replicates of each biological replicate were used for the qRT-PCR analysis. Based on the sequences, a set of genespecific primer pairs was designed using Primer Premier 5.0 software, which were shown in Additional file 1: Table S1.

3.3. GO classification of unigenes Based on sequence homology, 85,585 unigenes were classified against GO terms to classify the predicted function genes of Ci. aromaticum. As shown in Fig. 4, the functionally annotated unigenes were categorized into 57 functional group in the three categories of biological processes, cellular components, and molecular function (Fig. 4. Wherein the categories of biological processes with the largest amount of annotated genes were “cellular process” 47,114 unigenes, “metabolic process” 44,257 unigenes, “single-organism process” 33,973 unigenes, and “biological regulation” 15583 unigenes, suggesting that the present study may be helpful in identifying new genes involved in the secondary metabolite synthesis pathway of Ci. aromaticum. Among the cell compositions, the main representative classes were “cells” (23, 915 unigenes), “cell parts” (23, 903 unigenes), “organelles” (16,104 unigenes), and “macromolecular complexes” 15,307 unigenes. In the molecular function classification, it mainly includes “binding” 52,769 unigenes, “catalytic activity” 39,815 unigenes, “transporter activity” 5,590 unigenes, and “nucleic acid binding transcription factor activity” 2,349 unigenes.

2.7. Statistical analysis All the data were analyzed using Microsoft Office Excel 2019 software (Microsoft Corporation, Redmond, WA, USA) with one-way ANOVA. Data were evaluated by the analysis of variance and significance by Duncan’s multiple range comparison test by using SPSS ver. 18.0 software (SPSS, Inc., Chicago, IL, USA). Differences were considered significant at P < 0.05 or P < 0.01.

3.4. KEGG pathways classification of unigenes 3. Results To identify the biological pathways that were activated in Ci. aromaticum, all unigenes were mapped in the KEGG database and 41,618 unigenes were attributed to 129 metabolic pathways. Among them, as shown in the Table 2, there were 9 metabolic pathways that are related to the biosynthesis of terpenoids, which were “terpenoid backbone biosynthesis” (469 unigenes), “monoterpene biosynthesis” (110 unigenes), “diterpenoid biosynthesis” (83 unigenes), “sesquiterpene and triterpenoid biosynthesis” (228 unigenes), “ubiquinone and other terpenoid-quinone biosynthesis” (210 unigenes), “zeatin biosynthesis” (133 unigenes), “carotenoid biosynthesis” (212 unigenes), “steroid biosynthesis” (179 unigenes), and “N-glycan biosynthesis” (237 unigenes), respectively. The above steroid related unigenes accounted for 4.47 % of all annotated unigene e ratios. In addition, many unigenes had been annotated to other secondary metabolic pathways such as flavonoid biosynthesis -related pathways, phenylpropanoid-related pathways, and alkaloid-related pathways. The annotations of these unigenes showed a significant transcriptional complexity and provided

3.1. Transcriptome sequencing and de novo assembly of Ci. aromaticum RNA-seq on Ci. aromaticum control with different treatments was performed using the Illumina HiSeq™ 4000 sequencing platform, which resulted into three samples named Ctrl, MJ4, and MJ24 with a total of 9 sequencing libraries. The results in Table 1 showed that at least 48,775,846 Raw Reads were obtained in each library, and at least 47,681,166 Clean Reads and 7.15 G Clean Bases were obtained for each library after removing the sequencing linker, low-quality Reads, and Reads with N ratio greater than 10 %. The percentage of Clean Reads in all nine libraries was over 97 %, the percentage of Q30 was at least 94.97 %, and the percentage of GC was about 44 %. The Trinity software was used to splice the Clean Reads of Ci. aromaticum in 9 libraries, leading to 212,966 transcripts. The average read length and N50 length were 1036 and 1656 bp, respectively. Trinity was continued to be used to splice those transcripts together to 3

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Table 1 Summary of Illumina HiSeq™ 4000 sequencing data. Samples

Raw reads

Clean reads

Ctrl_1 Ctrl_2 Ctrl_3 MJ4_1 MJ4_2 MJ4_3 MJ24_1 MJ24_2 MJ24_3

65,329,272 53,648,624 61,556,944 58,908,116 48,775,846 66,693,480 50,654,828 58,032,074 51,285,780

63,888,772 52,446,220 60,259,008 57,669,110 47,681,166 65,286,622 49,537,918 56,715,224 50,176,564

a b c

(97.80 (97.76 (97.89 (97.80 (97.73 (97.84 (97.90 (97.76 (97.89

%) %) %) %) %) %) %) %) %)

a

Clean bases (bp)

Error (%)

Q20

9.58G 7.87G 9.04G 8.65G 7.15 G 9.79G 7.43G 8.51G 7.53G

0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02

98.49 98.45 98.45 98.35 98.47 98.36 98.46 98.41 98.51

(%)

Q30

b

(%)

95.32 95.23 95.24 94.97 95.28 95 95.24 95.13 95.37

GC content

c

(%)

44.74 45.18 45.05 43.89 43.89 44.52 44.4 44.13 44.07

Q20: Percentage of bases that account for more than 20 of Phred values. Q30: Percentage of bases that account for more than 30 of Phred values. GC content: The sum of base G and C accounts for the percentage of total base number.

Fig. 3. Species distribution of the BLASTX result against the Nr database. Fig. 1. Length distribution of assembled Transcripts and unigenes.

3.6. Differentially expressed genes (DEGs) analysis To identify DEGs among the MeJA-treated Ci. aromaticum samples, we compared them with each other and identified significant DEGs with P-adj < 0.05 and |log2FoldChange| > 1, and the hierarchical clustering was used to gain a global view of DEGs (Fig. 5). The results showed that the MJ4vsCtrl comparison group identified 750 DEGs, of which 537 were up-regulated and 213 were down-regulated; the MJ24vsCtrl comparison group identified 650 DEGs, with 424 up-regulated and 226 down-regulated; the MJ24vsMJ4 comparison group revealed 313 DEG, and 150 were up-regulated and 163 were down-regulated (Fig. 5A–C). As can be seen from Fig. 5D, there were 9 DEGs shared by the three comparison groups, and all of them were up-regulated. Among the 9 DEGs, three genes were annotated as secondary metabolism related enzyme genes, namely, UDP-glycosyltransferase (UDPG), β-fructofuranosidase, and aurone synthase. There were 197 Shared DEGs in the MJ4vsCtrl and MJ24vsCtrl comparison groups, and 486 and 398 exclusive DEGs in the MJ4vsCtrl and MJ24vsCtrl comparison groups, respectively. In order to visualize the DEGs of Ci. aromaticum under the MeJA processing, we generated a heat map by hierarchical clustering (Fig. 5E). The heat map showed that the up- and down-regulated DEGs in the MJ4 and MJ24 sample groups had either the same or different expression trend, indicating that Ci. aromaticum with the MeJA treatment for 4 h and 24 h had both common and unique gene expression patterns.

Fig. 2. Functional annotations of Ci. aromaticum in seven databases.

information about the expression of valuable genes in the transcriptome of Ci. aromaticum.

3.5. Transcription factor function annotation A total of 4697 transcription factors annotated by searching the predicted CDS in the plant transcription factor database PlantTFdb belong to 79 transcription factor families, of which the majority are those annotated as the MYB family transcription factors (366, 7.79 %), followed by the AP2-EREBP family (268, 5.71 %), C3H family (247, 5.26 %), WRKY family (237, 5.05 %), bHLH family (227, 4.83 %), Orphans family (209, 4.45 %), HB family (182, 3.87 %), RWP-RK family (179, 3.81 %), C2H2 family (169, 3.6 %), and ABI3VP1 family transcription factors (148, 3.15 %) (Table 3).

3.7. GO enrichment analysis of DEGs To further explore the possible functions of DEGs, the GO enrichment analysis on DEGs was performed. Fig. 6 showed that 577 of the 750 DEGs, which included 1104 biological process terms, 259 cellular component terms, and 485 molecular functional terms, were enriched on the GO term in the MJ4vsCtrl comparison group. Similarly, in the 4

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 4. GO classification of all unigenes. Annotated sequences were classified into ‘Biological Process’, ‘Molecular Function’ and ‘Cellular. Component’ groups and 57 subgroups.

group, 236 of the 313 DEGs including 796 biological process terms, 200 cell-constituting terms, and 326 molecular functional terms were enriched on the GO terms. The expressed genes with significant difference in the MJ24vsMJ4 comparison group were particularly enriched in the carbohydrate metabolic process (33, 13.98 %), glucan metabolic process (6, 6.36 %), polysaccharide metabolic process (6, 6.36 %), and so on. In the MJ4vsCtrl, MJ24vsCtrl and MJ24vsMJ4 comparison groups, there were five GO terms associated with the synthesis and metabolism of terpenoids, which were terpenoid biosynthetic processes (4, 3, and 3, respectively), terpenoid metabolic process (4, 3, and 3, respectively), tetraterpenoid biosynthetic process (2, 1, and 1, respectively), tetraterpenoid metabolic process (2, 1, and 1, respectively), and terpene synthase activity (7, 10, and 5, respectively) (Table 4). In addition, in the MJ24vsMJ4 comparison group, there was one GO term separately associated with sesquiterpenoid biosynthetic process, sesquiterpenoid metabolic process, and sesquiterpene synthase activity.

Table 2 Overviews of several terpenoid biosynthetic pathways and unigenes distribution. Terpenoid related pathways

Unigene number

Percent (%)

Terpenoid backbone biosynthetic pathway Monoterpenoid biosynthetic pathway Diterpenoid biosynthetic pathway Sesquiterpenoid and triterpenoid biosynthetic pathway Ubiquinone and other terpenoid-quinone biosynthetic pathway Zeatin biosynthetic pathway Carotenoid biosynthetic pathway steroid biosynthetic pathway N-Glycan biosynthetic pathway Total

469 110 83 228

1.13 0.26 0.20 0.55

210

0.50

133 212 179 237 1861

0.32 0.509 0.43 0.57 4.47

Table 3 Top 10 transcription factors annotated by Ci. aromaticum transcriptome. Transcription factor families

Unigene number

Percent (%)

MYB AP2-EREBP C3H WRKY bHLH Orphans HB RWP-RK C2H2 ABI3VP1

366 268 247 237 227 209 182 179 169 148

7.79 5.71 5.26 5.05 4.83 4.45 3.87 3.81 3.6 3.15

3.8. KEGG enrichment analysis of DEGs All DEGs were mapped into the KEGG database to search for DEGs involved in metabolic or signal transduction pathways. From Fig. 7, it is obvious that the three comparison groups, i.e., MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4, were mainly enriched in 33 secondary metabolism related pathways, wherein the most enriched pathway for DEGs was the phenylpropanoid biosynthetic pathway with 18, 13, and 1 DEGs in each group, followed by the unsaturated fatty acid biosynthetic pathway (7, 16, 1) and Alpha-linoleic acid metabolic pathway (9, 4, 5). Among them, the biosynthetic pathways of terpenoids were enriched into five, which were terpenoid backbone biosynthetic (1, 5, 2), monoterpene biosynthetic (1, 1, 0), sesquiterpene and triterpenoid biosynthetic (0, 7, 3), diterpene biosynthetic (4, 1, 0), and ubiquinone and other terpenoid-quinone biosynthetic pathways (3, 2, 1). In addition, pathways related to some other secondary metabolites, including the flavonoids, carotenoids, zeatin and alkaloids biosynthetic pathways, acid and ester metabolic pathways, terpene degradation pathways, plant hormone signal transduction pathways, and cytochrome P450 metabolic pathway, were also enriched.

MJ24vsCtrl comparison group, 483 of the 650 DEGs including 939 biological process terms, 215 cellular component terms, and 434 molecular functional terms were enriched on the GO terms. Compared with the Ctrl group, the expressed genes with significant difference after the MeJA treatment for 4 h and 24 h were mainly enriched in catalytic activity (369, 63.95 %; 294, 60.87 %), single-organism metabolic process (231, 40.03 %; 167, 34.58 %), oxidation-reduction process (94, 16.29 %; 77, 15.94 %), and carbohydrate metabolic process (61, 10.57 %; 46, 9.52 %). Besides, in the MJ24vsMJ4 comparison 5

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 5. Analysis of DEGs under MeJA treatment. A: MJ4vsCtrl; B: MJ24vsCtrl; C: MJ24vsMJ4; D: Venn diagram of MJ4vsCtrl_MJ12vsCtrl_MJ12vsMJ4 DEGs; E: Cluster of DEGs.

Compared with the Ctrl group, difference observed in the MeJA treatment for 4 h was more notable, indicating that MeJA treatment for 4 h was more effective than 24 h treatment in inducing the expression of JA biosynthetic and signaling pathway genes in Ci. aromaticum, and the most significant DEGs were JAZ and jasmonate O-methyltransferase, which were up-regulated by 10.18-fold, 7.78-fold, 9.62-fold and 5.99fold in the MJ4 and MJ24 groups, respectively. The fold changes in differential expression for other DEGs were between 0.5 and 4 times.

3.9. DEGs involved in jasmonic acid (JA) biosynthetic and signal transduction pathways Exogenous MeJA is a major regulator of plant JA biosynthetic and JA signaling pathways. In the transcriptome data, 11 DEGs including lipoxygenase 2 (LOX2, 3), allene oxide cyclase (AOC, 1), acyl-CoA oxidase (ACX, 2), enoyl-CoA hydratase (MFP2, 1), acetyl-CoA acyltransferase 1 (ACAA1, 2), and jasmonate O-methyltransferase (2) related to the JA biosynthetic pathways were identified; 13 DEGs (JAZ (10) and MYC2 (3)) were recognized in the JA signaling pathways. Among all the 24 DEGs, 19 DEGs were up-regulated after the MeJA treatment, and 5 were down-regulated. From the heat map, as shown in Fig. 8, the expression trends of these 24 DEGs were basically the same.

3.10. DEGs involved in fatty acid metabolism pathways Among the fatty acid derivatives, saturated and unsaturated volatiles C6 and C9 aldehydes and alcohols are the main components of 6

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 6. GO enrichment analysis of DEGs. A: MJ4vsCtrl; B: MJ24vsCtrl; C: MJ24vsMJ4.

7

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Table 4 Distribution of DEGs enriched in terpenoid biosynthesis and metabolism related GO terms. Description

terpenoid biosynthetic process terpenoid metabolic process tetraterpenoid biosynthetic process tetraterpenoid metabolic process terpene synthase activity sesquiterpenoid biosynthetic process sesquiterpenoid metabolic process sesquiterpene synthase activity

Type

Biological process Biological process Biological process Biological process Molecular function Biological process Biological process Molecular function

MJ4vsCtrl

MJ24vsCtrl

MJ24vsMJ4

DEG/Total DEGs

Percent (%)

DEG/Total DEGs

Percent (%)

DEG/Total DEGs

Percent (%)

4/577 4/577 2/577 2/577 7/577 0/577 0/577 0/577

0.69 0.69 0.35 0.35 1.21 0 0 0

3/483 3/483 1/483 1/483 10/483 0/483 0/483 0/483

0.62 0.62 0.21 0.21 2.07 0 0 0

3/236 3/236 1/236 1/236 5/236 1/236 1/236 1/236

1.27 1.27 0.42 0.42 0.42 0.42 0.42 0.42

3.11. DEGs involved in terpenoids biosynthetic pathways

plant leaf aroma. In this study, the differential gene analysis of the transcriptome of Ci. aromaticum under the treatment of MeJA showed that there were 14 key enzyme genes annotated to the lipoxygenase pathway, namely LOX2 (3), hydroperoxidelyase, (HPL, 1), and alcohol dehydrogenase (ADH, 10). Their expression heat maps showed that all except one LOX2 and HPL were up-regulated after the MeJA treatment, and the differences were more significant in the MJ4 than in the MJ24 group. In addition, 16 fatty acid desaturase (FAD2) enzyme genes on the biosynthetic pathway of unsaturated fatty acids were also annotated. Fourteen of them were up-regulated after the MeJA treatment, and two were down-regulated. The expression patterns of all the FAD2 enzyme genes were consistent, and their induction was more significant in MJ24 than in MJ4group (Fig. 9). The results indicated that exogenous MeJA treatment could induce the expression of fatty acid metabolism pathway genes, and in turn produce more fatty acid derived volatile compounds.

Most of the chemicals identified from the surface secretion of Ci. aromaticum leaves were terpenes, mainly monoterpenes and sesquiterpenoids. Therefore, when conducting transcriptome data analysis, the DEGs on terpenoids biosynthetic pathways were the major concerns. Our results showed that there were 43 DEGs annotated to the terpenoids biosynthetic pathways. Five of these DEGs were involved in the steroid skeleton biosynthetic pathways, of which two were phosphomevalonate kinase (PMVK) genes annotated to the malonate pathway (MVA), and the remaining three genes, 2-C-methyl-d-erythritol-4-phosphatecyti-dyltransferase (ISPD), 2-C-methyl-d-erythritol 2,4-cyclo- diphosphate synthase (ISPF), and 1-hydroxy-2-methyl-2-(E)butenyl-4-diphosphate synthase (ISPG), were annotated to 2-C-methylD-erythritol-4-phosphate pathway (MEP). Compared with the Ctrl group, the longer the treatment time with the MeJA, the higher the

Fig. 7. DEGs analysis of KEGG enrichment to secondary metabolism related pathways. A: Terpene biosynthesis pathway; B: Phenylpropanoid related pathway; C: Flavonoid biosynthesis pathway; D: Other secondary metabolite biosynthesis pathway; E: JA related metabolic pathway; F: Acid and ester metabolic pathway; G: Terpenes degradation pathway; H: Plant hormones and post-modification enzymes metabolic pathway. 8

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 8. DEGs assigned to JA biosynthesis and signaling pathways under MeJA treatment (The expression heat maps from left to right were arranged in the following order: MJ4vsCtrl、MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples.). A: JA biosynthesis pathway; B: JA signaling pathways.

expression of PMVK gene. The ISPD, ISPF, and ISPG genes were upregulated in the MJ4 group, but down-regulated in the MJ24 group, indicating that MJ24 could induce the expression of the MVA pathway enzyme genes, and MJ4 tends to induce the expression of the MEP pathway enzyme genes to generate the terpenoids synthesis precursors (isopentenyldiphosphate, IPP and dimethylallyldiphosphate, DMAPP). Besides, two important genes, geranyl diphosphate synthase (GPS) and farnesyl pyrophosphate synthase (FPS), were also annotated, which were involved in monoterpenes synthesis direct acting precursor geranyl pyrophosphate (GPP) and sesquiterpene and triterpenoids synthesis direct acting precursor farnesyl pyrophosphate (FPP). As shown in Fig. 9, the GPS and FPS gene expression was continuously up-regulated as the processing time increased, thereby facilitating the generation of GPP and FPP. In addition, we also annotated 10 monoterpene synthase genes, including R-linalool synthase (R-LIS, 4), linalool synthase (LIS, 1), (E)-β-ocimene synthase (3), (+)-neomenthol dehydrogenase (1), and β-myrcene synthase (1). All these genes were increased in expression after the MeJA treatment, of which, the genes for (+)-neomenthol dehydrogenase (8.99- and 6.65-fold) and (E)-β- myrcene synthase (6.97-, 4.41-fold) exhibited the most pronounced changes in gene expression. There were 10 sesquiterpene synthase DEGs, including βcaryophyllene synthase (3), cascarilladiene synthase (2), β-amyrin synthase (2), E-β-farnesene synthase (1), and sesquiterpene synthase TPS (2). Seven of them showed expression proportional to the time of MeJA treatment. The genes with the most significant expression change was cascarilladiene synthase (5.46- and 6.9-fold). The expression of CYP82G1(4), described in the diterpene biosynthetic pathway, as well as CYP73A (1), homogentisate phytylprenyltransferase (HPT, 2), and tyrosine aminotransferase (TAT, 1), which belong to ubiquinone and other terpenoid-quinones biosynthetic pathway, was also induced by MeJA, and MJ4 was better than that of MJ24. Moreover, the genes for the strictosidine synthase gene (2), the steroid synthase (3), the carotenoid synthase (1), the zeatin synthase (1), and the N-glycan synthase (1) were also annotated (Fig. 10).

Fig. 9. DEGs assigned to fatty acid metabolism pathway under MeJA treatment (the expression heatmaps from left to right were arranged in the following order: MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples). Red: Lipoxygenase pathway; Blue: Biosynthesis of unsaturated fatty acids.

3.12. DEGs involved in phenylpropanoid metabolic pathways The phenylpropanoid metabolic pathways are important aroma synthesis pathways in plants. The benzene ring and side chains 9

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 10. DEGs assigned to Terpenoids biosynthesis pathway under MeJA treatment (The expression heatmaps from left to right were arranged in the following order: MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples). A: Monoterponoid biosythesis; B: Sesequiterpenoid and triterpenoid biosythesis; C: Steroid biosythesis; D: N-Glycan biosythesis; E: Carotenoid biosythesis; F: Ubiquinone and other terpenoid-quinone biosythesis; G: Diterpenoid biosythesis; H: Indole diterpene alkaloid biosythesis; I: Zeatin biosythesis.

pathways, including anthocyanidin 3-O-glucoside 5-O-glucosyltransferase (4), anthocyanidin 5,3-O-glucosyl- transferase (3), anthocyanin 5-aromatic acyltransferase (2), and anthocyanin malonyltransferase (1). The heat map analysis of these 53 DEGs showed that 49 DEGs were up-regulated and only 4 DEGs were down-regulated under the MeJA treatment, and most of the genes were more differentially expressed in MJ4 group than in MJ24 group (Fig. 11), indicating that MeJA treatment could stimulate the expression of genes related to the phenylpropanoid metabolism pathway of Ci. aromaticum, and the effect of MeJA treatment of 4 h was more significant than that of 24 h.

contained in phenylpropanoid compounds can be modified to form volatile aldehydes, alcohols, ether or esters after hydroxylation, methylation, and acylation. In our transcriptome data analysis, 53 DEGs were annotated to the phenylpropanoid metabolic pathways, of which 21 are involved in the phenylpropanoid biosynthetic pathway, including cinnamic acid-4-hydroxylase (C4H, 1), β-glucosidase (11), shikimate O-hydroxycinnamoyl transferase (7), and peroxidase (2). A total of 14 DEGs participated in phenylalanine, tyrosine, and tryptophan biosynthetic pathways, including 3-deoxy-7-phosphoheptulonate synthase (7), 3-dehydroquinoline dehydratase (3), anthranilate synthase (1), arogenate dehydratase (1), and tryptophan synthase genes (2). There were two DEGs which participated in the phenylalanine metabolism pathway: amidases (1) and aminotransferase (1). In addition, 6 DEGs including flavone synthase II (1), flavonol synthase (2), flavonol sulfotransferase (1), dihydroflavonol 4-reductase (1), and isoflavone reductase (1) were found to be involved in flavonoid biosynthetic pathways, and 10 DEGs participated in anthocyanin biosynthetic

3.13. DEGs encoding post-plant modification enzymes The terpenoids in plants can be hydroxylated, glycosylated, methylated, isomerized epoxidized, added/reduced, halogenated, or rearranged in their framework structures by the action of post-plant modification enzymes, which in turn can greatly contribute to structure 10

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 11. DEGs assigned to Phenylpropanoid metabolism pathway under MeJA treatment (the expression heatmaps from left to right were arranged in the following order: MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples). Red: Phenylpropanoid biosynthesis pathway; Blue: Phenylalanine, tyrosine and tryptophan biosynthesis pathway; Green: Phenylalanine metabolism pathway; Yellow: Flavonoid biosynthesis pathway; Purple: Anthocyanin biosynthesis pathway.

diversity of species and biological activity of terpenoids (Cheng et al., 2007). Cytochrome P450 s (CYP450 s) are hemoglobin enzymes involved in the biosynthesis of some secondary metabolites such as phenylpropanoids, flavonoids, steroids, sterols, and terpenoids (Rai et al., 2015). The results of DEGs annotation showed that there were 23 CYP450 s, of which the 16 DEGs up-regulated after the MeJA treatment, there are 4 CYP450 monooxygenases. The heat map analysis showed that their expression patterns were consistent, and the differential expression was more significant in MJ4 than in MJ24 group. SAM-dependent methyltransferase utilizes methyl groups from SAM and catalyzes the methylation of substrate hydroxyl groups to produce methyl secondary metabolites, such as phenylpropanoids (A.). In present study, a total of 7 SAM-dependent methyltransferase DEGs were annotated with 6 up-regulation and 1 down-regulation. Typically, glycosylation usually occurs at the end of secondary metabolite biosynthesis to increase the stability, water solubility, and biological activity of the substance, and the glycosylation reaction is generally catalyzed by UDPGs in nature. As shown in Fig. 12, among the 13 UDPG DEGs annotated, 12 of them were up-regulated after the MeJA treatment. The heat expression maps of the three post-plant modification enzymes showed that their relevant DEG expression patterns were consistent, and the differential expression was more significant in MJ4 than in MJ24 group (Fig. 12).

Fig. 12. DEGs assigned to post-modification enzymes under MeJA treatment (the expression heatmaps from left to right were arranged in the following order: MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples).

11

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

studies on its functional genomics. In some public databases including NCBI, only a few complete or partial coding nucleotide sequences have been published. The EST or genome survey sequences (GSS) of Ci. aromaticum has not been found in the GenBank database. Therefore, the present study analyzed and reported the transcriptome sequencing date of Ci. aromaticum treated with MeJA for the first time, which could provide some useful information for further genetic study on this species. Transcriptome sequencing is an effective tool for large-scale gene identification and expression profiling (Bushmanova et al., 2019; Dong et al., 2018). In this study, the transcriptome of the MeJA-treated Ci. aromaticum was identified by using a high-throughput RNA-seq technology. A total of 176,091 unigenes were obtained after assembly. In addition, numerous DEGs involved in the biosynthetic and metabolism of flavor volatile compounds in Ci. aromaticum were also obtained, which will provide information in further understanding the biosynthetic pathways of the flavor volatiles and other secondary metabolites of Ci. aromaticum. The synthesis of JA is regulated by the feedback of JA and its synthetic pathway intermediates (Creelman and Mullet, 1997). Exogenous MeJA can induce the expression of JA biosynthetic pathway related genes and promote the synthesis of JA in plants (Men et al., 2013). Our transcriptome data confirmed the similar results. Endogenous MeJA was considered as a major regulator of plant JA signaling pathways. In general, the level of JA in plants is extremely low, and the inhibitory protein JAZ of JA signaling pathway binds to a series of transcription factors or signal transduction proteins such as MYC2 to inhibit the JA response. When stimulated by exogenous MeJA, JA levels in plants increased rapidly because JAZ was degraded by SCFCOI1 ubiquitin proteasome, resulting in the release of MYC2, thereby activating JAregulated gene expression (Wasternack and Song, 2017). Our results showed that after the MeJA treatment, 9 JAZ proteins and 3 MYC2 transcription factors were all up-regulated. These results indicated that the application of exogenous MeJA could mediate the JA biosynthetic and signaling pathways, thereby regulating a range of downstream genes. Terpenoids are the main aromatic components of Ci. aromaticum (Gao et al., 2017). Generally, in higher plants, terpenoids are mainly synthesized by the MVA/MEP pathway (Aharoni et al., 2004; Rohmer, 2003). In our transcriptome data, the differentially expressed key enzyme genes which were annotated to the MVA and MEP pathways, respectively, and their expression trends were consistent, indicating that the biosynthesis of terpenoids in Ci. aromaticum plants were accomplished through the MVA and MEP pathways. On the other hand, we found that after the MeJA treatment, the GPS enzyme gene for the production of the monoterpene substance synthesis substrate GPP and the FPS enzyme gene for the sesquiterpene and triterpenoid synthesis substrate FPP were significantly up-regulated by 2.57- and 2.14-folds, respectively, and the results were similar to some previous studies on Camellia sinensis (Shi et al., 2015) and Taxus chinensis (Li et al., 2012). In addition, a large number of monoterpene synthases and sesquiterpene synthases were found to be significantly up-regulated after the MeJA treatment. This observation agreed with our previous study, in which the content of monoterpenes and sesquiterpenoids in Ci. aromaticum after the MeJA treatment were significantly increased (Gao et al., 2017). Therefore, it is speculated that the content of these precursors and the expression of terpene synthases are the crucial factors for the release of aromatic volatiles from Ci. aromaticum. In addition to terpenes, some alcohols, esters, aldehydes, and ketones are also important aromatic volatiles of Ci. aromaticum. These substances are mainly synthesized by phenylpropanoid and fatty acid metabolic pathways (Dudareva et al., 2004; Schwab et al., 2008). In present study, there were 36 DEGs associated to phenylpropanoid metabolic pathway and 28 DEGs associated to the fatty acid metabolic pathway, all of which were up-regulated after the MeJA treatment, leading to the increase in the contents of aromatic chemicals of Ci.

Fig. 13. DEGs assigned to transcription factors under MeJA treatment (the expression heatmaps from left to right were arranged in the following order: MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4. The z-scores were applied to represent the fold changes among samples). A: Statistics of differentially expressed TFs; B: The expression heatmaps of TFs.

3.14. DEGs encoding transcription factors In plants, JAs regulate a variety of biological processes by stimulating transcription factors to regulate downstream pathway genes. In our study, many DEGs were predicted to be transcription factors. A total of 59 DEGs belonging to 22 transcription factor families was annotated, of which, Tify (16), AP2-EREBP (6), MYB (5), and bHLH (4) families account for a large proportion. These differentially expressed transcription factors are primarily involved in plant secondary metabolism or stress responses (Fig. 13). With aid of heat map analysis of the differentially expressed AP2-EREBP, MYB, and bHLH transcription factors that regulate secondary metabolism, it was found that all bHLH transcription factors were up-regulated after the MeJA treatment, while the AP2-EREBP and MYB transcription factors were in irregular changes. Most DEGs had the same expression pattern, and the MeJA-mediated induction was more obvious in MJ4 than in MJ24 group. 3.15. QRT-PCR validation of DEGs To further verify the gene expression profiles obtained by RNA-Seq and the reliability of important DEGs obtained from the assembled transcriptome, we selected eight key enzyme genes (i.e., HMGR, PMVK, DXS, DXR, ISPF, ISPG, GPS, and FPS) in terpenoid backbone biosynthetic pathways, three terpene synthase genes (i.e., LIS, (E)-β-ocimene synthase, β-caryophyllene synthase), one modifying CYPP450 enzyme, and four transcription factors (i.e., MYC2, AP2-EREBP, bHLH, and MYB) that regulate the synthesis of terpenoids for qRT-PCR. The results in Fig. 14 showed that although the expression level of these genes was different between the qRT-PCR results and the RNA-Seq results, which may be due to the differences in sensitivity, specificity, algorithms of qRT-PCR or sequencing techniques, the expression trend was almost the same. The results obtained from the qRT-PCR analysis confirmed that unigenes obtained from the assembled transcriptomes and the gene expression profiles from RNA-Seq data are reliable. 4. Discussion As a natural aromatic plant, Ci. aromaticum is an important resource of aromatic plant in the genus Chrysanthemum (He et al., 2016). However, previous researches on Ci. aromaticum mainly focused on the extraction of essential oils (Fan et al., 2018) and the analysis of aromatic constituents (Wang et al., 2018), and there have been only a few 12

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Fig. 14. qRT-PCR validation of sixteen DEGs involved in terpenoids biosynthesis treated with MeJA compared with the Ctrl. The Log2FC values obtained from RNASeq data are indicated on the top of each graph.

regulated after MeJA treatment, which may have played active roles in the regulation of JA-responsive transcriptional regulation of terpenoids. In contrast, we also found that 3 AP2/ERF and 3 MYB transcription factors were significantly down-regulated, similar to the results of transcriptome analyses found in Mentha canadensis L (Qi et al., 2018) and Taxus yunnanensis cultivar (He et al., 2018). Therefore, it is reasonable that those transcription factors in Ci. aromaticum have played negative regulatory roles in the synthesis of terpenoids.

aromaticum. Cytochrome CYP450 s, SAM-dependent methyltransferases, and UDPGs are mainly used for the modification of terpenoids in the later stage of their syntheses to enrich the species of terpenoids (Yamamura et al., 2017). We screened 16 CYP450 s, 6 SAM-dependent methyltransferases, and 12 UDPGs, and their expression patterns were consistent with the enzyme genes in the biosynthetic pathway of terpenoids. All of them were significantly up-regulated under the MeJA treatment. These selected post-plant modification enzymes could be used as candidate genes for genetic improvement for Ci. aromaticum compounds. In plants, a major regulatory mechanism for secondary metabolite production is to control biosynthetic genes by modifying the expression of transcription factors (Keasling, 2012; Yanagisawa, 1998). At present, several transcription factors including AP2/ERF (Chen et al., 2012; Yu et al., 2012), MYC2 (Shen et al., 2016), bHLH (Hartmann et al., 2005), MYB (Chezem and Clay, 2016) and WRKY(Ma et al., 2009; Xu et al., 2004), which regulate the biosynthesis of terpenes, have been identified from plants. Our transcriptome results exhibited that AP2-EREBP, MYC2, bHLH, and MYB transcription factors were significantly up-

5. Conclusions In this study, transcriptome sequencing analysis of MeJA-induced transcriptional changes in Ci. aromaticum was performed in order to identify candidate genes involved in the biosynthesis of secondary metabolites, especially genes for the terpenoids biosynthesis. A total of 176,091 unigenes were obtained, of which 69.94 % could be functionally annotated in at least one database. In the MJ4vsCtrl, MJ24vsCtrl, and MJ24vsMJ4 comparison groups, we obtained 750, 650, and 313 DEGs, respectively. Except the transcripts related to the JA biosynthesis and signal transduction, we found that, in the KEGG 13

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

pathway, most DEGs participated in metabolic pathways, particularly in the terpenoids, phenylpropanoids, fatty acid pathways, and are genes for secondary metabolism-related post-plant modification enzymes, and transcription factors. An increase in the level of transcription of key genes during the JA biosynthesis indicated that JA had an autoregulatory loop. Many DEGs associated with JA signaling were up-regulated, suggesting that exogenous MeJA may regulate the JA signaling pathway in Ci. aromaticum. The data also showed that the MeJA-responsive MEP and MEV pathways provided a 5C building block for the biosynthesis of different terpene metabolites. The up-regulated expression of PMVK, ISPD, ISPF, ISPG, GPS, and FPS enzyme genes increased the content of precursor substances for terpenoids biosynthesis. In addition, it was found that the expression of terpene synthetase genes under the MeJA induction was remarkably different, which substantially changed the content of terpenes in Ci. aromaticum. This provides the first piece of evidence confirming that MeJA induced the up regulation of the biosynthetic pathway of terpenoid backbone and terpenoid synthase genes in Ci. aromaticum. In addition, we also identified some key enzyme genes in the phenylpropanoid and fatty acid metabolic pathways, providing valuable resources for study of elucidation of the MeJA-regulated secondary metabolic processes. The differentially expressed transcription factors identified in the present study is valuable to understand their relative regulatory significance in different metabolic processes of plants and to conduct metabolic engineering studies. The present study may allow a comprehensive understanding of the response of Ci. aromaticum to the MeJA treatment and help explain the changes of terpenoids in Ci. aromaticum at the molecular level. The obtained data provided a large amount of genetic resources for terpenoids biosynthesis of Ci. aromaticum. This is the first report on the transcriptome of Ci. aromaticum, which provides useful genomic information for further exploitation.

Chezem, W.R., Clay, N.K., 2016. Regulation of plant secondary metabolism and associated specialized cell development by MYBs and bHLHs. Phytochemistry 131, 26–43. Cho, H.Y., Son, S.Y., Rhee, H.S., Yoon, S.Y.H., Lee-Parsons, C.W., Park, J.M., 2008. Synergistic effects of sequential treatment with methyl jasmonate, salicylic acid and yeast extract on benzophenanthridine alkaloid accumulation and protein expression in Eschscholtzia californica suspension cultures. J. Biotechnol. 135, 117–122. Cohen, S., Flescher, E., 2009. Methyl jasmonate: a plant stress hormone as an anti-cancer drug. Phytochemistry 70, 1600–1609. Creelman, R.A., Mullet, J.E., 1997. Biosynthesis and action of jasmonates in plants. Annu. Rev. Plant Biol. 48, 355–381. Davidson, N.M., Oshlack, A., 2014. Corset: enabling differential gene expression analysis for de novoassembled transcriptomes. Genome Biol. 15, 410. Dong, L., Huang, Z., Liu, D., Zhu, P., Lv, S., Li, N., Mao, H., 2018. Transcriptome analysis of chrysanthemum in responses to white rust. Sci. Hortic. 233, 421–430. Dudareva, N., Pichersky, E., Gershenzon, J., 2004. Biochemistry of plant volatiles. Plant Physiol. 135, 1893–1902. Fan, S., Chang, J., Zong, Y., Hu, G., Jia, J., 2018. GC-MS analysis of the composition of the essential oil from Dendranthema indicum Var. Aromaticum using three extraction methods and two columns. Molecules 23, 576. Gao, W., Wang, X., He, M., Zhou, Y., 2017. Effects of methyl jasmonate on the trichomes and secretions of Chrysanthemum indicum var. aromaticum in leaves. J. Beijing For. Univ. 39, 79–85. Gao, W., Wang, X., Purente, N., Muhammad, L., Zhou, Y., He, M., 2018. A 1-deoxy-dxylulose 5-phosphate reductoisomerase gene probably involved in the synthesis of terpenoids in Chrysanthemum indicum var. aromaticum. Can. J. Plant Sci. 98, 1254–1264. Ge, X., Wu, J., 2005. Induction and potentiation of diterpenoid tanshinone accumulation in Salvia miltiorrhiza hairy roots by β-aminobutyric acid. Appl. Microbiol. Biot. 68, 183–188. Hartmann, U., Sagasser, M., Mehrtens, F., Stracke, R., Weisshaar, B., 2005. Differential combinatorial interactions of cis-acting elements recognized by R2R3-MYB, BZIP, and BHLH factors control light-responsive and tissue-specific activation of phenylpropanoid biosynthesis genes. Plant Mol. Biol. 57, 155–171. He, M., Gao, W., Gao, Y., Liu, Y., Yang, X., Jiao, H., Zhou, Y., 2016. Polyploidy induced by colchicine in Dendranthema indicum var. aromaticum, a scented chrysanthemum. Eur. J. Hortic. Sci. 81, 219–226. He, C.T., Li, Z.L., Zhou, Q., Shen, C., Huang, Y.Y., Mubeen, S., Yang, J.Z., Yuan, J.G., Yang, Z.Y., 2018. Transcriptome profiling reveals specific patterns of paclitaxel synthesis in a new Taxus yunnanensis cultivar. Plant Physiol. Biochem. 122, 10–18. Jin, L., Yang, Y., Gao, W., Gong, M., Wang, J., Anderson, N., He, M., 2017. Establishment of callus induction and cell suspension cultures of Dendrathema indicum var. aromaticum a scented chrysanthemum. J. Plant Stud. 6, 38–44. Keasling, J.D., 2012. Synthetic biology and the development of tools for metabolic engineering. Metab. Eng. 14, 189–195. Li, S.T., Zhang, P., Zhang, M., Fu, C.H., Zhao, C.F., Dong, Y.S., Guo, A.Y., Yu, L.J., 2012. Transcriptional profile of Taxus chinensis cells in response to methyl jasmonate. BMC Genomics 13, 295. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408. Love, M.I., Huber, W., Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. Lubbe, A., Verpoorte, R., 2011. Cultivation of medicinal and aromatic plants for specialty industrial materials. Ind. Crop. Prod. 34, 785–801. Ma, D., Pu, G., Lei, C., Ma, L., Wang, H., Guo, Y., Chen, J., Du, Z., Wang, H., Li, G., 2009. Isolation and characterization of AaWRKY1, an Artemisia annua transcription factor that regulates the amorpha-4, 11-diene synthase gene, a key gene of artemisinin biosynthesis. Plant Cell Physiol. 50, 2146–2161. Men, L., Yan, S., Liu, G., 2013. De novo characterization of Larix gmelinii (Rupr.) Rupr. transcriptome and analysis of its gene expression induced by jasmonates. BMC Genomics 14, 548. Nkembo, K.M., Lee, J.B., Hayashi, T., 2005. Selective enhancement of scopadulcic acid B production in the cultured tissues of Scoparia dulcis by methyl jasmonate. Chem. Pharm. Bull. 53, 780–782. Qi, X., Fang, H., Yu, X., Xu, D., Li, L., Liang, C., Lu, H., Li, W., Chen, Y., Chen, Z., 2018. Transcriptome analysis of ja signal transduction, transcription factors, and monoterpene biosynthesis pathway in response to methyl jasmonate elicitation in Mentha canadensis L. Int. J. Mol. Sci. 19, 2364. Qihong, L., Shufan, Z., 1983. A new variety of Dendranthema Gaertn. from shennongjia of Hubei. J. Wuhan Bot. Res. 2, 237–238. Rai, A., Singh, R., Shirke, P.A., Tripathi, R.D., Trivedi, P.K., Chakrabarty, D., 2015. Expression of rice CYP450-like gene (Os08g01480) in Arabidopsis modulates regulatory network leading to heavy metal and other abiotic stress tolerance. PLoS One 10, e0138574. Ramya, H., Palanimuthu, V., Rachna, S., 2013. An introduction to patchouli (Pogostemon cablin Benth.)–a medicinal and aromatic plant: It’s importance to mankind. Agric Eng Int: CIGR J. 15, 243–250. Rohmer, M., 2003. Mevalonate-independent methylerythritol phosphate pathway for isoprenoid biosynthesis. Elucidation and distribution. Pure Appl. Chem. 75, 375–388. Schwab, W., Davidovich‐Rikanati, R., Lewinsohn, E., 2008. Biosynthesis of plant‐derived flavor compounds. Plant J. 54, 712–732. Shen, Q., Lu, X., Yan, T., Fu, X., Lv, Z., Zhang, F., Pan, Q., Wang, G., Sun, X., Tang, K., 2016. The jasmonate‐responsive AaMYC2 transcription factor positively regulates artemisinin biosynthesis in Artemisia annua. New Phytol. 210, 1269–1281. Shi, J., Ma, C., Qi, D., Lv, H., Yang, T., Peng, Q., Chen, Z., Lin, Z., 2015. Transcriptional responses and flavor volatiles biosynthesis in methyl jasmonate-treated tea leaves.

Contributions Yunwei Zhou, and Miao He designed the experiments; Wenjie Gao analyzed the data and prepared the manuscript; Qingran Meng helped carried out the experiment; and Hong Luo and Feng Chen reviewed and revised the manuscript. All authors read and approved the final manuscript version. Funding This research was supported by the National Natural Science Foundation of China (31870687) and Fundamental research funds for the central university (2572017EA04). Declaration of Competing Interest The authors have no conflicts of interest to declare. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.indcrop.2020.112254. References Aharoni, A., Giri, A.P., Verstappen, F.W., Bertea, C.M., Sevenier, R., Sun, Z., Jongsma, M.A., Schwab, W., Bouwmeester, H.J., 2004. Gain and loss of fruit flavor compounds produced by wild and cultivated strawberry species. Plant Cell 16, 3110–3131. Bushmanova, E., Antipov, D., Lapidus, A., Przhibelskiy, A.D., 2019. rnaSPAdes: a de novo transcriptome assembler and its application to RNA-Seq data. GigaScience 8 giz100. Chen, Y.Y., Wang, L.F., Dai, L.J., Yang, S.G., Tian, W.M., 2012. Characterization of HbEREBP1, a wound-responsive transcription factor gene in laticifers of Hevea brasiliensis Muell. Arg. Mol. Biol. Rep 39, 3713–3719. Cheng, A.X., Lou, Y.G., Mao, Y.B., Lu, S., Wang, L.J., Chen, X.Y., 2007. Plant terpenoids: biosynthesis and ecological functions. J. Integr. Plant Biol. 49, 179–186.

14

Industrial Crops & Products 147 (2020) 112254

W. Gao, et al.

Yanagisawa, S., 1998. Transcription factors in plants: physiological functions and regulation of expression. J. Plant Res. 111, 363–371. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14. Yu, Z.X., Li, J.X., Yang, C.Q., Hu, W.L., Wang, L.J., Chen, X.Y., 2012. The jasmonateresponsive AP2/ERF transcription factors AaERF1 and AaERF2 positively regulate artemisinin biosynthesis in Artemisia annua L. Mol. Plant 5, 353–365. Zang, Y., Zheng, W., He, Y., Hong, S.B., Zhu, Z., 2015. Global analysis of transcriptional response of Chinese cabbage to methyl jasmonate reveals JA signaling on enhancement of secondary metabolism pathways. Sci. Hortic. 189, 159–167.

BMC Plant Biol. 15, 233. Wang, Y., Li, X., Jiang, Q., Sun, H., Jiang, J., Chen, S., Guan, Z., Fang, W., Chen, F., 2018. GC-MS analysis of the volatile constituents in the leaves of 14 compositae plants. Molecules 23, 166. Wasternack, C., Song, S., 2017. Jasmonates: biosynthesis, metabolism, and signaling by proteins activating and repressing transcription. J. Exp. Bot. 68, 1303–1321. Xu, Y.H., Wang, J.W., Wang, S., Wang, J.Y., Chen, X.Y., 2004. Characterization of GaWRKY1, a cotton transcription factor that regulates the sesquiterpene synthase gene (+)-δ-cadinene synthase-A. Plant Physiol. 135, 507–515. Yamamura, Y., Kurosaki, F., Lee, J.B., 2017. Elucidation of terpenoid metabolism in Scoparia dulcis by RNA-seq analysis. Sci. Rep. 7, 43311.

15