Transcriptome Analysis of JA Signal Transduction, Transcription Factors, and Monoterpene Biosynthesis Pathway in Response to Methyl Jasmonate Elicitation in Mentha canadensis L.

Mentha canadensis L. has important economic value for its abundance in essential oils. Menthol is the main component of M. canadensis essential oils, which is certainly the best-known monoterpene for its simple structure and wide applications. However, the regulation of menthol biosynthesis remains elusive in M. canadensis. In this study, transcriptome sequencing of M. canadensis with MeJA treatment was applied to illustrate the transcriptional regulation of plant secondary metabolites, especially menthol biosynthesis. Six sequencing libraries were constructed including three replicates for both control check (CK) and methyl jasmonate (MeJA) treatment and at least 8 Gb clean bases was produced for each library. After assembly, a total of 81,843 unigenes were obtained with an average length of 724 bp. Functional annotation indicated that 64.55% of unigenes could be annotated in at least one database. Additionally, 4430 differentially expressed genes (DEGs) with 2383 up-regulated and 2047 down-regulated transcripts were identified under MeJA treatment. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment indicated that “Monoterpenoid biosynthesis” was one of the most significantly enriched pathways in metabolism. Subsequently, DEGs involved in JA signal transduction, transcription factors, and monoterpene biosynthesis were analyzed. 9 orthologous genes involved in menthol biosynthesis were also identified. This is the first report of a transcriptome study of M. canadensis and will facilitate the studies of monoterpene biosynthesis in the genus Mentha.


Introduction
The genus Mentha has important economic value for its abundance in essential oils that are widely used in the flavor, fragrance, and aromatherapy industries [1]. Monoterpenes are the major constituents of essential oils, which represent a large and diverse class of volatile C10 isoprenoids. The biosynthesis of monoterpenes originated in the upstream methylerythritol phosphate (MEP) pathway, then

Transcriptome Sequencing and De Novo Assembly of M. canadensis
Transcriptome sequencing was performed for M. canadensis under control check (CK) and MeJA treatment using Illumina HiSeq™ 4000. Six sequencing libraries were constructed including three replicates for both CK and MeJA treatment. As a result, each library produced at least 54 Mb clean reads and 8 Gb clean bases. The percentage of clean reads was more than 98% and Q20 percentage was at least 98.69% for the six libraries (Table 1). Then, clean reads of the six libraries were de novo assembled using Trinity program and finally 81,843 unigenes were generated with a GC percentage of 43.92%. The average length and N50 length for the assembled unigenes were 724 and 1126 bp, respectively ( Table 2). The length distribution of the unigenes is indicated in Supplementary Figure S1 and 70.63% of all unigenes show lengths longer than 300 bp.

Functional Annotation of Unigenes
For functional annotation analysis, all the assembled unigenes were searched against the public databases using BLAST. The results show that 52,700 (64.39%) unigenes could be annotated in Nr (non-redundant protein) database, 34,565 (42.23%) in Swissprot, 29,536 (36.09%) in KOG (Clusters of eukaryotic Orthologous Group), and 19,013 (23.23%) in KEGG (Kyoto Encyclopedia of Genes and Genomes). Taken together, there were 52,826 (64.55%) unigenes that could be annotated in at least one database (Table 3). GO (Gene Ontology), KOG and KEGG assignments were used to classify the function of the assembled unigenes of M. canadensis. For GO annotation, 4396, 4422 and 2566 unigenes were assigned to "biological process" category, "molecular function" category and "cellular component" category, respectively (Supplementary Figure S2). In the "molecular function" category, genes assigned to "binding" (2260) and "catalytic activity" (3326) accounted for the vast majority of this category. In the "biological process" category, "cellular process" (2891), "metabolic processes" (3447) and "single-organism process" (2344) were the main subcategories. For KOG classification, 29,536 unigenes were categorized into 25 KOG functional groups. Among them, "General functional prediction only" (9523) was the largest group, followed by "Signal transduction mechanisms" (5796) and "Posttranslational modification, protein turnover, chaperones" (4901) (Supplementary Figure S3). To further understand the biological functions of M. canadensis unigenes, KEGG pathway analysis was performed. As a result, a total of 10,365 unigenes could be mapped to 133 metabolic pathways. "Metabolism" (5104) was the largest category (Supplementary Figure S4).

Identification and Analysis of Differentially Expressed Genes (DEGs)
To identify differentially expressed genes under MeJA treatment, reads were mapped to the unigenes and reads per kb per million reads (RPKM) was used to measure the transcriptional levels. The Pearson correlation analysis indicated that there were high correlations between three biological replicates of both CK and MeJA treatment (Supplementary Figure S5). Then using a two-fold change of RPKM with False Discovery Rate (FDR) <0.05 as the cut-off criteria, 4430 differentially expressed genes (DEGs) were identified between CK and MeJA-treated samples in M. canadensis. Among them, 2383 DEGs were up-regulated and 2047 DEGs were down-regulated under MeJA treatment ( Figure 1). To further explore the possible roles of the DEGs, GO and KEGG enrichment were conducted. For GO enrichment, 88, 27 and 6 GO terms were significantly enriched in "biological process", "molecular function" and "cellular component" categories, respectively (Supplementary Table S1). "Terpene synthase activity" and "oxidoreductase activity" were the most significantly enriched terms in "molecular function". In "biological process" category, biosynthetic processes involved in lipid, small molecule, dicarboxylic acid, phospholipid, terpene and sesquiterpene were the most significantly enriched terms. For KEGG pathway enrichment, 25 pathways were significantly enriched under MeJA treatment ( Figure 2 and Supplementary Table S2). Of these, pathways involved in secondary metabolism accounted for a large part, including biosynthetic pathways of phenylpropanoid, terpenoid, alkaloid, flavonoid, glucosinolate and so on.

Identification and Expression Verification of JA Signal Pathway Genes under MeJA Treatment
In Arabidopsis, the JA signal is perceived and transduced via the SCF COI1 -JAZ co-receptor complex. The SCF COI1 complex includes COI1, ASK2, CULLIN1, Rbx, and E2. The JAZ repressors recruit the co-repressor complex consisting of Novel Interactor of JAZ (NINJA), TOPLESS (TPL), and histone deacetylase (HDA), interact with and repress the transcription activator bHLHzip transcription factor MYC2 (MYC2). The perception of JA-Ile by SCF COI1 complex leads to degradation of JAZs via the 26S proteasome, then activate the downstream transcription factors of JA responses [17,31]. In this study, 24 DEGs associated with JA signal transduction were identified in M. canadensis. Of these, 23 DEGs showed up-regulation under MeJA treatment, which encode putative JAZ (9), MYC2 (7), NINJA (6), and S-phase kinase-associated protein 1 (SKP1) (1) (Figure 3). 12 of the 24 DEGs were selected for quantitative real-time PCR (qRT-PCR) to verify the reliability of the RNA-Seq data. The results showed that all the expression patterns of qRT-PCR were quite a good match with the RPKM results of RNA-Seq data ( Figure 4). These validation results of qRT-PCR indicated that the RNA-Seq data were quite reliable.

Identification and Analysis of Differentially Expressed TFs under MeJA Treatment
To further understand the transcriptional regulation of MeJA-treated M. canadensis, differentially expressed TFs were identified between CK and MeJA-treated samples. As a result, 1535 TFs belonged to 55 families were identified by annotating in PlantTFDB database. bHLH (basic/helix-loop-helix), ERF (ethylene-responsive factor), MYB_related (myeloblastosis DNA-binding related protein), C2H2 (C2H2 zinc-finger protein), MYB (myeloblastosis DNA-binding protein), WRKY (WRKY-type DNA binding protein), C3H (Cys3His zinc finger protein), NAC (NAC domain protein), bZIP (basic region/leucine zipper motif), and GRAS (Gibberellin-insensitive (GAI), Repressor of GA (RGA), and SCARECROW (SCR) protein) families were the top 10 largest TF families (Supplementary Table S3). Of the 1535 TFs, 157 belonged to 32 families were differentially expressed. WRKY, bHLH, ERF, MYB, MYB_related, NAC, bZIP, C2H2, HD-ZIP (Homeodomain-leucine zipper protein), and Trihelix families were the top 10 families that had the largest number of DEGs. Each of these 10 families had more than 5 DEGs and the WRKY and bHLH family were the largest differentially expressed families that both contained 21 DEGs ( Figure 5 and Supplementary Table S3). Comparing the expression levels of these DEGs, 102 of the 157 DEGs were up-regulated and the other 55 were down-regulated. Considering the family distribution of the up-regulated and down-regulated DEGs, 19 families showed that the up-regulated unigenes were more abundant than the down-regulated ones, especially in the WRKY family, 19 unigenes were up-regulated and only 2 were down-regulated. On the contrary, there were 11 families that contained more down-regulated unigenes than up-regulated and 2 families contained the same number of up-regulated unigenes and down-regulated ( Figure 5).

Identification and Expression Verification of Monoterpenoids and Menthol Biosynthetic Genes under MeJA Treatment
In M. canadensis, monoterpenes are the major constituents of essential oils. In this study, the enrichment results of the KEGG database indicated that "Monoterpenoid biosynthesis" was one of the most significantly enriched pathways under MeJA treatment. 20 DEGs associated with monoterpenoid biosynthesis were identified, which consisted of genes encoding 10 terpene synthases (TPSs), 1 (−)-isopiperitenone reductase (IPR), 1 menthol dehydrogenase (MR), and 8 neomenthol dehydrogenases (NMRs). TPS is a multi-gene family that is widespread in plants, whose diversity and substrate complexity lead to the wide variety of terpenoids [32]. In this study, 10 differentially expressed TPSs were identified under MeJA treatment, 8 of which were up-regulated and the other 2 were down-regulated. The remaining 10 DEGs in "Monoterpenoid biosynthesis" were menthol/neomenthol biosynthesis related genes. Interestingly, the expression levels of the 10 DEGs were almost all up-regulated except for an NMR (Unigene0051701) under MeJA treatment ( Figure 6). Menthol is the main constituent of monoterpenoids in M. canadensis, which is synthesized by a series of enzymatic reactions in the peltate glandular trichomes on the aerial surfaces of the plant. The biosynthetic pathway and enzyme catalysis mechanism of menthol has been extensively studied in peppermint and spearmint, two other species of genus Mentha [2][3][4][5][6]. Using the reference genes of peppermint and spearmint as queries, 9 orthologous genes were identified in M. canadensis (Supplementary Figure S6), including the previously reported GPPS-l, GPPS-s, and LS [15]. The other 6 genes that catalyze the biosynthesis of (−)-Menthol from (−)-Limonene were identified in M. canadensis for the first time. Considering the transcriptional levels of the 9 genes, although their RPKMs showed a certain level of increase after MeJA treatment, they did not reach significant levels above controls. Then, qRT-PCR was used to verify the expression of the 9 menthol biosynthetic genes (Figure 7). The results indicated that the expression levels of most genes were not significantly different between MeJA treatment and CK groups, which was consistent with the RNA-seq results. Among these genes, only Unigene0030907 (MFS) showed significant down-regulated expression after MeJA treatment. Unigene0038587 (GPPS-s) and Unigene0033747 (iPD) showed a certain degree of up-regulation in both qRT-PCR and RNA-Seq data.

Discussion
Transcriptome sequencing is an effective tool for large-scope gene identification and expression profiling analysis. In this study, high-throughput RNA-seq was applied to characterize the transcriptomes of M. canadensis treated with MeJA. Six sequencing libraries were constructed including three replicates for both CK and MeJA treatment and at least 8 Gb clean data was produced for each library. After assembly, a total of 81,843 unigenes were obtained with an average length of 724 bp. Functional annotation indicated that 64.55% of unigenes could be annotated in at least one database. This is the first reported transcriptome study of M. canadensis and will provide useful information for further study of this species.
Endogenous MeJA is believed to be a primary regulator of the JA signal pathway in plants. In Arabidopsis, the JA signal is perceived and transduced via the SCF COI1 -JAZ co-receptor complex. In the resting state, JAZ proteins interact with and repress downstream MYC2 to suppress of JA responses. In response to JA signal, JAZs are degraded by SCF COI1 -ubiquitin-proteasome, then MYC2 is released from the repressor complex to regulate gene expression [17,31]. In this transcriptome dataset, 24 DEGs associated with JA signal transduction were identified in M. canadensis. Of these, 23 DEGs showed up-regulation under MeJA treatment, encoding putative JAZ (9), MYC2 (7), TPL (1), NINJA (6), and SKP1 (1). qRT-PCR was also used to verify the reliability of the RNA-Seq data. These results suggested that application of exogenous MeJA may regulate the JA signaling pathway in M. canadensis.
TFs are sequence specific DNA-binding proteins that interact with the regulatory regions of the target genes and regulate their expression. Several transcriptome profiling studies have indicated that MeJA treatment triggers an extensive transcriptional reprogramming of metabolism [33][34][35]. A lot of TFs are elicited by MeJA and participate in regulation of specific metabolic processes [36]. Here, 157 differentially expressed TFs (102 up-regulated and 55 down-regulated) belonged to 32 families were identified by annotating in PlantTFDB database. It is worth noting that the WRKY family was the largest differentially expressed family containing 21 DEGs. It is also reported that more than 25% of WRKYs are induced in response to jasmonate in A. thaliana and Catharanthus roseus [37]. Of the 21 differentially expressed WRKYs, the number of up-regulated genes was much higher than that of down-regulated ones (19 vs. 2). Similar observations have been found in other plants. For example, 16 up-regulated and 3 down-regulated WRKY TFs were identified under MeJA treatment in A. annua [27]. In Lycoris aurea, 32 differentially expressed WRKY TFs were identified, of which 26 were up-regulated and only 6 were down-regulated [38]. The WRKY family is an important class of JA-responsive TFs that regulate plant secondary metabolism [39]. In A. annua, a WRKY TF AaWRKY1 was strongly induced by MeJA and regulated the expression of ADS, a key gene of artemisinin biosynthesis [21]. In Panax quinquefolius, a MeJA-responsive WRKY TF PqWRKY1 was isolated, which could positively regulate the biosynthesis of triterpene ginsenoside [40]. The large number of up-regulated WRKY TFs in M. canadensis suggests that they may function in JA-responsive transcriptional regulation of secondary metabolism.
Monoterpenes are the major constituents of M. canadensis essential oils. In this study, the KEGG pathway enrichment of DEGs indicated that "Monoterpenoid biosynthesis" was the most significantly enriched pathway under the treatment of MeJA. 20 DEGs associated with monoterpenoid biosynthesis were identified, which consisted of genes encoding TPS (10), IPR (1), MR (1), and NMR (8). Of these, genes encoding IPR, MR, and NMR were involved in menthol and neomenthol biosynthesis and interestingly, the expression levels of the 10 DEGs were almost all up-regulated under MeJA treatment. These results suggest that the biosynthesis of monoterpenes including menthol and neomenthol might be elicited by the MeJA signal in M. canadensis. The biosynthesis pathway of menthol has been extensively studied in peppermint and spearmint [2][3][4][5][6]. In this study, 9 orthologous genes were identified in M. canadensis including 6 new reported genes. However, expression analysis of these 9 genes in CK and MeJA-treated samples indicated that most of the gene expression levels did not change significantly after treatment. qRT-PCR results indicated that GPPS-s and iPD showed a certain degree of up-regulation and MFS was down-regulated. We speculate that only part of the menthol biosynthetic genes responds to the JA signal in M. canadensis. The response mechanism of menthol biosynthetic genes to JA signaling needs further study in M. canadensis, especially the possible DEGs including GPPS-s, iPD, and MFS. Another explanation might be that different secondary metabolic pathways are controlled by different regulatory modules. For example, in Medicago truncatula, genes involved in phenylpropanoid biosynthesis were transiently induced after application of MeJA at low concentrations (0.5-5 µmol/L), but triterpene biosynthetic genes were induced after 12-24 h at high concentrations (5-500 µmol/L) [41]. Many M. canadensis menthol biosynthetic genes might not respond to the MeJA signal under our treatment conditions (200 µmol/L, 24 h).

RNA Extraction, cDNA Library Construction and Illumina Sequencing
Total RNA of the M. canadensis samples was extracted using RNAiso Plus (Takara, Dalian, China) according to the manufacturer's instructions. The quality and concentration of RNAs were measured using a ND-2000 UV spectrophotometer (Nanodrop Technologies, Wilmington, DE, USA) and Agilent 2100 (Agilent, Santa Clara, CA, USA). cDNA library construction and Illumina sequencing were performed by Gene Denovo Biotechnology Co. (Guangzhou, China). Briefly, mRNA was enriched from total RNA by Oligo(dT) beads and then fragmented into short fragments. The fragmented mRNA was reverse transcribed into cDNA with random primers and second-strand cDNA was subsequently synthesized. Then the cDNA fragments were purified, end repaired, poly(A) added, and ligated to sequencing adapters. After size selection by agarose gel electrophoresis, the ligation products were PCR amplified and sequenced using Illumina HiSeq TM 4000 (Illumina, San Diego, USA).

Transcriptome Assembly and Annotation
The raw reads obtained from RNA sequencing were filtered by trimming adapters and removing low quality reads (reads containing more than 10% of unknown nucleotides (N) or 40% of low quality (Q-value ≤ 10) bases) to obtain the high-quality clean reads. All data generated in this study have been deposited in the National Center for Biotechnology Information (NCBI) and can be accessed in the Short Read Archive (SRA) Sequence Database under accession number SRP132644. Then, Trinity program was used to perform de novo assembly of the clean data [51]. After assembly, unigenes were obtained and functional annotation was carried out using BLASTx program (http://www.ncbi. nlm.nih.gov/BLAST/) with an E-value threshold of 1 × 10 −5 . Public databases including the NCBI non-redundant protein database (Nr) (http://www.ncbi.nlm.nih.gov), the SwissProt database (http: //www.expasy.ch/sprot), the KEGG (http://www.genome.jp/kegg), and the COG/KOG database (http://www.ncbi.nlm.nih.gov/COG) were used to annotate the M. canadensis unigenes. TFs were identified by aligning unigenes to the PlantTFDB database [52].

Differentially Expressed Genes Analysis
The abundances of unigenes were calculated and normalized to RPKM [53]. To identify DEGs between MeJA-treated samples and controls, the edgeR package (http://www.r-project.org/) was used. Genes with a fold change ≥2 and a FDR <0.05 in a comparison were identified as significant DEGs. Pathway enrichment analysis was conducted by comparing gene numbers of each pathway in DEGs to genome background. The significance test of enriched pathways was determined by calculating the p-value and FDR correction. Heatmaps of DEGs were generated using HemI [54].

Quantitative Real-Time PCR
qRT-PCR was conducted to verify the expression of selected genes in M. canadensis. Goldenstar™ RT6 cDNA Synthesis Kit (TsingKe Biotech, Nanjing, China) was used with 1 µg of total RNA to synthesize first-strand cDNA. The qRT-PCR reactions were carried out using the qTOWER2.2 Real-Time PCR Systems (Analytik, Jena, Germany). 2 × T5 Fast qPCR Mix Kit (SYBR Green I) (TsingKe Biotech, Nanjing, China) was used to prepare qRT-PCR reactions with 2 µL of diluted cDNA as a template. The reaction systems and steps were performed according to the manufacturer's instructions. The M. canadensis actin gene was used as a control to normalize the relative expression levels of target genes. All results were representative of three independent experiments. Primers used for qRT-PCR were listed in Supplementary Table S4.

Conclusions
In this study, high-throughput RNA-seq was applied to characterize the transcriptomes of M. canadensis treated with MeJA. A total of 81,843 unigenes were obtained and 64.55% of which could be functionally annotated in at least one database. Additionally, 4430 DEGs with 2383 up-regulated and 2047 down-regulated transcripts were identified under MeJA treatment. A lot of unigenes associated with JA signal transduction were up-regulated, which suggested that application of exogenous MeJA may regulate the JA signaling pathway in M. canadensis. KEGG enrichment indicated that "Monoterpenoid biosynthesis" was one of the most significantly enriched pathways in metabolism. 9 orthologous genes involved in menthol biosynthesis were identified in M. canadensis and the response mechanism of menthol biosynthetic genes to JA signaling needs further study, especially the possible DEGs including GPPS-s, iPD, and MFS. 157 differentially expressed TFs belonged to 32 families were identified and the WRKY family was the largest differentially expressed family. The number of up-regulated WRKY TFs was much higher than that of down-regulated ones (19 vs. 2). The large number of up-regulated WRKY TFs in M. canadensis suggests that they may play important roles in JA-responsive transcriptional regulation of secondary metabolism. However, the regulation mechanism of JA signaling on development and metabolism of M. canadensis still requires further characterization and putative DEGs identified in this study might be important targets. This is the first reported transcriptome study of M. canadensis and will provide useful information for further study of this species.