Differential Gene Expression Analysis in Polygonum minus Leaf upon 24 h of Methyl Jasmonate Elicitation

Polygonum minus is an herbal plant that grows in Southeast Asian countries and traditionally used as medicine. This plant produces diverse secondary metabolites such as phenolic compounds and their derivatives, which are known to have roles in plant abiotic and biotic stress responses. Methyl jasmonate (MeJA) is a plant signaling molecule that triggers transcriptional reprogramming in secondary metabolism and activation of defense responses against many biotic and abiotic stresses. However, the effect of MeJA elicitation on the genome-wide expression profile in the leaf tissue of P. minus has not been well-studied due to the limited genetic information. Hence, we performed Illumina paired-end RNA-seq for de novo reconstruction of P. minus leaf transcriptome to identify differentially expressed genes (DEGs) in response to MeJA elicitation. A total of 182,111 unique transcripts (UTs) were obtained by de novo assembly of 191.57 million paired-end clean reads using Trinity analysis pipeline. A total of 2374 UTs were identified to be significantly up-/down-regulated 24 h after MeJA treatment. These UTs comprising many genes related to plant secondary metabolite biosynthesis, defense and stress responses. To validate our sequencing results, we analyzed the expression of 21 selected DEGs by quantitative real-time PCR and found a good correlation between the two analyses. The single time-point analysis in this work not only provides a useful genomic resource for P. minus but also gives insights on molecular mechanisms of stress responses in P. minus.


INTRODUCTION
Due to the sessile lifestyle, plants have to cope with different changes of biotic/abiotic factors in their surrounding environments to survive in nature. However, this defense response is costly and often setback by the repression of growth (Yang et al., 2012b). In response to stresses, plants convert resource allocation from growth to the biosynthesis of defensive compounds which is energetically demanding (Attaran et al., 2014). Plants produce a wide range of secondary metabolites for coordinated plant response to various internal and external cues (Davies and Schwinn, 2003). Plant secondary metabolites comprise a large group of lipophilic volatile organic compounds (VOCs) with high vapor pressure and low molecular weight. VOCs serve as signal molecules to mediate plant communication with their environment for herbivore deterrence, attraction of pollinators, seed dispersers, and protection against different stresses (Dudareva et al., 2013).
Beside their main role in defending plants from environmental stresses, investigations on the biosynthesis of VOCs have also attracted much interest for their applications as pharmaceuticals, pesticides, fragrance, and flavoring (Crozier et al., 2008). However, due to low abundance of these secondary metabolites (<1% dry weight) and high industrial demands, some plants in the wild which produce valuable compounds are overharvested and exposed to endangerment (Yang et al., 2014). Therefore, a comprehensive knowledge of metabolic pathways and regulatory mechanisms involved in the biosynthesis of VOCs is necessary to understand their role in defense responses and to enhance their yields for various applications.
Increased accumulation of volatile secondary metabolites in plants leads to enhanced stress tolerance and the activation of defense mechanisms (Aranega-Bou et al., 2014). A multilevel network of biosynthesis pathways and regulation of secondary metabolites converges on complex hormone signaling cascades in which jasmonates (JAs) play a crucial role in transcriptional control of plant defense gene expression and metabolism (Attaran et al., 2014).
JAs consisting jasmonic acid (JA) and its cyclopentanon derivatives such as its volatile methyl ester (methyl jasmonate, MeJA), are plant-specific endogenous signaling phytohormones that have long been observed to be potent regulators of elicitor signals for the biosynthesis of plant secondary metabolites (Zhao et al., 2005). This endogenous hormone regulates a variety of physiological plant processes such as root elongation, vegetative growth, production of viable pollen, senescence, cell cycle regulation, fruit ripening, production of specialized metabolites, plant response to wounding and defenses against pathogens, insects, and notorious pests (Farmer et al., 2003). It is wellstudied that elicitation is a common method of inducing the accumulation of plant defensive secondary metabolites in whole plant or cell culture (Memelink et al., 2001). Previous genomewide transcript profiling studies revealed that the addition of JA elicitors such as MeJA induces an extensive transcriptional reprogramming (De Geyter et al., 2012) leading to activation of several metabolic pathways including terpenoids (Misra et al., 2014), phenylpropanoids (Cocetta et al., 2015), and alkaloids (Kang et al., 2004). In addition to defense induction, the exogenous MeJA elicitation affects other plant functions including ontogeny, vegetative growth, and photosynthesis (Moreira et al., 2009). Previous studies showed that exogenous application of MeJA induces defense responses and reduces growth in several species (Heijari et al., 2005;Nabity et al., 2012;Yang et al., 2012b;Noir et al., 2013;Attaran et al., 2014).
Polygonum minus Huds (syn. Persicaria minor), commonly known as "kesum" in Malaysia, is an aromatic herbal plant in the Polygonaceae family which is widely used as a flavoring ingredient in local cuisine "asam laksa" and traditionally used to treat digestive disorders (Christapher et al., 2015). P. minus leaf is rich in VOCs comprising terpenoids (Baharum et al., 2010), aldehydes (Yaacob, 1990), and phenolic compounds (Maizura et al., 2011). P. minus leaf extracts possess high free radical scavenging activity and reducing power due to the high content of total phenolic compounds (Qader et al., 2011). This plant was shown to have antiviral, antimicrobial and anticancer activity (Vikram et al., 2014). Hence, it has a great potential for various applications in the pharmaceutical and perfumery industries for its essential oil.
Previously, an elicitation experiment using MeJA was performed on P. minus leaves through cDNA-amplified fragment length polymorphism (AFLP) transcript profiling approach with reported effects on some of the genes involved in secondary metabolite biosynthesis (Ee et al., 2013). However, detailed understanding of the biological mechanism of MeJA-mediated production of secondary metabolites and whole transcriptome changes was incomplete due to limitations in the gel-based selection of cDNA-AFLP transcripts for sequencing.
To date, there are only 3352 expressed sequence tags (ESTs) of P. minus available in the NCBI database from previous study (Roslan et al., 2012). In the current study, we aim to identify differentially expressed genes (DEGs) in the MeJAelicited P. minus leaf at 24 h compared with mock treatment. To gain an in-depth knowledge of gene expression changes in MeJA-elicited leaves, RNA-seq was employed to greatly expand the genomic resources (Loke et al., 2016). Over 192 million short pair-end reads were generated by Illumina HiSeq TM 2000 platform and used for de novo transcriptome assembly through Trinity analysis pipeline. The assembled UTs were analyzed to identify DEGs for functional annotation and downstream analyses. This is the first report on genome-wide transcriptional response in P. minus leaf elicited by MeJA with a focus on the biosynthesis of secondary metabolites. Furthermore, understanding the molecular mechanisms underlying MeJA elicitation will assist in genetic engineering of targeted secondary metabolites in P. minus.

Plant Materials and MeJA Treatment
P. minus stem cuttings were collected from Genting Highland (3 • 25 ′ 42.18 ′′ N, 101 • 47 ′ 21.45 ′′ E) and propagated in controlled environment chambers (A1000, Conviron, Canada) at 22/16 • C day/night temperatures under 12 h light/dark photoperiod with light intensity of 170 ± 20 µmol m −2 s −1 at ∼75% RH. After 45 days (Khairudin et al., 2014), an aqueous solution of 150 µM MeJA (Sigma-Aldrich) (Ismail et al., 2011) and 0.01% (v/v) Tween 20 was sprayed on all plant leaves to the point of runoff. The control (mock-treated) plant leaves were sprayed with only 0.01% (v/v) Tween 20. After elicitation, treated and control plants were placed separately in different growth chambers. Expanded young leaves from apical parts of plants (Ee et al., 2013) were harvested at 24 h after treatment and immediately frozen in liquid nitrogen before storage at −80 • C. Two biological replicates from independent control and treated plants were prepared for RNA sequencing. qRT-PCR was performed using three biological and three technical replicates.
to high phenolic compounds and polysaccharides in P. minus. RNA purity was quantified by A260:A280 ratios using ND-1000 Nanodrop spectrophotometer (Thermo Scientific) and the quality of RNA was visualized through agarose gel electrophoresis. RNA integrity numbers (RIN) of RNA samples for RNA-seq were determined by Agilent 2100 Bioanalyzer. RNA samples with RIN > 7 were sent to BGI (Shenzhen, China) for cDNA library preparation and RNA sequencing. Briefly, transcriptome sequencing of four cDNA paired-end libraries were constructed using TrueSeq RNA kit (Illumina) following the manufacturer's protocols, starting with 4 µg of total RNA.

RNA Sequencing and Assembly
Approximately 13.9 Gb of 192.17 million 90 bp paired-end reads were generated on Illumina HiSeq TM 2000 platform. The raw sequence data of each sample condition (Treated and Control) which associated with this study were deposited in GenBank under BioProject accession numbers PRJNA279327 and PRJNA208436 respectively. The raw reads were processed to trim off primer/adaptor sequences and filtered by quality (QV > 25) using Trimmomatic (Bolger et al., 2014). Subsequently, all the clean reads were de novo assembled using the wellestablished short-read assembling Trinity pipeline (http:// trinityrnaseq.github.io/) (Grabherr et al., 2011). Briefly, clean reads with a specified length of overlap were combined to produce longer contiguous sequences (contigs), and then these reads were mapped back onto the contigs with the paired-end information. The distance and relation among these contigs could be calculated based on paired-end reads, which allowed detection of contigs from the same transcript and also to calculate the distances among these contigs. Finally, the reconstructed contigs were further assembled and the constructed sequence that could not be extended on either end were considered as unique transcripts (UTs) .

Functional Annotation
The assembled UTs were annotated using BLASTX alignment against protein databases NCBI non-redundant (nr) (http:// www.ncbi.nlm.nih.gov) and SwissProt (http://www.expasy.ch/ sprot) (Magrane and Consortium, 2011), with an e-value cut-off < 10 −5 . For functional analysis of DEGs, the BLAST results based on hit score were imported into Blast2GO (Conesa et al., 2005) program to retrieve Gene Ontology (GO) annotation of DEGs for describing biological process, cellular component and molecular function categories. The distribution of DEGs functions was further classified with WEGO (Ye et al., 2006) software at the macro level. The DEGs sequences were also aligned to the Cluster of Orthologous Groups (COG) protein database (http:// www.ncbi.nlm.nih.gov/COG/) to predict and classify possible functions. Furthermore, in order to elucidate the biological pathways represented in P. minus DEGs, the sequences obtained from Blast2GO were searched against Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway maps database (Kanehisa and Goto, 2000).

Differential Gene Expression Analysis
RNA-seq by Expectation Maximization (RSEM) software package (Li and Dewey, 2011) was used to estimate transcript abundance. For each sample, the relative gene expression levels were normalized and expressed as FPKM (Fragments per Kilobase of Exon per Million reads Mapped) values. A Bioconductor package edgeR (Robinson et al., 2010) was used to analyze the differentially expressed genes by evaluating the dispersion of the entire dataset. UTs showing P ≤ 0.001, Benjamin-Hochberg False Discovery Rate (FDR) ≤ 0.05 and log2|fold change|> 2 were considered as significant differentially expressed genes (DEGs).

Functional Network Analysis and KEGG Pathway Enrichment Analysis
To assess major functional ontologies which are statistically overrepresented in the DEGs, global biological network analysis was performed using BiNGO (Biological Networks GO Ontology; v3.0.2), a Cytoscape plugin enrichment tool (Maere et al., 2005). GO-full terms within the BiNGO namespace with P ≤ 0.001 and Benjamini & Hochberg FDR correction were applied to depict the GO functional enrichment network and to show enriched functional categories in the up-and down-regulated DEGs.
To annotate putative pathways and biological functions involved in DEGs, the pathway enrichment analysis with hypergeometric test was carried out using "Annotate" and "Identify" programs in KEGG Orthology Based Annotation System (KOBAS 2.0) with Benjamini-Hochberg FDR correction (Xie et al., 2011). Up-and down-regulated DEGs were separately annotated against Arabidopsis database for enrichment analysis.

Identification of Transcription Factors
To identify the transcription factors (TFs) in P. minus, the assembled reads were searched against known TFs, as grouped in Plant Transcription Factor Database (PlnTFDB), using default parameters and cut-off E-value of 1e −5 . PlnTFDB is an integrative library of plant TFs which offers complete lists of TFs families that have fully sequenced genomes. Online protein sequences data for all species genes listed in the PlnTFDB database were downloaded from http://plntfdb.bio.uni-potsdam.de/v3.0/ website to align annotated UTs using local BLASTX. Moreover, a heatmap depicting the general trend of the differential expression of P. minus TF families in response to MeJA elicitation was created using MeV (http://mev.tm4.org/#/welcome) tool.

qRT-PCR Validation
The expression of 21 candidate genes was evaluated in relation to phenylpropanoid biosynthesis, JAs biosynthesis, JA signaling pathway and related TFs. Using PrimerBlast software, primers were designed to amplify short regions for each target and reference genes, ranging in size from 75 to 200 nt. Primers of 21 candidate UTs and reference gene are listed in Supplementary  Table 1. For each sample, 1 µg DNase-treated (DNA-free TM DNase; Ambion, Huntingdon, UK) RNA was reverse transcribed using iScript TM cDNA Synthesis Kit (Bio-Rad, CA) according to the manufacturer's instruction. qRT-PCR was carried out by using iTaq Universal SYBR R Green SuperMix kit (Bio-Rad). The amplification was executed with the following cycling program: 3 min at 95 • C, 40 cycles of 10 s at 95 • C, 30 s at 60 • C, and 30 s at 72 • C; and 0.06 s at 65 • C for reading. Real-time PCR was performed in the iQ TM 5 Real-Time PCR detection System (Bio-Rad, Hercules, CA). To find out the differences of relative fold for each sample, the CT value of each candidate genes was normalized to the CT values of reference gene CDPK (Calcium-Dependent Protein Kinase, Gene ID: comp58469_c0_seq3) and was calculated based upon the comparative CT (2 − Ct ) method (Schmittgen and Livak, 2008). The relative fold change for all candidate genes was set for both control and treated plants. Relative expression in RNA-seq was confirmed by three independent biological replicates and negative control of treated and mock-treated samples along with their three technical replicates.

RNA Sequencing, De novo Assembly and Annotation of Polygonum minus Leaf Transcriptome
We created a P. minus transcriptome dataset by sequencing RNA extracted from MeJA-treated and mock-treated leaf samples with two biological replicates each. The sequencing of four cDNA libraries on the Illumina HiSeq TM 2000 platform produced ∼13.9 G bases of total nucleotides (192,167,972 raw short reads; Supplementary Table 2). After trimming the adaptors and lowquality reads, we obtained 191,565,800 high quality paired-end reads from four samples.
To study MeJA-regulated gene expression, the trimmed reads were de novo assembled using Trinity pipeline, generating 182,111 UTs with an average length of 537 bp, an N50 of 1387 bp and a GC percentage of 43.64% (Supplementary Table 3).
For annotation of assembled UTs, BLASTX annotation algorithm was used with an E-value threshold of 10 −5 . The sequence homology search which is based on sequence similarities to the public protein databases, was carried out against the nr and Swiss-Prot protein databases. Of all the UTs, 36.43% showed significant homology with nr database; and the E-value distributions of the top hits in the Swiss-Prot database indicated that 29.45% of the mapped sequences showed significant similarity.

Differentially Expressed Genes in Response to MeJA
The abundance of UTs was estimated by mapping back the trimmed reads against the UTs with RSEM and their expression levels were represented as FPKM values. The two biological replicates of both control and MeJA-treated samples showed high correlation (R 2 > 0.79) in FPKM values (Supplementary Figure S1). Among 182,111 assembled UTs, a total of 2374 DEGs were found to be significantly differentially expressed at P ≤ 0.001, FDR ≤ 0.05, and log2|fold change|> 2 in response to MeJA treatment, of which 1419 were upregulated and 955 were down-regulated. This suggests that the MeJA-elicited response involves more gene activation than suppression.

Functional Annotation and Gene Ontology Classification of DEGs
For the prediction and classification of possible functions, all DEGs were aligned to the COG database (Tatusov et al., 2003). A total of 884 DEGs were classified into 23 COG functional categories, among which the category "Signal transduction mechanisms" (113, 12.78%) was predominant cluster that suggesting the differentially regulation of signal transduction pathways, particularly the transduction of plant hormones. In addition, a high percentage of genes were also assigned to "Amino acid transport and metabolism" (104, 11.76%) followed by "Carbohydrate transport and metabolism" (73, 8.26%) and "General function prediction only" (62, 7.01%). Only a few genes were assigned to the categories "Intracellular trafficking, secretion, and vesicular transport, " "Cytoskeleton, " "Cell motility, " and "Chromatin structure and dynamics." Analysis of COG classification indicated that the identified genes are involved in different biological processes. For example, 54 and 35 genes were respectively annotated to "Defense mechanisms" and "Secondary metabolites biosynthesis, transport and catabolism, " and are therefore associated with defense and biosynthesis of secondary metabolism pathways. Hence, the identification of genes in these pathways is required for the elucidation of the molecular mechanisms orchestrating stress response induction. The results of COG functional annotation of the P. minus DEGs are depicted in Figure 1. Further, DEGs with nr annotation was also annotated and classified under GO analysis to classify the functions of the predicted differential expressed genes in P. minus.
Using Blast2GO (Conesa et al., 2005), a total of 1178 DEGs were assigned with 205 GO terms. The majority of the top hits belonged to protein sequences of Beta vulgaris, followed by Vitis vinifera, Theobroma cacao, Nelumbo nucifera, and Jatropha curcas (Supplementary Figure S2). As shown in Figure 2, the assignments to biological processes (109, 53.17%) were mostly metabolic and cellular processes (710 and

Go-Based Functional Network Analysis of DEGs
The derived graphs from BINGO represent specific functions which were significantly enriched in three gene ontology categories. From the enrichment analysis, we identified biologically informative GO terms for comparison between the up-and down-regulated DEGs. The representative graph for the  up-regulated DEGs (Supplementary Figure S4A) indicated that most of the GO terms were significantly overrepresented in the defense and protective functions for plant response to chemical stimulus, stress and wounding. Secondary metabolic processes such as biosynthesis of phenylpropanoid, flavonoid, amino acid derivative, and aromatic compound were overrepresented. Concomitantly, under the molecular function category, DEGs related to "transferase, " "transporter, " "oxidoreductase, " "lyase, " and "transcription repressor" activities were significantly up-regulated 24 h after MeJA elicitation.
In contrast, majority GO terms of the down-regulated DEGs were significantly enriched in biological process related to photosynthesis, response to abiotic and biotic stimuli, respectively to light and bacterium, and generation of precursor metabolites and energy (Supplementary Figure S4B). In the molecular function category, GO terms related to "catalytic activity" of oxidoreductase and transferase, and "binding" of chlorophyll were represented in down-regulated DEGs.

KEGG Pathway-Based Gene Classification
According to sequence homology searches using BLASTX against the KEGG pathway database, 366 of 2374 DEGs were assigned to 3 main categories of 82 KEGG pathways representing compound biosynthesis, metabolism and degradation. In addition, 415 sequences of P. minus leaf were mapped to 193 enzymes (Supplementary Table 4). As shown in Table 1, the category with the greatest number of mapped entry was metabolism with the most in carbohydrate metabolism (89), followed by amino acid metabolism (49), energy metabolism (40), biosynthesis of other secondary metabolites (29), and lipid metabolism (27). Notably, the most represented pathways in secondary metabolism were those of phenylpropanoid and flavonoid biosynthesis (Supplementary Table 4). Only a few DEGs were mapped to translation (2) and signal transduction (1) categories.

KEGG Pathway Enrichment Analysis
KEGG pathway enrichment analysis was conducted using KOBAS to identify statistically significantly enriched pathways in response to MeJA elicitation. We identified 289 up-regulated and 180 down-regulated DEGs which were significantly enriched in eight and nine KEGG pathways respectively. The most significant with a majority number of up regulated DEGs involved in protective activities pathways such as biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, and phenylalanine metabolism ( Table 2). This is consistent with the finding from KEGG pathway mapping that greater number of up-regulated DEGs were found in phenylalanine, tyrosine and tryptophan biosynthesis (15), phenylalanine metabolism (13), phenylpropanoid biosynthesis (13), flavonoid biosynthesis (8), and α-linolenic acid metabolism (4) (Supplementary Table 5).
Conversely, photosynthesis, photosynthesis-antenna proteins as well as glyoxylate and dicarboxylate metabolism were amongst the significantly enriched pathways for down-regulated DEGs ( Table 2). In the KEGG pathway mapping, more down-regulated DEGs were also mapped to carbon fixation in photosynthetic organisms (13), fructose, and mannose metabolism (9), glycolysis/gluconeogenesis (9), pentose phosphate pathway (9), starch and sucrose metabolism (8), glyoxylate and dicarboxylate metabolism (7), and porphyrin and chlorophyll metabolism (4) (Supplementary Table 5). Such growth-defense tradeoffs presumably contribute to improvement in plant fitness through the diversion of saved energy from down-regulated photosynthesis to pertinent defense responses in changing environments (Rojas et al., 2014).

Identification of Differentially Expressed MeJA-Responsive Transcription Factors
A total of 2374 DEGs were searched against the Plant Transcription Factor Database (PlnTFDB) and 472 transcripts upon MeJA elicitation were significantly annotated as differentially expressed TFs. These TFs were grouped into 54 families including one orphan family containing 14 TFs which could not be grouped in any TF family. Most of the differentially expressed TFs are involved in plant stress responses, secondary metabolism and development processes, such as up-regulated expression of AP2/ERF, MYB, MADS, and NAC families, while bHLH and WRKY families were down-regulated (Figure 3). Members of these TF families have been studied and noted as candidate genes regulating the stress responses in different species (Van Verk et al., 2009;Yang et al., 2012a). The top 20 annotated differentially expressed TF families of P. minus leaf in response to MeJA treatment are shown in Supplementary Figure S5. In addition, our results indicated that among the up-regulated TF families (Supplementary Table 6), the MYB family contained the largest number of DEGs (19 members). MYB TFs represent one of the largest plant families that have important functions in various physiological processes under different conditions. These TFs are involved in primary and secondary metabolism, hormone signal transduction, defense response to environmental stresses and development process, etc. In recent years, several MYB genes were shown to be involved in the regulation of phenylpropanoid metabolism (Md-Mustafa et al., 2014;Liu et al., 2015). The Arabidopsis AtMB12 protein was identified as a specific regulator of phenylpropanoid biosynthesis of flavonol (Mehrtens et al., 2005). In tobacco MeJA signal transduction, NtMYBJS1 protein was shown to regulate some of early phenylpropanoid-related genes leading to the enhanced accumulation of phenylpropanoid conjugates under stress (Gális et al., 2006), and these results are consistent with our RNA-seq data. The most abundant down-regulated TF (Supplementary Table 6) belonged to FAR1 (FAR-RED IMPAIRED RESPONSE1) family. This family is regarded as plant-specific positive regulators of phytochrome A signal transduction pathway that able to absorb far-red (FR) light (Hudson et al., 1999). Interestingly, FR light related-responses genes were significantly down-regulated in MeJA-elicited P. minus leaf (Supplementary Figure S4).
Differentially expressed TFs which respond to MeJA treatment demonstrated that transcription regulation played a main role in MeJA-induced response network. Moreover, these results support previous studies which revealed the interaction of stress responses and developmental process (Chen et al., 2002;Cooper et al., 2003).

Validation of RNA-Seq Analysis by qRT-PCR
To evaluate our transcriptome result, we performed qRT-PCR to confirm the expression levels of 21 DEGs involved in the biosynthesis of phenylpropanoids and jasmonates, as well as JA signaling pathways and related TFs (Figure 4). qRT-PCR analysis revealed considerable activation for most of the candidate genes (except WRKY and NAC) in P. minus leaves upon elicitation which is consistent and correlates well (R 2 > 0.82) with the RNA-seq data (Figure 5).

DISCUSSION
Since plants produce various volatile secondary metabolites to cope with changing environments against oxidative and other stresses, plant chemical elicitors such as MeJA was applied to induce a defense response. Understanding the molecular mechanism underlying MeJA that triggers plant defense response is key to enhancing secondary metabolite production. P. minus remains understudied despite being one of the most commonly used aromatic and medicinal plants in Southeast Asia.
In a previous study, using a cDNA-AFLP approach, Ee et al. (2013) identified significant changes in gene expression of 52 transcript-derived fragments (TDFs) of P. minus leaf resulting from MeJA treatment, including those involved in stress response (6) and secondary metabolic processes (1). However, cDNA-AFLP study is hindered by the limited number of sequenced fragments. RNA-seq approach allows global and more accurate quantification of gene expression to obtain a better picture of gene regulation with more identified DEGs.
Leaf samples harvested 24 h after MeJA-elicitation were chosen for RNA-seq analysis based on results from metabolite profiling, which showed the greatest response (data not shown). We obtained 192 million clean reads which were assembled de novo into 182,111 UTs. The number of assembled UTs in P. minus is comparable to de novo transcriptome assembly of cotton, perennial ryegrass and Camellia sinensis using the same sequencing platform and assembler Farrell et al., 2014;Wu et al., 2014). The similarity searches of UTs against public protein databases found that majority (63.57%) UTs without functional annotation due to the lack of matching homologous sequences. These UTs with unknown function could be candidates of novel genes for future studies. FIGURE 3 | Hierarchical cluster analysis of differentially expressed putative TFs in P. minus leaf in response to MeJA. Several transcription factor families indicating differential expression are shown on the right side. The x-axis represents each replicate of control and MeJA treated conditions. The y-axis refers to TF genes expression levels. C, conrol; T, treated.
DEG analysis found 2374 UTs as statistically significant DEGs in MeJA-elicited P. minus leaf. DEGs were assigned to a broad range of COG classifications and GO categories suggesting that DEGs represent a variety of transcripts in P. minus leaf transcriptome. Among the 23 functional categories identified, we found the category "signal transduction mechanisms" to be the most represented in COG classification accounting for its need to induce the expression of specific subset of defense genes that result in the activation of the overall defense reaction (Pérez-Clemente et al., 2012). A large number of DEGs were annotated to the categories "Amino acid transport and metabolism, " "Carbohydrate transport and metabolism, " and "General function prediction only." Pathways involving in transport and metabolism of amino acid and carbohydrate are involved in plant growth and response to various stresses (Gupta and Kaur, 2005;Less and Galili, 2008). In addition amino acids are tightly linked to carbohydrate metabolism and serve as precursors for the synthesis of various classes of secondary metabolites involved in defense mechanism (Pratelli and Pilot, 2014).
GO-based functional network analysis of DEGs showed GO terms related to "response to stress" involving the biosynthesis of primary and secondary metabolites were enriched in the up-regulated DEGs, whereas terms related to "cell structure" and "growth" such as photosynthesis, response to radiation, light reaction and photosynthetic electron transport chain were enriched in the down-regulated DEGs (Supplementary Figure  S4). Light provides energy needed for photosynthesis and plays a role in photoregulation of plant growth and development (Diffey, 2013). These results support the idea that down-regulation of photosynthetic genes and their corresponding proteins as well as reduction in electron transport chain is a conserved feature of plant responses to many environmental stresses, however, the physiological significance of this phenomenon remains unclear and no experimental evidence is available to explain why it happens (Attaran et al., 2014;Rojas et al., 2014).
The KEGG database was searched to identify the biological pathways operating in MeJA-elicited P. minus leaves. Genes in the KEGG metabolism categories including starch and sucrose metabolism and pentose phosphate pathways with important functions in plant growth and development, as well as flavonoid and phenylpropanoid biosynthetic pathways important for stress responses were enriched. These secondary metabolites were previously reported to accumulate in high levels in sweet basil after MeJA treatment (Ku and Juvik, 2013;Misra et al., 2014).
The pathway enrichment analysis results of differentially MeJA-expressed genes can give helpful information to study specific functions, mechanisms and pathways that have important roles in MeJA response in P. minus leaves. Pathwaybased analysis highlights that genes within similar pathways commonly cooperate to operate certain biological functions (Wenping et al., 2011). This allows us to gain insight into the functions of these genes in the leaf transcriptome of P. minus. We found that genes related to the defense response were up-regulated, while the genes associated with growth and development processes were down-regulated upon MeJA elicitation. (Supplementary Figure S3). These collective findings from different functional annotation (KEGG, GO, and COG) of DEGs also support the view that the triggering defense responses by allocating resources from growth to defense lead to the negative effect on photosynthesis, and consistent with findings that plant stress hormone JAs plays a central role in controlling resource allocation between the competing processes of defense and growth (Attaran et al., 2014).
Most of the genes associated with lipid metabolism, such as alpha-linolenic acid metabolism, biosynthesis of unsaturated fatty acids, fatty acid biosynthesis/degradation, and linoleic acid metabolism displayed differential expression in response to MeJA. Four DEGs were found to be associated with the alpha-linolenic acid metabolism (Supplementary Table 4) that results in JA biosynthesis, including oxidase (EC: 1.3.3.6) and 13S-lipoxygenase (EC: 1.13.11.12) (Supplementary Figure S6). The unsaturated fatty acid alpha-linolenic acid is conserved in jasmonates production which is involved in transcriptional regulation of defense genes in response to environmental changes (Conconi et al., 1996). This result suggests the presence of a regulatory positive feedback mechanism for the JA biosynthesis (Sasaki et al., 2001).
Exogenous MeJA was reported by many studies to induce the expression of genes encoding JA biosynthetic enzymes and the JA signaling pathway in some plant species such as tomato, Arabidopsis and tobacco (Maucher et al., 2000;Chen et al., 2006;Kazan and Manners, 2008). In JA signaling pathway, the F-box protein CORONATINE INSENSITIVE 1 (COI1) bind to members of the JA ZIM domain (JAZ) repressor protein in the presence of JA-isoleucine (JA-Ile), and marks the complex for degradation by the 26S proteasome which releases the MYC2 TFs to activate transcription of JA-inducible genes (Wasternack and Hause, 2013).
Phytohormones such as JAs play key roles in plants responding to environmental cues to regulate plant growth and development under different stresses (Miao et al., 2015). Response of cellular functions to environmental cues requires differential gene expression, which is regulated by TFs. In plant cells, external signals modulate the levels of these regulatory elements which in turn regulate the transcription of defense genes in response to stress (Vom Endt et al., 2002), e.g., the NAC TFs have been shown to regulates JA-induced expression of defense genes in Arabidopsis thaliana (Bu et al., 2008). It has been suggested that stress-responsive TFs might be involved in stress tolerance (Saibo et al., 2009).
Since TFs play significant roles in stress signaling and plant defensive responses by translating signals from stresses into changes in gene expression (Lindemose et al., 2013), we executed comprehensive analysis of the differentially expressed TFs that regulate biosynthesis of MeJA-responsive genes. Based on the annotated UTs, 2.27% of the DEGs in the P. minus leaf encode TFs. As JA is essential in mediating plants defense responses and the induction of secondary metabolite synthesis (Shoji and Hashimoto, 2013), we hypothesized that those TFs involved in the regulation of defensive secondary metabolite biosynthesis can also induce JA genes. Interestingly, majority of P. minus TF families induced by MeJA elicitation mainly associated with either a specific hormones (e.g., jasmonic acid and salicylic acid) pathway such as AP2-EREBP family or a general stress response of MYB, NAC, MYC, bHLH, and WRKY families. These TF families were previously shown as jasmonate responsive transcriptional regulators of plant defensive gene biosynthesis (Chatel et al., 2003;Gális et al., 2006;Bu et al., 2008;Yu et al., 2012;Schluttenhofer et al., 2014), and suggesting their possible roles in the regulation of secondary metabolite biosynthesis in P. minus leaf. The largest class of TFs down regulated by MeJA elicitation was the FAR1 TF family. Photosynthesis and related metabolism are affected under environmental stresses (Saibo et al., 2009) and FAR1 family is known as a regulator of chlorophyll biosynthesis that serve as major pigments in photosynthesis and chloroplast division which is the exclusive location of photosynthesis (Tang et al., 2013;Wang and Wang, 2015). Thus, the result has been proposed that FAR1 family may have a role in down-regulating photosynthetic genes in P. minus leaf under MeJA treatment. Additionally, in our data, the largest class of TFs up-regulated by MeJA stress was the MYB TF family, suggesting that this TF family plays a major role in MeJA stress tolerance besides its role in the regulation of genes related to phenylpropanoid biosynthesis.
Among the plant secondary metabolites, phenylpropanoids are a group of aromatic amino acid phenylalanine derived physiologically active secondary metabolites such as, anthocyanins and isoflavonoids (Zhao et al., 2013). Phenylpropanoids contribute to the scents/aromas and play an important role in all aspects of plant responses against biotic and abiotic stimuli (Van Moerkercke et al., 2012). In our DEGs study, treatment by MeJA affected the pathways contributing to phenylpropanoid biosynthesis, phenylalanine metabolism, and also phenylalanine, tyrosine, and tryptophan biosynthesis.

CONCLUSION
In this study of MeJA-elicited P. minus leaf transcriptome, 182,111 UTs were assembled with 2374 DEGs identified. The is the first RNA-seq study on the response to MeJA of a non-model plant from Polygonaceae family, which comprises many useful medicinal plants. Functional annotation of DEGs based on COG, GO and KEGG showed the upregulation of UTs involve in plant defense and protective functions, as well as secondary metabolic processes. Conversely, UTs related to photosynthesis were downregulated. This is in accordance to the TF families found to be differentially expressed. MeJA elicitation is found to upregulate endogenous JA biosynthesis. The molecular mechanism on how MeJA triggers the upregulation of phenylpropanoid pathway is proposed and validated by qRT-PCR analysis. These results support the physiological trade-offs between stress response and growth despite that the underlying mechanism remains unclear. Furthermore, current P. minus transcriptome provides an extensive sequence resource for gene discovery and functional studies in related species.

AUTHOR CONTRIBUTIONS
NN and HG conceived the project. NN supervised the research work and contributed reagents/materials/analysis tools. HG designed the experiments, supervised the analysis, and critically revised the manuscript. KL conducted bioinformatics for raw RNA-seq data. RR collected plant materials, performed stress treatment, prepared samples for RNA-seq, conducted the gene functional analysis and TF analysis, performed qPCR based expression analysis and wrote the manuscript.