Arabidopsis Transcriptome Analysis Reveals Key Roles of Melatonin in Plant Defense Systems

Melatonin is a ubiquitous molecule and exists across kingdoms including plant species. Studies on melatonin in plants have mainly focused on its physiological influence on growth and development, and on its biosynthesis. Much less attention has been drawn to its affect on genome-wide gene expression. To comprehensively investigate the role(s) of melatonin at the genomics level, we utilized mRNA-seq technology to analyze Arabidopsis plants subjected to a 16-hour 100 pM (low) and 1 mM (high) melatonin treatment. The expression profiles were analyzed to identify differentially expressed genes. 100 pM melatonin treatment significantly affected the expression of only 81 genes with 51 down-regulated and 30 up-regulated. However, 1 mM melatonin significantly altered 1308 genes with 566 up-regulated and 742 down-regulated. Not all genes altered by low melatonin were affected by high melatonin, indicating different roles of melatonin in regulation of plant growth and development under low and high concentrations. Furthermore, a large number of genes altered by melatonin were involved in plant stress defense. Transcript levels for many stress receptors, kinases, and stress-associated calcium signals were up-regulated. The majority of transcription factors identified were also involved in plant stress defense. Additionally, most identified genes in ABA, ET, SA and JA pathways were up-regulated, while genes pertaining to auxin responses and signaling, peroxidases, and those associated with cell wall synthesis and modifications were mostly down-regulated. Our results indicate critical roles of melatonin in plant defense against various environmental stresses, and provide a framework for functional analysis of genes in melatonin-mediated signaling pathways.

Since the relatively recent discovery of melatonin in plants, investigations to elucidate its function in plants has been driven by what is known in animals. One such area of focus is the involvement of melatonin in modulating circadian rhythms and photoperiod-dependent processes. While melatonin levels do appear to be affected by light/dark cycles in some plants, the pattern varies among species, tissues, and organs [15][16][17][18][19][20][21]. The nocturnal peak of melatonin characteristic of humans and animals has been observed in Chenopodium rubrum [16], grape skins [17] and Ulva sp [21]. Conversely, studies conducted in water hyacinth demonstrated a peak in melatonin levels late in the day [19], indicating its biosynthesis in light. Furthermore, melatonin biosynthesis occurred under constant light in senescent rice leaves and was nearly undetectable under constant darkness [22]. Other reports found no significant correlation with melatonin levels and day/night cycles [18]. Interestingly, developing sweet cherries exhibited a dual peak of melatonin levels, one nocturnal and one in late day [20]. Contradicting reports of melatonin levels in ripening fruits add to the variation observed among plant species; melatonin levels decreased in ripening cherries [20], but increased in ripening tomatoes [18]. The possible role of melatonin in regulating flowering has also been investigated [23][24][25]; however an unequivocal role of melatonin in photoperiod-dependent processes in plants has not yet been established Melatonin has been studied extensively as an antioxidant in mammals. Many studies demonstrate the ability of melatonin to protect against many human diseases, including those linked to oxidative stress [26][27]. Melatonin was able to attenuate paraquat-induced lung and liver damage in rats [28][29] and Parkinson's disease in mice [30]. Furthermore, exogenously applied melatonin can enhance the production of antioxidative enzymes such as glutathione peroxidase and superoxide dismutase [31]. Melatonin may similarly play a protective role against oxidative stress in plants. Oxidative stress is capable of inducing elevated melatonin levels in various plant species [17,[32][33][34]. Indeed, the daytime peak of melatonin levels found in sweet cherry was associated with high temperature and light intensity, suggesting melatonin was synthesized in response to oxidative stress [17]. Transgenic rice seedlings with elevated levels of melatonin were more resistant to herbicide induced oxidative stress than their wild type counterparts [35]. Furthermore, oxidative stress induced the expression of genes involved in melatonin biosynthesis, leading to increased melatonin production in both wild type and transgenic rice [35]. Melatonin also appears to protect plants against UV and ozone damage [36][37][38][39][40], attenuate photo-oxidation of the photosynthetic system, and, at moderate levels, protect chlorophyll during senescence [39][40][41][42]. In addition, melatonin can promote low temperature and osmotic stress tolerance [43][44][45][46][47][48], alleviate copper damage [49][50], and improve salt tolerance [51] and fungal disease resistance [52] in a diversity of plant species.
The structure of melatonin is another feature that has driven investigations into its function in plants. Melatonin is structurally similar to the plant hormone indole-3-acetic acid (IAA) and has many features that make it a candidate for a functional auxin [53][54]. In addition, melatonin and auxin biosynthetic pathways share the same precursor, tryptophan [55]. Since auxins play critical roles as growth regulators during plant development such as shoot elongation, lateral root formation, and cell expansion, much work has focused on the effect of melatonin on these processes [42,48,[56][57][58][59][60][61][62][63]. Investigations have shown that melatonin and its precursor serotonin affect root growth in a dose-dependent manner similar to auxin [56][57][58][59][60]. At low melatonin levels, lateral root growth is stimulated, while at higher levels, adventitious root formation occurs and lateral root growth is inhibited in a mechanism seemingly independent of auxin [64]. Furthermore, melatonin has been demonstrated to stimulate expansion of etiolated lupin cotyledons [62] and promote hypocotyl growth [61,63] similar to IAA. It is still unknown whether the auxin-like affects are due to the action of melatonin itself or if melatonin is converted into IAA [64]. Moreover, while mammalian systems have well documented receptor-mediated gene expression, melatonin receptors have not been identified in plants and evidence points to a chemical response rather than a receptor-dependent response [64].
While much of the work conducted on melatonin in plants has focused on its physiological influence on growth and development, and on its biosynthesis, little work has focused on its affect on gene expression. Microarray analysis using endogenous melatonin-rich transgenic rice identified several hundred genes that are up-or down-regulated by elevated melatonin levels [65]. Previously, in an effort to understand mechanisms on how melatonin promotes lateral root formation in cucumber, we conducted mRNA-seq analysis using cucumber root tissues and identified potential clusters of genes that may control melatonin-mediated lateral root formation [66]. In this study, using Arabidopsis as a model, we utilized next generation RNA sequencing technology (RNA-seq) to obtain a comprehensive analysis of genome-wide changes in responses to external application of melatonin. RNA-seq can detect changes in gene expression with more precision than a standard microarray, allowing for the potential to identify novel genes. As a model species, Arabidopsis has many advantages for both basic and applied research, including easy transformation and ample resources of available T-DNA lines. Systemic analysis of the effect melatonin has on genome-wide gene expression in Arabidopsis will provide us basic information to genetically dissect melatonin-mediated signaling pathway(s) in regulating plant growth and development.

Plant materials and growth condition
Arabidopsis thaliana ecotype Columbia (Col-0) was used for RNAseq experiments and oxidative stress experiments. All plants were grown in a growth-chamber under 23uC and 14/10 hours light/ dark cycle. Three-week-old seedlings were removed from soil, rinsed, and submerged in100 pM, or 1 mM melatonin for 16 hours with gentle shaking. Mock solution was used as control. All experiments were duplicated for statistical analysis.

RNA extraction, library construction and sequencing
After melatonin treatment, total RNA was extracted using the QIAGEN RNeasy Mini Kit according to the manufacturer's instructions (QIAGEN). Purification of mRNA from total RNA was conducted using an Oligotex mRNA Mini Kit (QIAGEN). The mRNA was then used to construct cDNA libraries using the mRNA-Seq Sample Preparation Kit TM (Illumina) following standard protocols. Briefly, the mRNA was fragmented by exposure to divalent cations at 94uC and the fragmented mRNA was converted into double stranded cDNA. The cDNA ends were polished, the 39-hydroxls extended with A bases, and ligated to Illumina-specific adapter-primers. The adaptor ligated DNA was amplified by 15 cycles of PCR followed by purification using Qiagen TM PCR purification kit to obtain the final library for sequencing. The DNA yield and fragment insert size distribution of the library were determined on the Agilent Bioanalyzer. Library quantifications were performed by qPCR assays using the KAPA Library Quant Kit TM following the manufacturer's instructions.
All six constructed libraries were loaded on the Illumina flow cell at the appropriate concentration and bridge amplified to create millions of individual clonal clusters. The flow cells were sequenced on the HiSeq2000 sequencing instrument using 50 b single end protocols at the Center for the Study of Biological Complexity of the Virginia Commonwealth University.

Post-sequence analysis
After high throughput sequencing performance, short prematurely terminated sequences and those with low quality, as well as reads with any ambiguous bases were removed from raw sequence reads. Clean reads were aligned to the Arabidopsis genome assembly TAIR10 (as a reference) using TopHat (v1.3.1) BWA software [67]. Alignment results were stored in SAM/BAM file format [68] to permit downstream analysis. Cufflinks (v1.0.3) software [69] was then applied to identify transcripts and their abundance levels in RNA-seq data, as well as to perform differential gene expression analysis among multiple samples. Cufflinks calculates expression levels in fragments per kilobase of sequence per million fragments mapped (FPKM), which were used to perform differential gene expression analysis among treatments [69].
Gene Ontology classification of all genes exhibiting significant changes in transcript levels was conducted using AmiGO online tool [70]. Those with only CUFF numbers assigned during automatic analysis were manually checked against their positions using TAIR resources. Gene IDs and their putative functions were assigned accordingly.
Validation of mRNA-seq data using qRT-PCR A total of 60 genes were selected for verification of mRNA-seq data using qRT-PCR. Primer pairs for each selected gene were presented in Table S1. One microgram RNase free DNase treated total RNA from seedlings treated or untreated with melatonin was used to synthesize first strand cDNA using Superscript III (Invitrogen). qRT-PCR was performed using SSoAdvanced SYBR Green Supermix (BioRad) with the following parameters: 95uC for 3 minutes followed by 40 cycles of 95uC for 10 seconds and 60uC for 30 seconds. Gene expression was normalized via the Livak method using Arabidopsis Elongation Factor 1 (EF1; AT5G60390) as a reference gene [71].

Paraquat-induced oxidative stress
The ability of melatonin to attenuate paraquat-induced oxidative stress was examined using detached Arabidopsis leaves. Seedlings were grown under short day conditions (10/14 hour light/dark photoperiod) to promote vegetative growth for four weeks. Leaves were then detached and incubated in 0 mM, 10 mM or 50 mM paraquat in the presence or absence of 1 mM melatonin under 16/8 hour light/dark photoperiod. After 48 hours, leaves were analyzed for oxidative stress as visualized by photobleaching.

Results and Discussion
Transcriptome sequencing and identification of up-and down-regulated genes by melatonin To elucidate the roles of melatonin in regulating genome-wide gene expression, we conducted mRNA-seq analysis using Arabidopsis after 16 hours melatonin treatment. Because melatonin may function as a hormone at low concentration, and an antioxidant at high concentrations [35,37,56], we chose 100 pM (low) and 1 mM (high) melatonin to treat Arabidopsis to capture what may be physiological levels in planta and elevated levels. Mock solution was used as control. Biological duplicates for each treatment and control were included. After removing short prematurely terminated sequences, low quality sequences and any ambiguous bases, a total of more than 10 million clean reads with at least 106transcriptome coverage for each library (controls, 100 pM melatonin and 1 mM melatonin treatments) were generated through RNA-seq sequencing (Table 1). Of these, 87.1% and 88.3% for control, 88.0% and 88.6% for 100 pM treatment, and 87.0% and 87.9% for 1 mM treatment were mapped to the Arabidopsis genome assembly TARI10 using TopHat software [67]. Cufflinks was then used to assemble the aligned reads into transcripts and estimate the abundance of the transcripts to analyze differential expression among treatments [69]. In 100 pM melatonin-treated seedlings, only 81 gene transcript levels changed significantly (p,0.05). Of those genes, 30 were up-regulated, 51 were down-regulated. In addition, 4 genes were completely suppressed and 5 genes were uniquely  expressed compared to controls. Of the genes that exhibited a change in transcript levels at least 2-fold, 26 were up-regulated and 38 were down-regulated. However, when treated with 1 mM melatonin, more genes were identified with altered expression. A total of 1308 genes exhibited significant changes (p,0.05) in transcript levels in response to 1 mM melatonin with 566 genes up-regulated, and 742 genes down-regulated. Five genes were identified that were uniquely expressed in melatonin treated seedlings while 11 were completely suppressed by 1 mM melatonin. Nearly 900 genes exhibited changes in gene expression of at least 2-fold in response to 1 mM melatonin (387 upregulated, 510 down-regulated). The lists of these genes are presented in additional files as Table S2 for low melatonin affected genes and Table S3 for high melatonin affected genes. All these genes with significant changes in expression levels in response to low or high concentrations of melatonin were also identified using TAIR database and classified according to their functions using the Gene Ontology GO slimmer.

Reproducibility of RNA-seq profiles
To examine the variability among our RNA-seq experiments, all clean reads from melatonin treated and control sets were plotted for all possible pairs of independent experiments. Scatter plots of these data for control plants, 100 pM and 1 mM melatonin-treated plants are shown in Figure 1 and Figure S1. These scatter plots show that the most genes exhibit less variation in expression between the biological duplicates ( Figure 1) when compared to the scatter plots between treatments ( Figure S1). To further confirm its reproducibility, a phylogenetic analysis using sequences from all 6 libraries was conducted. Again, the biologically duplicated treatments are more closely related than those of the different treatments ( Figure S2). Taken together, these results indicate that our experiments are highly reproducible.

Validation of sequencing data by qRT-PCR
Although the results indicated that our experiments are highly reproducible, we further examined the reliability of the observed changes between treatments. qRT-PCR experiments were performed for a total of 60 genes exhibiting expression changes in response to melatonin treatments in the mRNA-seq analysis. The transcript levels were measured using the same RNAs used for the RNA-seq transcriptome analysis. The full list of 60 genes and their qRT-PCR validation is included in Table S4. Comparisons between mRNA-seq data and qRT-PCR for the top 15 up-and down-regulated genes from the selected group are shown in Figure 2. As shown, most of the selected genes indicated as differentially expressed by RNA-seq transcriptome profiling were confirmed by qRT-PCR. Only 7 of the selected genes were found to not significantly change in response to 1 mM melatonin, although the mRNA-seq data found significant changes in their transcript levels. However, relative fold changes observed in mRNA-seq data are sometimes not comparable to qRT-PCR analysis. This is mainly due to limitations and sensitivity of RNAseq analysis and global normalization methods. Nevertheless, the qRT-PCR analysis fit well with the mRNA-seq data, indicating the reliability of the RNA-seq data.  Table 2. Gene ontology classification into biological process of all genes significantly (p,0.05) differentially expressed in response to 1 mM and 100 pM melatonin. Genes were classified using AmiGO GO slimmer. The proportion of genes assigned to each category was calculated by dividing the number of genes assigned to a category by the total number of differentially expressed genes for each treatment. Genes may be classified into one or more GO Slim Term. *Indicates biological process is unknown. doi:10.1371/journal.pone.0093462.t002 Annotation and classification of differentially expressed genes into functional categories The Gene Ontology (GO) classification was conducted for all genes exhibiting a significant change in transcript levels using the AmiGO Slimmer tool and TAIR database. The Slimmer tool assigned general parent (GO slim) terms to the genes according to biological process, cellular component, and molecular function. The most striking difference between up and down regulated genes are those involved in response to stress and responses to endogenous and biotic stimuli. Of the up-regulated genes that were classified to a GO slim term, approximately 42% were involved in response to stress, 25% in response to endogenous stimulus, 24% in response to biotic stimulus, and 22% in response to signal transduction ( Table 2). In down-regulated genes, only 17% were involved in response to stress, 9% in response to endogenous stimulus, and 7% in response to both biotic stimulus and signal transduction. Interestingly, this trend was not observed in response to 100 pM melatonin. The proportion of genes involved in response to stress and biotic stimulus following 100 pM melatonin treatment was approximately 30% and 14%, respectively, for both up-and down-regulated genes. The trend was reversed in cell communication and signal transduction as a greater proportion of genes involved in those processes were down regulated in response to 100 pM melatonin. Genes involved in photosynthesis also trended towards down-regulation in response to 1 mM melatonin but trended towards up-regulation in response to 100 pM. Furthermore, genes associated with the cell wall accounted for 7% of down-regulated genes but only 3% of upregulated genes. 9% of genes down-regulated and 6% of the upregulated genes with changes of at least 2-fold were associated with redox pathways.
The biological processes that trended towards down-regulation in response to 1 mM melatonin included biosynthetic processes, metabolism of carbohydrates and nucleobase-containing compounds, development, cellular organization, morphogenesis, photosynthesis, and generation of precursor metabolites and energy.
The cellular components associated with genes down-regulated in response to 1 mM melatonin included cytoplasm, extracellular region, and, consistent with the down-regulation of photosynthesis associated genes, the thylakoid and plastids ( Table 3). The cellular components assigned to genes that trended towards up-regulation in response to 1 mM melatonin were the nucleus, plasma membrane, and the Golgi apparatus. The molecular functions associated with genes up-regulated in response to 1 mM melatonin included transferase activity, protein binding, and kinase activity, consistent with the general trend of signaling ( Table 4). The molecular functions that trended towards down-regulation in response to 1 mM melatonin were hydrolase activity, nucleic acid (both DNA and RNA) binding, lipid binding, and structural molecule activity. Interestingly, expression of chlorophyllase (CLH1), a light regulated enzyme involved in chlorophyll degradation [72], was significantly down-regulated in response to 1 mM melatonin (Gene ID: AT1G19670 in Table S3). This is consistent with a study conducted on senescing apple leaves where exogenous melatonin inhibited transcript levels of pheide a oxygnease (PAO), another key enzyme involved in chlorophyll degradation [41]. These findings can thus explain the means by which melatonin preserves chlorophyll content in leaves [39][40][41][42], delays senescence [39][40][41], and enhances photosynthetic rates [41]. Indeed, this may provide another mechanism by which melatonin preserves chlorophyll content during senescence in addition to attenuating ROS.
Due to the classic function of melatonin as a ''clock'' hormone in animals, one might infer that it would play a central role in timing of events in plants, such as flowering. However, since less than 1% of the genes affected by melatonin were involved in flowering, but nearly 40% were involved in stress responses and signaling, it is tempting to deduce that the central role of melatonin in plants is more likely as an antioxidant involved in stress response and as a photo-protectant. However, since the melatonin treatment was a one-time event and flowering is a complex programmed event that is promoted over several days, we cannot completely rule out that fluctuations in melatonin levels over time may signal day length and promote flowering.
Much fewer genes were affected by low level (100 pM) of melatonin treatment with only 81 genes significantly altered their expression. In view of the biological processes, most of the identified genes involved in cellular process, cell communication, response to stress, transporters and signal transductions are down- regulated. Cell death associated genes identified were also mostly down-regulated by 100 pM melatonin. Comparative analysis was conducted on genes affected by 100 pM melatonin and those by 1 mM melatonin. Results showed that not all genes down-regulated by low (100 pM) melatonin were down regulated by high (1 mM) melatonin. In fact, out of 51 of 100 pM down-regulated genes, only 18 genes were down regulated by 1 mM melatonin treatment, and 13 genes were upregulated by 1 mM melatonin treatment. Similarly, out of 30 genes up-regulated by 100 pM melatonin, only 5 were also upregulated by 1 mM melatonin and 8 were down regulated by 1 mM melatonin. Seventeen genes were uniquely up-regulated and 20 were uniquely down-regulated by low (100 pM) melatonin treatment (Figure 3). These results suggest that melatonin may play significantly different roles in control of plant growth and development under low and high concentrations. Overall, the results indicate that melatonin plays important physiological roles in many biological processes during plant growth and development such as stress defense, photosynthesis, cell wall modification and redox homeostasis.

Identification of novel transcriptionally active regions affected by melatonin
One of the advantages of mRNA-seq technique over microarray is that it can discover novel genes that were not previously annotated [73,74]. By mapping RNA-seq reads (contigs) against TAIR, we identified 13 contigs that were located in regions of the Arabidopsis genome where there were not predicted genes. Most of these contigs very strong FPKM signals, indicating they are

Melatonin altered many genes involved in plant defense
Given that a large group of genes altered by melatonin are involved in plant stress tolerance, we further categorized these genes according to their specific roles in plant defense. A schematic of the plants response to biotic and abiotic stresses to include key events and players involved in perception, signaling, and counterattack that were affected by melatonin is depicted in Figure 4. During biotic and abiotic stresses, stress perception involves Table 5. Receptors, kinases, and genes involved in calcium signaling whose expression levels were affected by 1 mM melatonin by at least 2 fold.  A closer examination of stress-responsive genes indicates that melatonin altered the expression of stress-response genes involved in all the steps along the way: from receptors through transcription factors. In addition, hormonal signaling was also examined and included in the schematic (Figure 4) due to the structural similarity between melatonin and IAA and the cross-talk among hormonal pathways and their involvement in signaling during stress events.

Melatonin-altered receptors, kinases, and calcium signaling genes pertaining to stress defense
When facing environmental stresses, plants first perceive stress signal and initiate a signal transduction cascade to defend against the pronounced stresses [75][76][77]. In order to understand the role of melatonin in plant stress defense, we identified all possible genes encoding receptors and kinases, and genes involved in calcium signaling with at least 2 fold changes in RNA-seq data. As shown in Table 5, all but 6 stress receptors identified exhibiting a change in transcript levels of at least 2 fold in response to melatonin were up-regulated. The plant immune response continued to trend towards up-regulation downstream of the receptors. The majority of kinases identified, including two MAPK and MAPKKK, were also up-regulated.
Calcium signals are also a core regulator for plant cellular responses to environmental stimuli [78][79][80]. Interestingly, all but two of the 14 genes identified involving calcium-dependent signaling in our RNA-seq data were up-regulated in response to melatonin. These genes encode calcium binding protein, calmodulin like protein, CBL-interacting protein kinases and calcium dependent phosphotriesterase. The involvement of calcium signaling in melatonin-mediated pathways is extensively documented in mammalian system [81][82][83]. The discovery of similar components of this pathway in plants indicates that melatonin may use similar signaling to implement its critical role(s) in both plants and mammals.

Transcription factors affected by melatonin
A total of 53 transcription factors were identified with at least 2 fold changes by melatonin treatment (Table 6). Interestingly, all but one of the 29 transcription factors identified as up-regulated genes are stress related transcription factors. These genes include 8 WRKY transcription factors, 5 NAC domain containing proteins and 5 zinc finger related transcription factors. The ERF transcription factor, CEJ1, is also up-regulated. CEJ1 mediates immune responses in Arabidopsis by repressing DREB [84]. In the melatonin down-regulated transcription factor group, 16 of the 24 genes are related to stress responses including 5 ERF (ethylene response factor) transcription factors. Two bZIP transcription factors are also partially suppressed by melatonin.

Melatonin altered gene expression along stress-induced hormone signaling pathways
Plant hormones and their cross-talk play important roles in both biotic and abiotic stress defense. Of these, abscisic acid (ABA) and ethylene (ET) are known for their abiotic stress defense functions [85][86][87], while salicylic acid (SA) and jasmonic acid (JA) are critical for plant biotic stress defense [88][89][90]. In addition, indole acetic acid (IAA) or auxin, also functions in plant stress defense systems [91][92]. In light of these stress-induced hormone signals, we identified 183 genes involved in hormone signaling exhibiting at least a 2-fold change in transcript levels in response to melatonin (Table S5). Of these, 52 genes pertaining to auxin responses and signaling were altered by melatonin with 29 down-regulated and 23 up-regulated. Most of the auxin responsive genes that were down-regulated in response to melatonin are involved in Auxin transport and homeostasis. Furthermore, one of the up-regulated genes encodes a GH3 protein, an IAA-amino synthase which conjugates amino acids to IAA, thus inactivating it [93]. These results may suggest that the Arabidopsis seedlings were responding to high levels of melatonin as they would respond to excess auxin.
Two ACC synthases were also up-regulated in response to melatonin, suggesting that melatonin may induce ethylene biosynthesis. One of the ACC synthases identified is ACS8, an auxin inducible ACC synthase. Indeed, this response is similar to what was induced by 2,4-D [94] and consistent with defense responses in plants. However, none of the known genes on auxin biosynthesis pathways [95] were significantly altered in their expression in response to melatonin (Table S3). This result is in agreement with findings that auxin inducible DR5:GUS reporter gene cannot be induced by melatonin and the function of melatonin is independent of auxin [60].
While the majority of auxin-responsive genes were downregulated in response to 1 mM melatonin, most genes on ABA, SA, JA and ET pathways were up-regulated (Table S5). On the ABA pathway, 36 out of 50 genes were up-regulated, 70 out of 92 were up-regulated on the SA pathway, 53 out of 67 were upregulated on the JA pathway, and 32 out of 42 were up-regulated on the ET pathway. As expected, some of these genes are induced by multi-hormone signals confirming important roles of crosstalk among hormones in plant defense systems. Many of the SA, JA, ABA and ET responsive genes induced by melatonin are also induced in response to biotic and abiotic stresses. Consistent with the changes in upstream gene expression levels, many downstream stress-associated genes were also affected, as shown in Table S6. These results further confirm the critical roles of melatonin in defense against both biotic and abiotic stresses in plants.  Table 7. Cell wall related genes with at least a 2 fold change in expression levels in response to 1 mM Melatonin.

Most cell wall associated genes were down-regulated by melatonin
Genes with functions involving the cell wall, including modification and growth, were largely affected by melatonin. Of the 60 genes associated with the cell wall, 45 were down-regulated and 14 were up-regulated at least 2-fold, and one gene was completely suppressed (Table 7). The down-regulated genes include 8 expansions and 4 pectin lyases or pectin methyltransferases. Two xyloglucan endotransglucosylases (XTH) were identified among up-regulated genes.
The results in the current study are significantly different from our previous study in cucumber roots showing that many of cell wall related genes, particularly pectinesterase genes, were upregulated by melatonin [66]. These discrepancies could be explained by different species and tissues being used or the dissimilar methods of melatonin treatment between the two studies. Therefore, it is necessary to further examine the role(s) of these cell wall related genes in melatonin-mediated signaling pathway. In general, cell wall modification and strengthening is the first line of defense plants have against pathogens [96][97]. Plant growth involving cell wall expansion would make cells more vulnerable to pathogen attack. Furthermore, XTHs are involved in strengthening the cell wall [98]. Thus, in the event of biotic stress, it is beneficial for plants to suppress cell wall expansion while strengthening the cell wall via XTHs. As a result, the cell wall will be more difficult for pathogens to breach. These discoveries are consistent with the studies reporting the role of melatonin in plant disease resistance [52].

Melatonin alters expression of redox associated genes
Genes and their protein products involved in various redox pathways are also used to protect cells from oxidative damage [99][100]. Many genes related to redox homeostasis were identified among the transcripts that were significantly differentially expressed when treated with melatonin (Table 8). In total, we identified 59 genes involved in redox homeostasis with significant changes in expression. Of these, 36 were down-regulated, and 23 were up-regulated by melatonin treatment. Most genes (53 of 59) identified were related to biotic and/or abiotic stress responses. Several members of the same families were differently affected by melatonin. For example, 7 of the 10 cytochromes, 6 of the 9 reductases, 4 of the 5 oxygenases, and 3 of the 5 oxidases were down-regulated by melatonin.
Peroxidases comprise another group of genes related to oxidative stress and are used as biochemical markers for plant resistance to bacterial and fungal diseases [101][102][103]. Thirteen genes encoding peroxidases were identified with significantly different expression levels following melatonin treatment ( Table 8). Eight of these genes were down-regulated by melatonin and five were up-regulated. All identified peroxidase genes were reported to be associated with plant stress defense. Of these, 8 (At1g05260, At1g14540, At3g21770, At3g49120, At4g08770, At4g08780, At4g12290, and At5g39580) were linked with plant biotic stress defense [104][105][106] and 5 (At1g05260, At1g14550, At1g49570, At4g30170, and At4g33420) were involved in plant abiotic stress defense [107][108]. The trend towards down-regulation of peroxidase genes in the current study is also contradictory to our previous study using cucumber roots [66] where most peroxidase genes were up-regulated by melatonin. Comparative analysis between these two studies indicates that most peroxidase genes identified in these studies belong to different members of the same gene family, with only two genes [At1G49570 (PER10), and At4G11290 (PER39)] mutual to both studies. Peroxidases belong to a large gene family with more than 70 members in plant species [109]. Their expression levels and patterns in different tissues vary among tissues [109] indicating they play important roles in a diversity of developmental processes [110][111].

Melatonin alleviates paraquat-induced photobleaching
Since the majority of genes with altered expression levels are involved in plant stress defense and melatonin is a potent antioxidant, we further examined the potential of melatonin to alleviate paraquat induced oxidative stress. Four-week old detached Arabidopsis leaves were incubated in 0 mM, 10 mM or 50 mM paraquat in the presence or absence of 1 mM melatonin. Leaves treated with paraquat in the absence of melatonin were completely photobleached after 48 hours under 16/8 hour light/dark photoperiod ( Figure 5). Leaves treated with 1 mM melatonin remained green, similar to leaves in the absence of paraquat. This result clearly demonstrated that melatonin can play critical roles in attenuating oxidative stress in plants. Since  melatonin is able to alleviate paraquat-induced disease in mammals [28][29][30], a similar role(s) of melatonin in defense against oxidative stresses may exist in both animals and plants. We can further utilize Arabidopsis as a model to elucidate the mode of action melatonin employs to protect against oxidative stressinduced diseases in humans and to further our understanding of its unique role(s) in plant species.

Conclusion
In conclusion, we report here the first comprehensive analysis of the effect of melatonin on genome-wide gene expression in Arabidopsis seedlings using RNA-seq technology. Given that Arabidopsis is an established plant model species and forward and reverse genetics methodologies are readily available, these datasets will provide fundamental information and serve as new tools to genetically dissect melatonin-mediated pathway(s) either common to both plants and animals, or unique to plants. Our transcriptome analysis reveals broader roles of melatonin in regulating plant growth and development. However, more importantly, melatonin may play critical role(s) in plant defense systems. Out of nearly 900 genes that were significantly up-or down-regulated by melatonin with at least 2 fold changes, almost 40% of the genes were related to plant stress defense, including many stress receptors, kinases,  and transcription factors, as well as downstream genes encoding end products that were directly used for stress defense. Furthermore, the expression of many genes involved in different hormone signaling pathways such as auxin, ABA, SA, ET and JA, and linked to plant stress defense, was also altered in response to melatonin treatment. Concurrently, expression of many cell wall associated genes, and genes involved in redox pathways, particularly peroxidases were significantly changed by melatonin treatments. Taken together, our results suggest that melatonin plays a critical role in plant defense against environmental stresses, including both biotic and abiotic stresses. Further dissection of the melatonin mediated pathway may lead to the development of novel strategies for crop improvement in the face of ubiquitous environmental stresses. Figure S1 Scatter plots between treatments.

Supporting Information
(TIF) Figure S2 Phylogenetic analysis of all clean RNA-seq data from six constructed cDNA libraries. (TIF) Table S1 List of genes and their primer pairs used for qRT-PCR validation. (DOCX)