## Histological and Quantitative Analysis of GFP

Histological analysis of GFP in rice tissues was detected and photographed under fluorescence microscope. Various tissues of T0 transgenic-positive transformants (root, leaf, sheath, panicle, stem, and mature seed) were sampled and observed under a fluorescence microscope (Leica MZ16F) using GFP filter sets and Leica Application Suite software.

The relative expression levels of GFP in rice tissues were detected by quantitative real-time PCR (qRT-PCR). Total RNAs of different rice tissues were extracted and reverse-transcribed as described previously (Wang et al., 2015a), and qRT-PCR was performed according to the same reference. The primers of GFP were GFP-F: 5′ -ATCCGCCACAACATCGAGGA-3′ and GFP-R: 5′ -TCGTCCATGCCGAGAGTGAT-3′ , and the primers

TABLE 2 | Polymerase chain reaction (PCR) primers used in this study.

*aThe underlined letters indicate the restriction enzyme sites.*

of GAPDH (the endogenous control) were GAPDH-F: 5′ - CTGCAACTCAGAAGACCGTTG-3′ and GAPDH-R: 5′ -CCT GTTGTCACCCTGGAAGTC-3′ . Relative expression levels were determined using 2−11C T method (Livak and Schmittgen, 2001).

### Melatonin Treatment

Melatonin (N-acetyl-5-methoxytryptamine), which is known as an indispensable hormone related to many physiological activities in animals, has also been identified as an important signaling molecule in response to many stresses in plants. In order to test the response of the four bidirectional promoters to melatonin, melatonin treatment was performed on the transgenic plants according to the procedure described by Shi and Chan (2014).

#### 5 ′ and 3′ Deletion Analysis of BIP1

Among the four bidirectional promoters above, BIP1 showed the highest expression efficiency in both 5′ and 3′ orientations. Therefore, it was selected for 5′ and 3′ deletion analysis in order to find the regulatory regions responsible for bidirectional expression activity (**Figure 1**). The specific primers used for PCR amplification to generate different 5′ and 3′ truncated fragments are shown in **Table 2**. Vector construction, callus culture and transformation, histochemical and fluorometric analysis of GUS activity, histological and quantitative analysis of GFP were performed as described above.

### Conservation Analysis of the Four Bidirectional Promoters and Bioinformatic Identification of Their Conserved Sequences

The conserved arrangements of the four gene pairs in six gramineous plants (O. sativa, Sorghum bicolor, Setaria italica, Brachypodium distachyon, Zea mays, and Triticum aestivum) were identified with the information from the Ensembl Plants database (http://plants.ensembl.org/index.html). The bidirectional genes whose homologous genes in other

species were still arranged in a bidirectional architecture were considered to be regulated by conserved bidirectional promoters (c-BIP); otherwise, they were considered to be regulated by non-conserved bidirectional promoters (n-BIP).

The conserved sequences in the four bidirectional promoters were discovered by MEME (http://meme-suite.org//tools/meme) and their frequencies in the intergenic regions between divergent gene pairs in rice genome were identified by FIMO (p < 1E-8) using the reference of random promoters (http://meme-suite.org//tools/fimo).

### RESULTS

### Selection of Four Novel Bidirectional Promoters in Rice Genome and Their Functional Characterization in Transgenic Plants

Based on RNA-seq and microarray data, we chose four divergent gene pairs (**Table 1**) according to the criteria in Section Methods and isolated their intergenic regions for functional identification. The four fragments were respectively cloned to a dual reporter vector pDX2181 and transformed into rice variety Zhonghua 11.

According to the results of GUS and GFP assays of the transgenic plants, all the intergenic regions showed bidirectional expression activity in various tissues. Histological GUS and GFP analysis of the transgenic plants showed that four novel bidirectional promoters (BIP1, BIP2, BIP3, and BIP4) were successfully identified in our work. Among them, BIP1, BIP2, and BIP3 showed bidirectional constitutive expression patterns and BIP4 showed bidirectional seed-specific expression pattern (**Figure 2**). Analysis of GUS fluorometric activities in various tissues of BIP1 transgenic plants (**Figure 3**) showed that the expression efficiency of BIP1 toward 3′ was the highest in the seed, which showed a GUS enzymatic activity of 17806 ± 2108 pmol 4-MU/min/mg protein, followed by in the root, which exhibited a GUS enzymatic activity of 14769 ± 1782 pmol 4- MU/min/mg protein, while the GUS enzymatic activities in the stem, sheath, panicle, and leaf were 9825 ± 1510, 8681 ± 834,

FIGURE 2 | Histological analysis of GFP and GUS expression in various tissues of the transgenic plants containing different GFP/bidirectional promoter/GUS fusions. (A–D), plants containing *GFP::BIP1*-*BIP4::GUS*. Localization of GFP is shown at the left area; Localization of GUS is shown at the right area.

7380 ± 895, and 6092 ± 875 pmol 4-MU/min/mg protein, respectively. Analysis of GFP expression in BIP1 transgenic plants (**Figure 3**) revealed that the expression efficiency of BIP1 toward 5′ was the highest in the panicle, which was 2.4-fold higher than that in the leaf. While the expression levels of GFP in the root, sheath, stem, and seed were 1.8-, 1.5-, 0.4-, and 0.1-fold higher than that in the leaf, respectively. GUS assays of BIP2 transgenic plants (**Figure 3**) indicated that the expression efficiency of BIP2 toward 3′ was the highest in the leaf, which showed a GUS enzymatic activity of 3713 ± 445 pmol 4-MU/min/mg protein, and the GUS activities in the sheath, panicle, stem, and seed were 2548 ± 339, 1755 ± 220, 1063 ± 145, and 872 ± 173 pmol 4-MU/min/mg protein, respectively, while GUS activity was hardly detected in the root. Analysis of GFP expression in BIP2 transgenic plants (**Figure 3**) showed that the expression efficiency of BIP2 toward 5′ was the highest in the stem, which was 8.6-fold higher than that in the panicle, while the expression levels of GFP in the root, sheath, leaf, and seed were 3.9-, 2.6-, 2.3-, and 1.6-fold higher than that in the panicle, respectively. BIP3 transgenic plants showed GUS enzymatic activities of 1126 ± 128, 841 ± 117, 806 ± 85, and 734 ± 71 pmol 4-MU/min/mg protein in the sheath, seed, stem, and leaf, respectively, while GUS activity was hardly detected in the panicle and root. The expression efficiency of BIP3 toward 5 ′ was the highest in the panicle, which was 5.7-fold higher than that in the leaf, while the expression levels of GFP in the seed, stem, root, and sheath were 3.6-, 2.9-, 1.3-, and 0.9-fold higher than that in the leaf, respectively. BIP4 showed a bidirectional seed-specific expression pattern, as a GUS enzymatic activity of 1659 ± 211 pmol 4-MU/min/mg protein was detected in the seed while almost no GUS activity was detected in other tissues, and a 23-fold higher expression level of GFP was observed in the seed compared with in the leaf. In addition, these results indicate that all the bidirectional promoters identified here direct gene expression in an orientation-independent manner; namely, the expression patterns in opposite directions of these bidirectional promoters are similar, which is consistent with the co-expression characteristics of the adjacent genes.

Melatonin is one of the most important hormones in plant and animal. In order to test the response of the four bidirectional promoters to melatonin, melatonin treatment was performed on the BIPs transgenic plants. The results of GUS and GFP assays indicated that the four bidirectional promoters were not induced by melatonin (**Figure 4**). Hence, it can be inferred that no melatonin-responsive cis-element was harbored in these promoters.

### Identification of the Expression Regulatory Regions in BIP1

In order to identify the regulatory regions in BIP1, 5′ and 3 ′ deletion analysis of this promoter was performed. A series of truncated BIP1 were respectively cloned to pDX2181 and transformed into Zhonghua 11. Transgenic plants carrying BIP1- 1F2R showed much lower GUS activity than BIP1 transgenic plants in various tissues, especially in the root, stem, seed and panicle, whose GUS activities were lower than 10% of that in the corresponding tissues of BIP1 transgenic plants (**Figures 5**, **6**). It thus could be inferred that region 1 could greatly increase the transcriptional activity of BIP1 toward 3′ . Meanwhile, BIP1- 1F2R transgenic plants also showed an obvious decrease of GFP expression in various tissues compared with BIP1 transgenic plants (**Figure 6**). These results could be integrated to reveal that region 1 is a bidirectional transcription-enhancing region of BIP1. Further truncating in 3′ of BIP1 led to complete abolishment of GUS activity in BIP1-1F3R transgenic plants, while the GFP expression level was not obviously reduced in BIP1-1F3R transgenic plants compared with in BIP1-1F2R transgenic plants. These results suggest that region 2 is the essential region responsible for the basic expression activity of 3′ but not for that of 5′ . Transgenic plants carrying BIP1- 2F1R or BIP1-2F2R showed no expression of GFP, indicating that truncating 5′ of BIP1 will completely abolish 5′ expression activity of the promoter. GUS assays in BIP1-2F1R transgenic plants revealed that truncating 5′ of BIP1 caused slight changes of 3′ expression activity in most tissues except for the root, which showed obviously decreased GUS activity compared with that of BIP1 transgenic plants (**Figure 6**). These results suggest that region 3 is the essential region responsible for the basic expression activity of 5′ but not for that of 3′ ; however, it can positively regulate the expression activity of 3′ in the root.

### Conserved Arrangement of the Four Gene Pairs Regulated by Bidirectional Promoters

The sequences of promoters are known to be variable (Müller et al., 2007). Therefore, in order to analyze the conservation of the four bidirectional promoters in different species, we investigated the conservation of the four gene pairs regulated by these promoters.

The conserved arrangement of the four gene pairs in six gramineous plants was identified with information from the Ensembl Plants database. The bidirectional genes whose homologous genes in other species were still arranged in a bidirectional architecture were considered to be regulated by c-BIP; otherwise, they were considered to be regulated by n-BIP (**Table 3**). It was found that BIP1 and BIP3 were the most c-BIP in the six gramineous plants. BIP1 was conserved in O. sativa, S. bicolor, B. distachyon, and Z. mays; BIP3 was conserved in O. sativa, S. bicolor, S. italic, and B. distachyon; BIP4 was conserved in O. sativa, S. bicolor, and S. italica; while BIP2 was only conserved in O. sativa and T. aestivum.

### Potential cis-Sequences Involved in Bidirectional Expression

By MEME, two conserved cis-sequences in the four bidirectional promoters were identified (**Figure 7**). Cis-sequence 1 was a G/C-rich sequence, while cis-sequence 2 was an A/T-rich sequence. Subsequently, the frequencies of these cis-sequences in the intergenic regions between divergent gene pairs in rice genome were analyzed by FIMO using the reference of random promoters. Consistent with the expectation, the two cissequences conserved in the four bidirectional promoters both showed overrepresentation in potential bidirectional promoters in rice genome compared with random promoters. This result further reveals that the two novel cis-sequences are probably involved in bidirectional expression.

### DISCUSSION

In this study, we initiatively combined RNA-seq data and cDNA microarray data to discover potential bidirectional promoters in rice. Four adjacent and oppositely transcribed gene pairs were selected based on their expression levels and correlations.

*GFP::BIP1*-*BIP4::GUS*. a: no significant difference. Error bars indicate SE based on five independent biological replicates.

The intergenic regions between the four gene pairs were successfully cloned, and all of them were identified to be bidirectional promoters, confirming the feasibility of our method for discovering bidirectional promoters. This is the first study to clone and identify bidirectional promoters using two reporter genes simultaneously with stable transformation in rice. Among the four identified bidirectional promoters, BIP1 shows high expression efficiency in various tissues, and thus has a high application potential in genetic engineering, such as driving two resistant genes simultaneously in transgenic breeding against pest/disease stress, which can confer more strong, broad, and durable resistance in rice (Du et al., 2009; Shah et al., 2009; Yang et al., 2011; Wang et al., 2015a). Rice is one of the most important food crops in the world, and its seed is the edible part

version/GUS fusions. 1F2R, plants containing *GFP::BIP1*-*1F2R::GUS*; 1F3R, 2F1R and 2F2R follow the same pattern. a, b, c, d: significant difference (*P* < 0.05). Error bars indicate SE based on five independent biological replicates.

*If 5*′ *gene/3*′ *gene of the bidirectional promoters had homologous gene in another species (C/C) and the homologous genes were still arranged in bidirectional architecture, they were considered to be regulated by conserved bidirectional promoters (c-BIP). Otherwise, they were considered to be regulated by non-conserved bidirectional promoters (n-BIP).*

consumed by human. Therefore, it is highly necessary to improve the nutrient quality of the seed (Ha et al., 2010; Li Y. et al., 2011, 2014; Ogo et al., 2013). Efficient and specific expression of multiple target genes for seed improvement in rice could hardly be realized without seed-specific promoter. In this work, BIP4 shows a bidirectional seed-specific expression pattern, indicating its high application potential in the improvement of seed quality by specifically driving multiple genes. The results of 5′ and 3 ′ deletion analysis reveal that region 1 is the bidirectional transcription-enhancing region of BIP1; region 2 is the essential region specifically responsible for the basic expression activity of 3′ ; region 3 is the essential region responsible for the basic expression activity of 5′ but not for that of 3′ , while it can positively regulate the 3′ expression activity in the root.

Conservation analysis of the four bidirectional promoters in gramineous plants reveals the possible co-evolution of adjacent genes regulated by these promoters. The bidirectional arrangement of LOC\_Os02g42314 and LOC\_Os02g42320,

which are regulated by BIP1, is conserved in four species of six gramineous plants, suggesting that they are relatively conserved during co-evolution. The functional annotations of these two genes, which are "ubiquitin-conjugating enzyme" and "peptidase," show that both of them are structural genes conserved during evolution. Moreover, the functional relationship between the two genes further supports their co-evolutionary conservation. LOC\_Os02g47000 and LOC\_Os02g47010, which are regulated by BIP3, also show conserved arrangement in four species of six gramineous plants. Although the functional annotation of LOC\_Os02g47000 is unclear, considering that LOC\_Os02g47010 is annotated to encode "secretory carrier-associated membrane protein," we speculate that LOC\_Os02g47000 may encode a structural protein related to secretory pathway.

 A previous report also has suggested that several

### REFERENCES

known cis-sequences might be related to bidirectional expression (Dhadi et al., 2009). Here, we used the experimentally verified bidirectional promoters to predict two cis-sequences related to bidirectional expression which had not been identified in rice genome. Subsequently, overrepresentation of the two novel cis-sequences in the intergenic regions between divergent gene pairs further reveals their involvement in bidirectional expression. Interestingly, cis-sequence 1 is a G/C-rich sequence, which is consistent with the characteristics of higher GC content in bidirectional promoters; however, cis-sequence 2 is an A/T-rich sequence, which might be a new finding in the sequence characteristics of bidirectional promoters. Overall, the novel bidirectional promoters identified using two reporter genes simultaneously with stable transformation in rice are expected to have high applicability in genetic engineering. Our study proposes a feasible method for selecting, cloning, and functionally identifying bidirectional promoters as well as for discovering their bidirectional regulatory regions and conserved sequences in rice.

### AUTHOR CONTRIBUTIONS

YL and RW conceived and designed the experiments. RW, YY, MZ, and MY performed the experiments. RW and YY performed the data analysis. RW and YL wrote the paper. FZ and HC revised the paper. YL secured the funds to support this research.

### ACKNOWLEDGMENTS

This research was supported by the National Program of Transgenic Variety Development of China (2014ZX08001-001), the National High Technology Research and Development Program of China (863 Program) and the National Natural Science Foundation of China. And we thank Prof. Zuoxiong Liu (Huazhong Agricultural University) for language improvement.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Wang, Yan, Zhu, Yang, Zhou, Chen and Lin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Comparative Transcriptional Profiling of Melatonin Synthesis and Catabolic Genes Indicates the Possible Role of Melatonin in Developmental and Stress Responses in Rice

Yunxie Wei1† , Hongqiu Zeng1† , Wei Hu2† , Lanzhen Chen3 , Chaozu He1 \* and Haitao Shi1 \*

1 Hainan Key Laboratory for Sustainable Utilization of Tropical Bioresources, College of Agriculture, Hainan University, Haikou, China, 2 Key Laboratory of Biology and Genetic Resources of Tropical Crops, Institute of Tropical Bioscience and Biotechnology, Chinese Academy of Tropical Agricultural Sciences, Haikou, China, 3 Institute of Apicultural Research, Chinese Academy of Agricultural Sciences, Beijing, China

### Edited by:

Jin-Gui Chen, Oak Ridge National Laboratory, USA

#### Reviewed by:

Ruediger Hardeland, University of Goettingen, Germany Kyoungwhan Back, Chonnam National University, South Korea

### \*Correspondence:

Haitao Shi<EMAIL_ADDRESS>Chaozu He<EMAIL_ADDRESS>
†These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Plant Cell Biology, a section of the journal Frontiers in Plant Science

Received: 23 February 2016 Accepted: 02 May 2016 Published: 18 May 2016

### Citation:

Wei Y, Zeng H, Hu W, Chen L, He C and Shi H (2016) Comparative Transcriptional Profiling of Melatonin Synthesis and Catabolic Genes Indicates the Possible Role of Melatonin in Developmental and Stress Responses in Rice. Front. Plant Sci. 7:676. doi: 10.3389/fpls.2016.00676 As a well-known animal hormone, melatonin (N-acetyl-5-methoxytryptamine) is also involved in multiple plant biological processes, especially in various stress responses. Rice is one of the most important crops, and melatonin is taken in by many people everyday from rice. However, the transcriptional profiling of melatonin-related genes in rice is largely unknown. In this study, the expression patterns of 11 melatonin related genes in rice in different periods, tissues, in response to different treatments were synthetically analyzed using published microarray data. These results suggest that the melatonin-related genes may play important and dual roles in rice developmental stages. We highlight the commonly regulation of rice melatonin-related genes by abscisic acid (ABA), jasmonic acid (JA), various abiotic stresses and pathogen infection, indicating the possible role of these genes in multiple stress responses and underlying crosstalks of plant hormones, especially ABA and JA. Taken together, this study may provide insight into the association among melatonin biosynthesis and catabolic pathway, plant development and stress responses in rice. The profile analysis identified candidate genes for further functional characterization in circadian rhythm and specific stress responses.

Keywords: melatonin, rice, gene expression, circadian rhythm, development, immunity, stress response

### INTRODUCTION

Melatonin (N-acetyl-5-methoxytryptamine) was first discovered in the cow's pineal gland (Lerner et al., 1958). Dubbels et al. (1995) and Hattori et al. (1995), melatonin was identified in plants by two research groups. Until now, melatonin has been found in multiple plant species, including alfalfa, almond, anise, apples, Arabidopsis, banana, beetroot, bermudagrass, black mustard, cabbage, celery,

**Abbreviations:** 2-ODD, 2-oxoglutarate-dependent dioxygenase; AANAT, arylalkylamine N-acetyltransferase; ABA, abscisic acid; ASMT, N-aceylserotonin methyltransferase; AXR3, Auxin Resistant 3; BL, indole-3-acetic acid; DAT, days after transplanting; GA, gibberellic acid; GEO, Gene Expression Omnibus; hpi, hour post inoculation; IAA, indole-3-acetic acid; JA, jasmonic acid; M2H, melatonin 2-hydroxylase; SNAT, serotonin N-acetyltransferase; T3S, type III secretion system; T5H, tryptamine 5-hydroxylase; TDC, tryptophan decarboxylase; tZ, trans-zeatin; Xoo, Xanthomonas oryzae pv. oryzae.

cherry, coriander, cucumber, fennel, fenugreek, flax, green cardamom, milk thistle, oranges, poppy, potato, rice, sunflower, tobacco, tomato, white mustard, wolf berry, etc. (Manchester et al., 2000; Zhao et al., 2013; Shi and Chan, 2014). In addition, the endogenous melatonin concentration can also be modulated through genetic transformation in tomato and rice (Okazaki and Ezura, 2009; Okazaki et al., 2009, 2010; Byeon et al., 2012, 2013, 2014; Byeon and Back, 2014a,b).

To date, the biosynthesis and metabolic pathways of melatonin in plants have been established (**Figure 1**). Melatonin in plants can be synthesized by four sequential enzymes from tryptophan (Kang et al., 2011), including TDC, T5H, SNAT, and N-aceylserotonin O-methyltransferase (ASMT) (Arnao and Hernández-Ruiz, 2014, 2015; Zuo et al., 2014). Thereafter, melatonin is catabolized by M2H into 2-hydroxymelatonin (Byeon and Back, 2015). In rice, gene families of TDC, T5H, SNAT, and ASMT contain 3, 1, 1, and 3 known members, respectively (Kang et al., 2007; Fujiwara et al., 2010; Kang et al., 2013; Park et al., 2013a). However, OsASMT3 is barely detectable in any of the plant organs (Park et al., 2013b). OsM2H genes belong to 2-ODD family and at least 4 of 2-ODD genes show M2H activities in rice (Byeon and Back, 2015).

Melatonin plays protective roles in the regulation of plant tolerance to abiotic stress and biotic stress (Yin et al., 2013; Lee et al., 2014, 2015; Zhan et al., 2015). Overexpression of OsTDC increases endogenous melatonin level and delays leaf senescence in rice (Kang et al., 2007, 2009; Byeon et al., 2014). The transcript of OsT5H can be induced by Magnaporthe grisea infection (Fujiwara et al., 2010). Exogenous application of serotonin, the penultimate substrate for melatonin biosynthesis, induces defense gene expression and increases resistance to rice blast infection (Fujiwara et al., 2010). Transgenic rice plants ectopically expressing the AANAT regulates cold stress resistance (Kang et al., 2010), seminal root elongation (Park and Back, 2012), oxidative stress resistance (Park et al., 2013b), and seedling growth (Byeon and Back, 2014a). The transcript of OsASMT can be induced by ABA and methyl JA treatments, and OsASMT overexpressing plants result in higher level of melatonin (Park et al., 2013b). Exogenous application of melatonin improved apple resistance to Marssonina apple blotch (Diplocarpon mali) (Yin et al., 2013), enhanced disease defense against Pseudomonas syringae DC3000 in Arabidopsis and tobacco (Lee et al., 2014; Shi et al., 2015a).

Rice is one of the most important crops around the world, and melatonin is also taken in from rice by many people everyday. Thus, it is very useful and important to dissect the distribution and regulation of endogenous melatonin in rice. Melatonin is widely involved in plant development, multiple

abiotic and biotic stress responses in Arabidopsis (Shi and Chan, 2014), and Bermudagrass (Shi et al., 2015b). However, transcriptional profiling of rice melatonin synthesis and catabolic genes has not been systematically carried out. In this study, we analyzed the expression profiling of 11 rice melatonin synthesis and catabolic genes in development, various tissues, and in response to hormone, pathogen infection, drought, salt, and cold stresses. These results may provide insight into the link among melatonin biosynthesis and catabolic pathway, plant development and stress responses in rice. Further functional characterization of identified candidate genes with potential involvement in circadian rhythm and stress responses through overexpressing, knocking down or knocking out will give more clues to melatonin-mediated signaling as well as underlying molecular mechanism.

### MATERIALS AND METHODS

### Plant materials and Growth Conditions

Rice (Oryza sativa L. ssp. japonica cv. Nipponbare) seeds were sown in germinating boxes. At 30 days after germination, the seedlings were transplanted in a paddy field under normal conditions of the cultivation season. Thereafter, 56 DAT, 58 DAT, and 90 DAT were considered as the stage of panicle initiation, the early stage of panicle development indicating a complete reproductive transition, the stages of flowering and early stages of seed development corresponding to the ripening-stage transition, respectively.

For hormone treatments, rice seeds were germinated, and grown hydroponically in a growth chamber at 28◦C under continuous light. Seven-day old seedlings were transferred in culture solution containing 50 µM ABA, or 10 µM GA, or 10 µM IAA, or 1 µM brassinolide (BL), or 1 µM tZ, or 100 µM JA, or in culture solution without hormone to serve as control (mock treatment). Samples were collected after 0, 0.25, 0.5, 1, 3, and 6 h incubation for root, and after 0, 1, 3, 6, and 12 h incubation for shoot.

### Pathogen Infection

Rice (O. sativa cv Nipponbare) plants grown in the greenhouse for 42 days were inoculated with Xoo T7174R, a wild-type strain, and 74HrcV::Km, a T3S-defective mutant by the leafclipping method. Plants treated with water were used as control. Leaf sections (3–5 mm) that included the inoculation site were collected at 3, 6, and 12 hpi and 1, 2, 4, 6 dpi.

For the blast fungus infection, three lines of rice cultivar Nipponbare carrying the blast resistance genes (Pia, Pish) were inoculated with two strains of Magnaporthe oryzae harboring AVR-Pia and AVR-Pish. Rice seedlings at the 4-leaf stage were placed in moist chamber and sprayed with a conidial suspension of M. oryzae (1 × 106 conidia/ml). The seedlings were incubated in a moist chamber at 25◦C for 24 h under dark condition, then grown in hydroponic culture under 14 h light (28◦C) and 10 h dark (24◦C). Leaf samples (4th leaf) from 3 individual experiments were harvested at 1, 2, 3, and 5 days post inoculation (dpi). Rice seedlings sprayed with water were used as control.

### Development- and Pathogen Infection-Related Data Analysis

The data of spatio-temporal transcript levels in various tissues or organs (RXP\_0001), leaf and root transcriptional profile in light (RXP\_003 and RXP\_007) and dark (RXP\_004 and RXP\_008) throughout entire growth in the field, diurnal, and circadian leaf (RXP\_002) and root (RXP\_009) transcriptional profile throughout entire growth, plant hormone profile (RXP\_001 to RXP\_012), Xoo-treated profile (RXP\_3002), and M. oryzae-treated profile (RXP\_3001) were downloaded from RiceXPro1 (Sato et al., 2011a,b, 2013). All samples were used for hybridization using the Agilent one-color (Cy3) microarraybased gene analysis system. As detailed described in Sato et al. (2013), all the above data were deposited in GEO through the following accession numbers: GSE21396, GSE21397, GSE36040, GSE36042, GSE36043, GSE36044, GSE39423, GSE39424, GSE39425, GSE39426, GSE39427, GSE39429, and GSE39432. All the raw data were downloaded and re-analyzed for cluster analysis of expression profile that shown as normalized data (log2).

### Abiotic Stress-Related Data Analysis

As described in Jain et al. (2007), 7-day-old light-grown rice seedlings were transferred to control condition and 200 mM NaCl solution as salt stress for 3 h, were dried on tissue paper as dehydration stress for 3 h, and were kept at 4◦C as cold stress for 3 h. Then the seedlings were sampled in triplicate. GEO series accession no. GPL2025 were used for microarray analysis as Jain et al. (2007) described. All the normalized data were obtained from Rice eFP Browser2 (Jain et al., 2007).

### Cluster Analysis

The original data from RiceXPro and Rice eFP Browser were listed in **Supplementary Table S1**. Hierarchical cluster analysis of transcriptional profile was performed using CLUSTER program3 (Larkin et al., 2007), and the heatmap was obtained using Java Treeview4 (Saldanha, 2004) according to the instructions.

### RESULTS

### The Spatio-Temporal Transcript Levels of Rice Melatonin Synthesis and Catabolic Genes in Various Tissues or Organs

To investigate the expression profiles of rice melatonin synthesis and catabolic genes in various tissues or organs, we analyzed the expression of these genes using published microarray data (Sato et al., 2013). Eleven of rice melatonin-related genes have the

1http://ricexpro.dna.affrc.go.jp/

2http://bar.utoronto.ca/efprice/cgi-bin/efpWeb.cgi

3http://bonsai.hgc.jp/∼mdehoon/software/cluster/software.htm 4http://jtreeview.sourceforge.net/

corresponding probe sets in the dataset (**Figure 1**). As shown in **Figure 2**, all genes showed different expression pattern in various tissues, indicating that these genes may play different roles in plant growth and development. Interestingly, the expression patterns of melatonin-related genes could be divided into two groups (**Figure 2**). One group contained six genes (OsT5H, Os2-ODD11, −19, −21,−33, and OsTDC1), and most of them showed high expression levels in endosperm tissues. The other group contained five genes (OsASMT1, −2, OsTDC2, −3, and OsSNAT), and most of them showed high expression levels in leaf blade and leaf sheath tissues. Moreover, Os2-ODD19 showed high expression level in four time points of anther, while other four genes (OsTDC1, −3, Os2-ODD11, and OsASMT1) exhibited lower expression level. Similarly, six genes (Os2-ODD11, −19, −21,−33, OsTDC1, and OsASMT1) showed high expression level in five time points of endosperm, while two genes (OsASMT2 and OsSNAT) with a relative low level of expression. These melatoninrelated genes showed high expression levels in a special tissue indicated their possible roles of melatonin in special tissue.

Moreover, three genes (Os2-ODD19, −21, and OsTDC2) showed different expression level in day and night at 9 tissues (leaf blade-vegetative, leaf blade-reproductive, leaf blade-ripening, leaf sheath-vegetative, leaf sheath-reproductive, root-vegetative, rootreproductive, stem-reproductive, and stem-ripening). The results indicate that Os2-ODD19, −21, and OsTDC2 may play some roles in circadian rhythm and may be used in further functional characterization.

### Transcriptional Profile of Rice Melatonin Synthesis and Catabolic Genes throughout Entire Growth in the Field

In rice leaves at day, the transcript levels of OsTDC3 and OsASMT2 were increased at vegetative stages, while that of OsT5H was decreased (**Figure 3A**). At reproductive stages, the transcripts of seven genes (OsTDC1, −3, OsASMT1, −2, Os2-ODD11, −33, and OsSNAT) and two genes (OsTDC2 and Os2-ODD21) showed up-regulation and down-regulation,

fpls-07-00676 May 14, 2016 Time: 12:18 # 4

respectively (**Figure 3A**). At ripening stages, the transcripts of three genes (Os2-ODD11, −33, and OsASMT1) and two genes (OsTDC2 and Os2-ODD21) showed up-regulation and downregulation, respectively (**Figure 3A**). In the leaves at night, the transcript of Os2-ODD11 showed up-regulation at all time points, while the transcripts of OsTDC2 and Os2-ODD19 were obviously down-regulated (**Figure 3B**). At vegetative stages, the transcripts of three genes (OsTDC3, OsT5H, and Os2-ODD11) and two genes (OsTDC2 and Os2-ODD19) showed up-regulation and down-regulation, respectively (**Figure 3B**). At reproductive stages, the transcripts of four genes (OsTDC3, OsT5H, Os2- ODD11,−33) and two genes (OsTDC2 and Os2-ODD19) showed up-regulation and down-regulation, respectively (**Figure 3B**). At ripening stages, the transcripts of two genes (Os2-ODD11, −33) and five genes (OsTDC1, −2, Os2-ODD19, OsSNAT, and OsASMT2) showed up-regulation and down-regulation, respectively (**Figure 3B**).

In rice roots at day, transcript of OsASMT2 increased in the entire growth stages in the field (**Figure 4A**). At vegetative stages, the transcripts of eight genes (Os2-ODD11, −19,−33, OsTDC1, −2, OsT5H, OsSNAT, and OsASMT2) showed upregulation (**Figure 4A**). At reproductive stages, the transcripts of two genes (OsASMT1, −2) and four genes (OsTDC1, −3, Os2-ODD11,−21) showed up-regulation and down-regulation, respectively (**Figure 4A**). At ripening stages, the transcripts of two genes (OsASMT1 and −2) and eight genes (OsTDC1, −2,−3, Os2-ODD19, −21,−33, OsT5H and OsSNAT) showed up-regulation and down-regulation, respectively (**Figure 4A**). Interestingly, the transcripts of six genes (OsTDC1, −2, Os2-ODD19,−33, OsT5H, and OsSNAT) were obviously upregulated at vegetative stages, but down-regulated at ripening stages (**Figure 4A**). In the roots at night, transcript of Os2- ODD21 decreased throughout entire growth stages in the field (**Figure 4B**). At vegetative stages, the transcripts of six genes

(OsTDC1, −2, OsT5H, OsSNAT, Os2-ODD33, and OsASMT2) and two genes (Os2-ODD11, −21) showed up-regulation and down-regulation, respectively (**Figure 4B**). Additionally, the transcripts of two genes (OsTDC3 and Os2-ODD21) and five genes (Os2-ODD11, −21,−33, OsTDC2, and OsT5H) showed down-regulation at reproductive and ripening stages, respectively (**Figure 4B**). The transcripts of three genes (OsTDC2, OsT5H, and Os2-ODD33) showed up-regulation at vegetative stages, but down-regulation at ripening stages (**Figure 4B**).

Some genes shared similar expression patterns at day and night at the same stage (**Figures 3** and **4**). The transcripts of Os2-ODD11 and -33 showed up-regulation at reproductive and ripening stages in leaves, the transcripts of six genes (OsTDC1, −2, OsT5H, OsSNAT, OsASMT2, and Os2-ODD33) showed upregulation at vegetative stages in roots, the transcripts of OsTDC1 and Os2-ODD21 shared down-regulation at reproductive and ripening stages of roots. However, there were also some genes shared different expression patterns. For example, the transcript of OsT5H shared down-regulation in day at vegetative stages in leaves, but showed up-regulation in night. Moreover, some genes also shared similar expression patterns in different tissues at the same stage. At day, the transcripts of OsASMT2 and OsASMT1 showed up-regulation at vegetative and ripening stages in leaves and roots, respectively, while the transcript of Os2- ODD21 shared down-regulation at reproductive and ripening stages. On the contrary, some genes shared different expression patterns in different tissues at the same stage. The transcripts of OsT5H showed down-regulation in day at vegetative stages in leaves, but showed up-regulation in roots. The transcripts of Os2-ODD11 and Os2-ODD33 showed up-regulation at ripening stages of leaves, but showed down-regulation at roots. These results suggest that the melatonin-related genes as well as endogenous melatonin may play important and dual roles in rice developmental stages.

### Diurnal and Circadian Transcriptional Profile of Rice Melatonin Synthesis and Catabolic Genes throughout Entire Growth

As shown in **Figure 5A**, Os2-ODD11 expression was induced throughout entire growth stages, and the transcripts of four genes (OsASMT1,−2 OsTDC1 and Os2-ODD33) were induced at most time points of growth stages, while those of Os2-ODD19 and Os2-ODD21 were intermittent. Before reproductive 1 stage, Os2-ODD21 showed significant induction at night. On the contrary, Os2-ODD19 was obviously down-regulated. This result suggested that Os2-ODD19 and Os2-ODD21 may play dual and important roles in the regulation of circadian rhythm. Moreover, the transcripts of four genes (Os2-ODD11, −33, OsASMT1, and −2) and two genes (OsASMT1 and Os2-ODD11) displayed significant up-regulation at four stages (vegetative 3, vegetativereproductive, reproductive 1, and reproductive-ripening stages) and two stages (ripening 1 and 2 stages), respectively. However, OsTDC2, OsTDC3, and OsSNAT expressions were repressed during ripening 1 and 2 stages. Additionally, the transcripts of most genes were induced during reproductive 1 and reproductive-ripening stages. At the last two stages (ripening 1 and 2 stages), most of genes were obviously down-regulated. Interestingly, the transcripts of five genes (OsTDC2, −3, OsT5H, OsSNAT, and OsASMT2) were induced during reproductive 1 and

heatmap (A) and line chart (B). The transcriptional expression data were obtained from RiceXPro (http://ricexpro.dna.affrc.go.jp/). Rice (O. sativa L. ssp. japonica cv. Nipponbare) seeds were sown in germinating boxes. Gene expression profile of rice plants grown under natural field conditions based on microarray analysis of leaf samples at various growth stages encompassing the vegetative, reproductive and ripening stages. Samples corresponding to the uppermost fully expanded leaves were collected in a 48-h period at 2-h intervals at eight different growth stages. The cluster analysis and normalized signal of expression profile for each gene in various tissues is shown as normalized data (log2).

reproductive-ripening stages, but exhibited down-regulation at ripening 1 and 2 stages. As shown in **Figure 5B**, the expression of Os2-ODD19 and Os2-ODD21 displayed obviously regular change throughout entire growth.

In roots, the transcript of Os2-ODD11 was up-regulated at all time-points, while those of OsTDC3, Os2-ODD21, and OsTDC2 were down-regulated (**Figure 6**). The transcriptional profile of Os2-ODD19 was intermittent, which was consistent with the result in leaves. During 15–17 DAT, the expressions of Os2- ODD21, OsTDC2, and OsASMT2 were repressed. The transcript of Os2-ODD11 was significantly induced during 43–45 DAT, whereas the transcripts of OsTDC1, −2, −3, OsT5H, and Os2- ODD21 were repressed.

### The Transcriptional Profile of Rice Melatonin Synthesis and Catabolic Genes in Response to Plant Hormones

In response to ABA and JA treatments, the transcripts of OsT5H, OsTDC2, −3, and Os2-ODD19 displayed significantly up-regulation in root or shoot (**Figure 7**). The transcript of Os2-ODD11 was significant up-regulated after IAA, BL, and JA treatments in root, but was strongly down-regulated after IAA, BL, and tZ treatments in shoot. The transcripts of OsTDC1 and OsASMT1 showed up-regulation after ABA, GA3, IAA, BL, and tZ treatments in shoot, while OsTDC3 expression was induced after ABA, IAA, BL, tZ, and JA treatments in shoot (**Figure 7**).

Generally, melatonin-related genes showed different expression profiles in root or shoot tissues for the same treatment (**Figure 7**). The transcripts of OsTDC1, -3, Os2- ODD33, and OsASMT1 were up-regulated in the shoots, but were not significantly regulated or down-regulated in the roots after ABA and tZ treatments. The expression of Os2-ODD11 was increased in roots, but was decreased in shoots after ABA and BL treatments. Although some melatonin-related genes were from the same family, they exhibited different responses to plant hormones treatments, such as OsTDC2 and OsTDC3 in roots, Os2-ODD11 and Os2-ODD21 in roots, OsASMT1 and OsASMT2 in shoots. Thus, the transcriptional response of melatonin-related genes to plant hormones treatments in roots and shoots may provide new insight into crosstalk between melatonin and plant hormones, as well as mechanism underlying melatonin-mediated signaling in rice.

### Gene Expression Profile in Whole Leaf of Rice Melatonin Synthesis and Catabolic Genes Inoculated with Pathogen Infection

Because melatonin plays important roles in response to pathogen infection (Yin et al., 2013; Lee et al., 2014, 2015; Reiter et al., 2015; Shi et al., 2015a; Zhao et al., 2015), so we analyzed the expression profile of rice melatonin synthesis and catabolic genes in response to pathogen inoculation to identify the candidate genes for further analysis.

After inoculation with Xoo, OsASMT2 expression was induced during almost all the time-points, while OsASMT1 transcript was decreased at these time points (**Figure 8**). The transcripts of OsT5H and Os2-ODD19 were induced during 1 to 6 dpi and 6 to 12 hpi, respectively (**Figure 8**). OsTDC1, −2, −3, Os2-ODD11, −21,−33, and OsSNAT, expressions were decreased during all the treated time points (**Figure 8**). Interestingly, the transcripts of OsASMT1 and OsASMT2 were increased after 6 dpi of wildtype strain T7114R, but were decreased after 6 dpi of 1ahrcV(III) strain. T3S is essential for XooT7174R conferred plant disease, and 1ahrcV(III) resulted in less plant disease in rice leaves (Sato et al., 2011b, 2013). Thus, the results indicate the possible role of OsASMT1 and OsASMT2 in immune response to Xoo.

After inoculation with the blast fungus (M. oryzae), the transcript levels of three genes (OsT5H, OsASMT1, and Os2- ODD11) showed up-regulation at most treated time points, while those of four genes (OsTDC1, −3, OsSNAT, and Os2-ODD21) showed down-regulation (**Supplementary Figure S1**). In response to inoculation with the three strains of M. oryzae harboring AVR-Pia and AVR-Pish (Pia/Pish × P91- 15B, Pish × P91-15B, and Pish × Kyu77-07A), the transcript levels of OsASMT1, Os2-ODD11, −19, −33, and OsT5H showed significant up-regulation, while those of OsTDC3, OsSNAT, and Os2-ODD21 were seriously down-regulated at all time points (**Supplementary Figure S1**). OsASMT2 expression showed no obvious trends in response to inoculation with the two former, while was induced at all the treated time-points after inoculation with Pish × Kyu77-07A strain (**Supplementary Figure S1**). In response to inoculation with the pish mutant of M. oryzae (1Pish × Kyu77-07A), the expressions of six genes (OsTDC1, OsT5H, OsASMT1, −2, Os2-ODD11, and 19) were significantly induced at all time points, while four genes (OsTDC2, −3, OsSNAT, and Os2-ODD33) were seriously downregulated (**Supplementary Figure S1**). Os2-ODD21 expression was strongly repressed at 2 dpi. Notably, Os2-ODD33 expression was induced at all the treated time-points after inoculation with the three strains of M. oryzae harboring AVR-Pia and AVR-Pish, but was repressed in response to inoculation with 1Pish × Kyu77-07A strain (**Supplementary Figure S1**).

### The Transcriptional Profile of Rice Melatonin Synthesis and Catabolic Genes in Response to Abiotic Stress Treatments

Melatonin is widely involved in plant stress responses (Shi et al., 2015b,d). Thus, investigation of the expression profiles of melatonin-related genes of rice in response to various abiotic stresses is needed. After drought treatment, the transcripts of OsTDC1, OsASMT1, and Os2-ODD19 were found to be upregulated between 1.2 and 1.8-folds, whereas those of OsTDC2, −3, Os2-ODD11, −21, OsT5H, and OsSNAT were strongly repressed in comparison to the control (**Figure 9**). After salt treatment, the transcripts of OsTDC1, −3, OsASMT1, and Os2-ODD19 were increased between 1.3 and 4.5-fold, whereas NaCl strongly repressed OsTDC2, OsT5H, OsSNAT, and Os2- ODD11 expressions (**Figure 9**). After 4◦C treatment, OsTDC1 and OsASMT1 transcript levels were slightly increased, whereas the expressions of six genes (Os2-ODD11, −19, −21, OsTDC3, OsT5H, and OsSNAT) were obviously down-regulated (**Figure 9**).

profile for each gene in various tissues is shown as normalized data (log2).

FIGURE 7 | The transcriptional profile of rice melatonin synthesis and catabolic genes in root (A) and shoot (B) in response to plant hormones. The transcriptional expression data were obtained from RiceXPro (http://ricexpro.dna.affrc.go.jp/). The cluster analysis of expression profile for each gene in various tissues is shown as normalized data (log2) in relative to the 0 h of treatment which was set as 0.

profile for each gene in various tissues is shown as normalized data (log2).

Notably, the transcript of OsTDC3 was significantly increased after salt treatment, but was seriously decreased in response to drought and cold treatments (**Figure 9**). The result indicated that OsTDC1 and OsTDC3 may be involved in salt stress response.

### DISCUSSION

To our knowledge, this is the first study systematically analyzed the diurnal and circadian transcriptional profile of melatonin synthesis and catabolic genes throughout the entire growth stages in rice. Our study confirmed that OsTDCs showed higher expression level at reproductive 1 and reproductive ripening stages then other stages (**Figure 5**). In cherry fruit (Prunus avium), the expression level of PaTDC was positively correlated to melatonin concentration throughout the entire period, and showed regularly circadian rhythm during a 24 h period with two peaks at 5:00 and 14:00 (Zhao et al., 2013), indicating that the melatonin concentration was higher in that two stages then others. The expression pattern of OsTDC3 also showed circadian

rhythm on 16 DAT (**Figure 6**). Similarly, this expression pattern was also existed at Os2-ODD19 during the entire development (**Figure 5**). These results suggested that OsTDC3 and Os2-ODD19 may be involved in modulating endogenous during the entire development in rice.

Melatonin is widely involved in plant growth and development, as well as stress responses (Bajwa et al., 2014; Meng et al., 2014; Wang et al., 2014; Zuo et al., 2014; Liang et al., 2015; Shi et al., 2015b). In apple (Malus prunifolia), the transcripts of melatonin synthesis genes (MdTDC1, MdAANAT2, MdT5H4, and MdASMT1) was induced after drought treatment (Li et al., 2014). The concentration of melatonin was increased in barley roots and lupin after cold, drought and salt treatments (Arnao and Hernández-Ruiz, 2009, 2013). Consistently, OsTDC1 and OsASMT1 transcript levels were increased after drought, salt and cold treatments (**Figure 9**), indicating their possible involvement in abiotic stress response.

Recently, melatonin was shown to function as positive modulator against plant pathogen infection (Yin et al., 2013; Lee et al., 2014, 2015; Reiter et al., 2015; Shi et al., 2015a; Zhao et al., 2015). Exogenous application of melatonin improved apple resistance to D. mali, the pathogen of Marssonina apple blotch (Yin et al., 2013), enhanced Arabidopsis and tobacco defense against P. syringae pv. tomato DC3000 (Lee et al., 2014). Moreover, Arabidopsis snat knockout mutants exhibited increased susceptibility to the avirulent pathogen P. syringae pv. tomato DC3000 with decreased SA levels and reduced defense genes expression compared with wild-type (Lee et al., 2015). However, whether OsSNAT also involves in the regulation of plant immunity remains unknown. Treatments with melatonin significantly enhanced antioxidant protection in rice, suggesting that melatonin plays a major role in regulating pathogen infection (Liang et al., 2015). In response to bacterial pathogen infection, some genes showed similar expression profiles (**Figure 8**). The differential response of melatoninrelated genes to pathogen infection in different kinds implied the dual mechanisms underlying melatonin-related genes mediated pathogens responses.

Previous studies have revealed that melatonin had significant effect in regulating ABA and GA4 in plant response to salinity and drought stress (Li et al., 2014; Zhang H.J. et al., 2014). Additionally, melatonin shared the common substrate (tryptophan) with IAA (Arnao and Hernández-Ruiz, 2014), and AXR3/IAA17 is involved in Arabidopsis melatonin signaling underlying senescence (Shi et al., 2015c). Thus, genome-wide expression analysis of melatonin-related genes in response to plant hormones may provide new insight into crosstalk between melatonin and plant hormones. Plant hormones, such as ABA, SA, and GA, related with most of the plant physiological responses, including water logging, drought and salt stress responses (Yang et al., 2004; Kim et al., 2011, 2014; Shimamura et al., 2014). Melatonin is also a regulator of ABA and GA4 in plant response to salinity and drought stress (Li et al., 2014; Zhang H.J. et al., 2014). Moreover, the expression of OsASMT2 and OsASMT3 were induced after ABA and JA treatments at 1-month-old detached rice leaves, while were down-regulated in response to ethephone, zeatin, and SA stress. OsASMT1 also showed up-regulation upon ABA and JA stress treatment, but

did not display obvious trends during ethephone, zeatin, and SA treatments (Park et al., 2013a). In this study, the transcripts of four genes (Os2-ODD11,−19, OsASMT2, and OsT5H) and five genes (OsASMT1, −2, OsT5H, Os2-ODD19, and OsTDC3) were increased in response to ABA stress during 3 h to 6 h treatment in root and during 3 to 12 h in shoot, respectively. Under JA stress, the transcripts of three genes (OsASMT2, OsT5H, and Os2- ODD11) and three genes (OsASMT1, −2, and Os2-ODD11) were significantly increased during 0.25 to 6 h in root and 3 h to 12 h in shoot, respectively. Thus, different transcriptional responses of melatonin-related genes in hormone specific manners, suggested the dual role and crosstalk between melatonin and various hormones.

It is widely known that ABA is the most important regulator of abiotic stress (Kim et al., 2011, 2014; Shimamura et al., 2014), and JA serves as the major defense hormone that are associated with pathogen infection (Van der Ent et al., 2009; Ballaré, 2011; Xie et al., 2011; Yamada et al., 2012; Yang et al., 2013; Campos et al., 2014). More recently, the crucial role of ABA in virulence of rice blast fungus M. oryzae is confirmed (Spence et al., 2015), the involvement of JA in abiotic stress response is also largely confirmed (Riemann et al., 2015; Wasternack and Strnad, 2016). We highlight the commonly regulation of rice melatonin-related genes by ABA, JA, pathogen infection and various abiotic stresses (**Figures 7–9**), indicating the possible role of these genes in multiple stress responses and underlying crosstalks of plant hormones, especially ABA and JA. Weeda et al. (2014) identified 1308 differentially expressed genes (566 up-regulated genes and 742 down-regulated genes) exhibiting at least of twofold change by exogenous melatonin treatment in Arabidopsis, and many of them are enriched in plant hormone signaling. These differentially expressed genes include 52 genes in auxin signaling, 50 genes in ABA signaling, 67 genes in JA pathway, and 42 genes in ET pathway. Our studies together with the data of Weeda et al. (2014) further indicate the interaction among melatonin, ABA and JA pathways.

We have to pointed out the possible limitation of this study, since the different changes in expression levels of the various genes do not always explain in a simple way why melatonin concentrations increase or decrease under the different conditions. On one hand, there may be difference between expression level and enzyme activity, such as the posttranslational regulation of AANAT in primates via phosphorylation/dephosphorylation and association/dissociation of a 14-3-3 protein, which is decisive for the melatonin rhythm in those organisms (Ganguly et al., 2001, 2005). On the other hand, incomplete knowledge of rate-limiting enzymes or isoenzymes may also lead to the difference. Further studies by other methods may give more clues.

Taken together, the expression patterns of 11 melatonin related genes from rice were synthetically analyzed at different periods and after different treatments in this study. These information may provide abundant resources for functional characterization of melatonin related genes. The differential expression patterns of melatonin related genes in different tissues throughout entire growth stages and stress responses will be useful to investigate in vivo role of specific gene in rice development and circadian rhythm. Thus, this study will contribute to better understand the melatonin biosynthesis and catabolic pathway as well as their association with development and stress responses in rice. Further functional analysis of identified candidate genes with potential involvement in circadian rhythm and stress responses will give shed more lights in melatonin-mediated signaling as well as underlying molecular mechanism.

### AUTHOR CONTRIBUTIONS

HS conceived and directed this study, analyzed the data, wrote, and revised the manuscript; YW and HZ performed the experiments, analyzed the data, wrote, and revised the manuscript; WH and LC analyzed the data and revised the manuscript; CH provided suggestions and revised the manuscript. All authors approved the manuscript and the version to be published, and agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

### ACKNOWLEDGMENTS

This research was supported by the National Natural Science Foundation of China (No.31570249 and No.31260418), the Special Fund for Agro-scientific Research in the Public Interest of China (No.201403075), the Department of Science and Technology of Hainan Province (No.ZDZX2013023), the scientific research foundation of higher education in Hainan province (the education curriculum reform program, title: Research on the overall optimization of the curriculum system and teaching content of Hainan province excellent course-gene engineering, No.Hnjg2016-10), a central financial support to enhance the comprehensive strength of the central and western colleges and universities, the startup funding and the scientific research foundation of Hainan University (No.kyqd1531) and the education curriculum reform program of Hainan University (No.hdjy1601) to HS.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.00676

TABLE S1 | The detailed data used for Figures 3–7 and Figure 9. The original data were obtained from RiceXPro (http://ricexpro.dna.affrc.go.jp/) and Rice eFP Browser (http://bar.utoronto.ca/efprice/cgi-bin/efpWeb.cgi).

FIGURE S1 | Gene expression profile in whole leaf of rice melatonin synthetic and degradated genes inoculated with the blast fungus, M. oryzae. The transcriptional expression data were obtained from RiceXPro (http://ricexpro.dna.affrc.go.jp/). Three lines of rice cultivar Nipponbare carrying the blast resistance genes (Pia, Pish) were inoculated with two strains of M. oryzae harboring AVR-Pia and AVR-Pish. Rice seedlings at the 4-leaf stage were inoculated with the conidial suspension of M. oryzae and the leaves were harvested at 1, 2, 3, and 5 dpi. The expression profile for each gene in various tissues is shown as normalized data (log2).

### REFERENCES

fpls-07-00676 May 14, 2016 Time: 12:18 # 13

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Wei, Zeng, Hu, Chen, He and Shi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Simple, Rapid Method for Determination of Melatonin in Plant Tissues by UPLC Coupled with High Resolution Orbitrap Mass Spectrometry

#### Tiantian Ye1 , Yan-Hong Hao1 , Lei Yu1 , Haitao Shi2 , Russel J. Reiter3 and Yu-Qi Feng1 \*

1 Key Laboratory of Analytical Chemistry for Biology and Medicine, Department of Chemistry, Ministry of Education, Wuhan University, Wuhan, China, 2 Hainan Key Laboratory for Sustainable Utilization of Tropical Bioresources, College of Agriculture, Hainan University, Haikou, China, 3 Department of Cellular and Structural Biology, The University of Texas Health Science Center, San Antonio, TX, USA

### Edited by:

Thomas Vogt, Leibniz Institute of Plant Biochemistry, Germany

### Reviewed by:

Mwafaq Ibdah, Agricultural Research Organization, Israel Kexuan Tang, Shanghai Jiao Tong University, China

> \*Correspondence: Yu-Qi Feng<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Plant Metabolism and Chemodiversity, a section of the journal Frontiers in Plant Science

Received: 02 November 2016 Accepted: 12 January 2017 Published: 25 January 2017

### Citation:

Ye T, Hao Y-H, Yu L, Shi H, Reiter RJ and Feng Y-Q (2017) A Simple, Rapid Method for Determination of Melatonin in Plant Tissues by UPLC Coupled with High Resolution Orbitrap Mass Spectrometry. Front. Plant Sci. 8:64. doi: 10.3389/fpls.2017.00064 Melatonin (MLT) was involved in regulating various stages of plant growth and development. However, due to the low concentration and complex matrixes of plant, the analysis of MLT is a challenging task. In this study, we developed a rapid and efficient method with simplified sample preparation by employing UPLC coupled with a high resolution Orbitrap mass spectrometry, and stable isotope-labeled MLT (MLT-d4) was first used as internal standard in the developed analytical method. In the developed method, we used one-step liquid–liquid extraction to purify the crude extracts both from shoot and root of rice for the analysis, which remarkably simplify the sample preparation process. The method exhibits high specificity and recovery yield (>96.4%). Good linearities were obtained for MLT ranging from 0.01 to 20 ng/ mL with determination coefficient (R 2 ) of 0.9991. The limit of detection for MLT was 0.03 pg. Reproducibility of the method was evaluated by intra-day and inter-day measurements and the results showed that relative standard deviations were less than 7.2%. Moreover, MLT quantification was accomplished by using only 100 mg fresh plant tissues. Additionally, the established method was successfully applied to investigate the spatiotemporal distributions of MLT in rice under cadmium (Cd) stress condition. We found that the content of MLT in shoot and root of rice increased under Cd stress, suggesting that MLT would play a crucial role in modulating the responses to Cd stress in different plant tissues.

Keywords: LC-MS, melatonin, Orbitrap, rice, stable isotope

## INTRODUCTION

Studies on melatonin (N-acetyl-5-methoxytryptamine, MLT) in plants have attracted great attention of researchers in recent years. MLT is an indoleamine that has been found to widely distribute in the animal kingdom (Arnao and Hernández-Ruiz, 2015). As an animal hormone, especially a neurohormone (Reiter, 1991), it plays a vital role in vertebrates. Since the first discovery and detection of MLT in plants by Dubbels et al. (1995) and Hattori et al. (1995), successive studies have been invested in its possible physiological functions in plants. In the past 20 years, MLT has been found in more than 140 plant species

(Hardeland et al., 2007; Zohar et al., 2011; Gomez et al., 2013; Arnao, 2014; Feng et al., 2014). Moreover, it was reported to be involved in almost all of the plant growth and development stages, from seeds germination to the senescence (Hardeland et al., 2007; Paredes et al., 2009; Arnao and Hernández-Ruiz, 2015; Reiter et al., 2015; Hardeland, 2016). Additionally, many studies have shown that MLT could effectively protect plant tissues from the damage caused by a variety of abiotic and biotic stresses, such as drought stress (Shi et al., 2015b,c), salt stress (Shi et al., 2015b,c), cold stress (Bajwa et al., 2014; Shi et al., 2015b,c), heat stress (Byeon and Back, 2014; Shi et al., 2015d), and pathogen infection (Yin et al., 2013; Shi et al., 2015a).

Accurate and sensitive determination of endogenous MLT enabled us to better understanding its physiological functions, biological synthesis and metabolic pathways, and the regulatory networks. Various quantitative approaches have been developed to analyze MLT content in plants. Immunoassays such as radioimmunoassay (RIA) and enzyme-linked immunosorbent assay (ELISA) have long been applied to MLT quantification (Dubbels et al., 1995; Hattori et al., 1995; Pape and Lüning, 2006; Gomez et al., 2013). ELISA kits are commercially available, along with its ease of operation make it a convenient and common tool for MLT analysis (Pape and Lüning, 2006). However, the cross-reactivity of antibodies is inevitable, which leads to reduced specificity and accuracy (Huang and Mazza, 2011a; Gomez et al., 2013). Gas chromatography coupled to mass spectrometry (GC-MS) has also been used to the qualitative and quantitative analysis of MLT, although it has improved sensitivity and specificity, the requirement for sample volatile derivatization presents a disadvantage for its application to MLT determination (Badria, 2002; Ragab et al., 2006). The sensitivity of capillary electrophoresis (CE) methods has also been reported to be similar to the high performance liquid chromatography (HPLC) approaches (Gomez et al., 2015). HPLC is the most commonly used technique for the analysis of MLT. In this case, different detectors have been used, such as fluorescence detector (FD) (Lavizzari et al., 2006; Pape and Lüning, 2006; Arnao and Hernández-Ruiz, 2007; Padumanonda et al., 2014), electrochemical detector (ECD) (Hernández-Ruiz et al., 2004, 2005; Reiter et al., 2005; Zettersten et al., 2009), and UV (Hosseinian et al., 2008; Ly et al., 2008). However, compared to the MS detector, these detectors exposed their shortcomings of low sensitivity and specificity. LC coupled with MS has been proved as a powerful method for MLT analysis because of its high sensitivity, reproducibility and accuracy (Cao et al., 2006; Hernández-Ruiz and Arnao, 2008; Huang and Mazza, 2011b; Gomez et al., 2012, 2013).

Unfortunately, it is unpractical to directly analyze endogenous MLT in plant tissues by LC-MS due to its low concentration and the complicated matrix of plant extracts. To overcome the difficulty, different sample purification techniques such as liquid–liquid extraction (LLE) (Arnao and Hernández-Ruiz, 2009; Huang and Mazza, 2011a,b; Byeon and Back, 2014), solid phase extraction (SPE) (Huang and Mazza, 2011a; Gomez et al., 2015), microextraction by packed sorbent (MEPS) (Mercolini et al., 2012), HPLC purification (Pape and Lüning, 2006) were commonly employed to enrich MLT from plant tissues. As the classical techniques, LLE and SPE are the most used methods for purification of MLT (Huang and Mazza, 2011a,b; Gomez et al., 2015). At the previous stage, a large volume of solvent and a large amount of plant tissues were needed when LLE was utilized. Thanks to the development of detection techniques, the extraction solvent volume for LLE and the plant amount needed could be decreased, and the LLE merits of fast and easy operation could be kept.

Moreover, matrix effects caused by endogenous interferents often occur during the analysis of low abundant metabolites in the complicated matrix of plant extracts using LC-MS. Nevertheless, these effects can be compensated by using isotope

of dilution techniques (Barkawi et al., 2010; Porfírio et al., 2016). Stable isotope-labeled compound(s) is (are) essential internal standard(s) for their physical and chemical similarities with the original analyte(s), providing correction for loss analyte during the experiment procedures and matrix effects during ionization (Barkawi et al., 2010; Porfírio et al., 2016). However, isotope dilution technique has not been reported to be applied to the analysis of MLT content in plants.

Additionally, a highly sensitive and selective MS method was crucial to identify analytes of intrest from complicated matrix of plant. The higher specificity of the technique also leads to a higher sensitivity due to the reduced background noise. Multi reaction monitoring (MRM) (Gomez et al., 2012, 2013; Kocadali et al., 2014) and selective reaction monitoring (SRM) modes (Huang and Mazza, 2011b) have been used previously to the quantification of MLT. And they have been reported to be powerful tools to improve the sensitivity and specificity of analytes. Ions at m/z 216 and 174 were chosen as MS/MS signatures for MLT and have proven particularly useful for quantification (Huang and Mazza, 2011b; Gomez et al., 2012; Gomez et al., 2013; Kocadali et al., 2014). However, the low mass

solvent on the extraction efficiencies and matrix effects of MLT. (A) Extraction recoveries of ACN and MeOH solvent for MLT. (B) MLT contents of rice shoot using ACN and MeOH as extraction solvent, respectively. (C) The matrix effects of ACN extract and MeOH extract for MLT analysis. The calculation of matrix factor (MF): the peak area response for MLT-d4 in the presence of the plant matrix before injection was compared with peak area response in the absence of matrix. Asterisk symbols ∗P < 0.05 and ∗∗P < 0.01 (t-test).

resolution used in these studies could result in false positives, specifically in case of trace metabolites like MLT. Therefore, highresolution LC-MS is of utmost importance in identification of trace metabolites like MLT or IAA allowing for high sensitivity and appropriate specificity even from crude plant extracts (Yu et al., 2014).

Herein, we developed a rapid, sensitive, and efficient method for the analysis of MLT content by employing a high resolution Orbitrap MS. It enables us to minimize the steps of sample preparation. For both shoot and root of rice, just one-step LLE was added to purify the crude extracts for the analysis. Furthermore, to correct the loss analyte and matrix effects, stable isotope-labeled MLT (MLT-d4) was added as internal standard to plant samples prior to extraction. Using the developed method, we successfully investigated the dynamic distributions of MLT in response to Cd stress in shoot and root of rice.

### MATERIALS AND METHODS

### Chemicals and Reagents

fpls-08-00064 January 23, 2017 Time: 15:30 # 4

Melatonin standard was purchased from Sigma Chemical (St. Louis, MO, USA). Stable isotope-labeled standard, [2H4] MLT, was purchased from C/D/N Isotopes Inc. (Pointe-Claire, QC, Canada), Chromatography grade acetonitrile (ACN) and methanol (MeOH) were obtained from TEDIA Co. (Fairfield, OH, USA). Ultra-pure water used in the study was purified with Milli-Q system (Milford, MA, USA). Hydrochloric acid (HCl) was purchased from Sinopharm Chemical Reagent (Shanghai, China).

### Plant Materials

Rice (Oryza sativa ssp. japonica cv. Nipponbare) seeds were germinated and then transplant under hydroponic conditions in Hoagland's nutrient solution in a growth chamber with 70–80% relative humidity under photoperiods of 16 h light (28◦C)/8 h dark (25◦C). After 10 days growing under normal conditions, the seedlings were subjected to Cd stress with a series concentration of CdNO3 (0, 100, 200, and 500 µM). Shoot and root of rice were harvested after 5, 10, 15 days stresses separately, weighted and then stored at –80◦C immediately after freezing in liquid nitrogen.

### Preparation of Plant Samples

Rice samples (100 mg F.W.) were frozen in liquid nitrogen and finely ground, followed by extraction with 1 mL ACN at 4◦C for 12 h in dark. [2H4] MLT (1.0 ng/g) were added to the samples as internal standards (IS). After centrifugation at 12,000 rpm and 4 ◦C for 10 min, the supernatants were sequentially evaporated under mild nitrogen stream. The evaporated samples followed by re-dissolving in 100 µL 0.1 M HCl, and then extracted with ether (1 mL). The ether phase was collected, dried under nitrogen gas and reconstituted in 100 µL H2O/MeOH (80/20, v/v) for further analysis."

### Instruments and Analytical Conditions

Analysis of MLT was performed on an UltiMate 3000 UHPLC System (Thermo-Dionex) equipped with LTQ Orbitrap MS (Thermo-Fisher Scientific, Waltham, MA, USA). The separation of MLT was achieved on an Acquity UPLC BEH phenyl column (2.1 mm × 100 mm, 1.7 µm, Waters) with a flowrate of 0.2 mL/min. The column oven temperature was set at 35◦C. Water (solvent A) and acetonitrile (solvent B) were employed as mobile phase. The gradient was started at an initial composition of 95% A and 5% B, 2–20 min, 40% B, then 23 min linear gradient to 90% B, held for 5 min. A return to the initial conditions was accomplished by a 2 min gradient to 95% B, it was held for 10 min, total chromatographic run time was 40 min.

TABLE 1 | Calibration characteristics of melatonin (MLT) standard spiked in solution of rice seedling samples.

extract of 100 mg rice shoot (A–D) and root (E–H) samples with or without purification of LLE.

The MS was set to acquire full MS scan in positive ion mode with a mass range of m/z 150–300 at a resolution of 120,000. Ion source conditions were as follows: heater temperature, 300 ◦C; capillary temperature, 350 ◦C; sheath gas flow, 35 arbitrary; auxiliary gas flow, 15 arbitrary; spray voltage, 3.5 kV; S-lens RF level, 60%. Data acquisition and analysis was performed using Xcalibur v3.0.63 (Thermo Fisher Scientific) and SIEVE v2.2 (ThermoFisher Scientific, USA).

### Determination of Electrolyte Leakage (EL)

Electrolyte leakage (EL) was determined according to the previously reported method (Ye et al., 2015). In brief, the initial conductivity was determined after gently shaking at room temperature for 6 h at 150 rpm using a conductivity meter (Leici-DDS-307A, Shanghai, China). The fully releasing conductivity was measured after boiling at 121 ◦C for 20 min using previous samples. The percentage of EL was determined as the ratio of the initial conductivity to fully releasing conductivity.

## Statistical Analysis

The results were shown as means ± standard errors ( n = 3), SPSS 13.0 software was used for statistical analysis. Analyses of variance (ANOVA) for variables from measurements were used for testing the species and treatment differences according to Duncan's method. Different letters above the columns in each figure indicate significant differences at P < 0.05.

## RESULTS AND DISCUSSION

## General Design of Strategy

Generally, each additional step would result in loss of analyte to some extent, especially for MLT with the antioxidant capability and light-sensitive characteristic. So, minimizing the number of sample preparation steps is essential. Moreover, high-resolution and accurate MS would benefit for simplifying procedures for sample preparation procedure and reducing analysis time, so an UPLC-high resolution Orbitrap MS was applied to the analysis of MLT in our study. Additionally, stable isotopelabeled MLT (MLT-d4) was added as internal standard, providing correction for loss analyte during the experiment procedures and matrix effects during ionization. The whole schematic diagram for the analysis of MLT in this study was shown in **Figure 1** .

### Comparison of Different Extraction Solvent

Melatonin is amphipathic substances, with antioxidant capabilities and poor ionization efficiency (Huang and Mazza, 2011a). And concentration of MLT in complicated plant matrix is fairly low, whereas numerous of other substances with far greater amounts are present. Therefore, it is important to choose an extraction method that yields good recovery. To achieve high extraction efficiency, two frequently used organic solvents, ACN and MeOH were chosen for the optimization of MLT extraction

from plant tissues. The retention time and accurate mass of MLT standard and MLT-d4 standard were shown in **Figure 2**. Our result showed the extraction efficiencies of MLT by ACN and MeOH were similar, and there is no significant different of MLT content in shoot (MLT signal was severely suppressed in crude extract of root in both solvents, and it will be discussed in next section) of rice using ACN (1.55 ng/g) or MeOH (1.47 ng/g), and the both solvents exhibited good extraction recoveries (>80%) (**Figure 3**).

In plant samples, matrix often play negative roles in MS signals. Due to the complicated matrix of plant extract, the abundant interferents from the plant extract might depress the MS response of MLT. So, the matrix effects were further assessed between ACN extract and MeOH extract. The MS response (peak area) for MLT-d4 in the presence of the plant matrix was compared to peak area in the absence of matrix to calculate matrix factor (MF) according to Yu et al. (2016). For shoot samples, MF from ACN extract was slightly higher than MeOH extract. But for root sample, a strong matrix effect was observed for MeOH extract (MF, 30%), while approximately more than two times lower than ACN extract (MF, 76%) (**Figure 3**). Therefore, ACN was selected as appropriate solvent for extracting MLT from rice tissues.

### Optimization of Sample Preparation Procedure

In the previous section, we mentioned that MLT signal was strongly suppressed in crude extract of root. To solve this problem, further purification step was needed. Taking advantage of the analyte solubility differences in different solvents, LLE is an effective way to purify the MLT in plants (Chen et al., 2012). Furthermore, LLE can improve the ionization efficiency by decreasing the influence of impurities (Chen et al., 2012). The evaporated sample followed by re-dissolving in 100 µL 0.1 M HCl, and then extracted with ether (1 mL). The results showed that with the purification of LLE, the signal intensity of MLT in root samples was significantly increased, while the intensity of MLT was also increased slightly in shoot samples (**Figure 4**). Therefore, LLE was adopted to eliminate the polar impurities and improve the ionization efficiency of MLT in rice.

### Method Validation

To evaluate the quantitative aspect of the developed method, calibration curve of MLT was established. The MLT standard (0.01–20 ng/mL) was spiked in shoot and root samples (100 mg F.W.). The internal standard MLT-d4 was added to the sample prior to extraction step. The calibration curves were plotted by the peak area ratio (MLT labeled with MLTd4) against the MLT concentrations. As shown in **Table 1**, good linearity was obtained for MLT with coefficient of determination (R 2 ) of 0.9991. The LODs and the LOQs were calculated as the concentration of the analyte at a Signal/Noise (S/N) ratio of 3 and 10, respectively. The results showed that LOD and LOQ for MLT were 0.03 and 0.11 pg, respectively.

The accuracy and precision of the developed method were further assessed, respectively, by the recoveries and by the relative standard deviations (RSDs) of intra- and inter-day. Both recoveries and intra- and inter-day RSDs were calculated with MLT standard spiked in rice samples at three different concentrations. The intra-day variation was evaluated by repeating the process for five times within one day and the interday variation was investigated on five successive days. For each concentration, five replicated measurements were performed in both shoot and root samples. The relative recoveries were calculated by comparing the peak area ratios of MLT from the spiked rice samples to those from the standard solutions. As shown in **Table 2**, good precisions were obtained with RSD values below 7.2%, and the relative recoveries were in the range of 96.4 to 105.1%, indicating the good reproducibility and accuracy of developed method.

Furthermore, to assess the general utility of proposed method for determination of MLT in other plant species or different tissues of plant, the method was used to detect MLT in whole seedlings and leaves of Arabidopsis (Col-0), rice leaves (Oryza sativa ssp. japonica cv. Nipponbare), cotton leaves and fibers (Gossypium hirsutum), rape leaves (Brassica napus L. cv. Zhongshuang 11). As shown in **Table 3**, the content of MLT varies in different plant species and tissues. The MLT content in leaves of one-month-old Arabidopsis and Brassica napus was similar. In rice leaves (10 days after transplant under hydroponic conditions), the MLT reached 0.58 ng/g F.W., which was about threefolds higher than in Arabidopsis leaves. The MLT content in Gossypium sp of maturity date was much higher than rice. The relative recoveries from different plant tissues were found to be in the range of 84.3 to 104.2%, indicating the universal applicability of the developed method.

Compared with previous LC-MS methods, our proposed method showed the following improvements. Firstly, high resolution MS was used for the detection of MLT. Accurate m/z

(5 ppm) of the target analytes was moniterd with high selectivity. The high specificity of the m/z value lead to the decreased matrix interference and improved sensitivity of analytes by reducing the background noise. Secondly, we employed a simple and effective LLE procedure to purify MLT in plant tissues, and the matrix interference got further decreased. Thirdly, the stable isotope dilution strategy was first used for absolute quantification of MLT in plant tissues, and worked well to correct the loss of analyte

fpls-08-00064 January 23, 2017 Time: 15:30 # 7

during the sample pretreatment procedures and signal variation in MS detection.

### MLT Distribution under Cd Stress in Shoot and Root of Rice

Cadmium is one of most toxic pollutions in the present. Unlike the other heavy metal pollutions, Cd is a non-essential nutrient, and seriously affects plant growth and development. Recently, Cd has been reported to regulate MLT content in rice leaves by coordinating the synthesis and degradation genes of MLT (Byeon et al., 2015). And then Hasan et al. found, exogenous application of MLT could mitigate phytotoxicity induced by Cd stress in Solanum lycopersicum L., by modulation of phytochelatins biosynthesis, vacuolar sequestration, and antioxidant potential (Hasan et al., 2015). However, our knowledge of MLT for regulating response of plant to Cd stress is still fragmentary (Arnao and Hernández-Ruiz, 2014).

In this study, we investigate the distributions of MLT both in shoot and root of rice under Cd stress condition using the established method. 10-day-old rice seedlings were subjected to Cd stress with a series concentration of CdNO3 (0, 100, 200, and 500 µM). When exposed to 15 day treatments, significant differences in the growth of rice were observed, suggesting that Cd stress caused serious damages (growth inhibition, wilted leaves) on rice (**Figure 5**). The shoot weight only reached 63.5 mg at 15 day under Cd condition of 500 µM, which were only 13% of control seedlings, while the growth of root was almost completely inhibited with the Cd treatment of 200 and 500 µM (**Figure 5**). Additionally, Cd stress progressively increased the EL with the duration of the treatment. The EL reached 91% at 15 day under Cd treatment of 500 µM (**Figure 5**). These result indicated Cd stress caused severe cell membrane damages and severely inhibited rice growth. Nevertheless, the MLT content in rice shoot was increased with the increasing of Cd concentration and the stress time. After 15 day treatment of 500 µM Cd, the MLT content in rice shoot was up to 21.0 ng/g, which was about 10-folds higher than control (**Figure 6**). Additionally, the MLT levels in rice root under Cd stress were detected as well (**Figure 6**). The results showed that MLT content was also gradually increased over time after Cd treatment in root. Notably, higher concentration of Cd treatment induced higher level of MLT. Under Cd condition of 500 µM at 15 day, the MLT content in root was up to 3.9 ng/g, which was over three times higher than control. Our results suggested that MLT might play a critical role in rice tolerance to Cd stress both in shoot and root tissues. The MLT content was also induced by Cd treatment in root, indicating MLT might be involved in preventing translocation of Cd from root to shoot, and might modulate the responses to Cd stress through the different mechanisms in shoot and root. Further studies are still needed to provide the evidence in support of the involvement of melatonin in Cd tolerance and detoxification, and to illuminate signaling pathway in different tissues of plants.

### CONCLUSION

In this study, we developed a rapid, sensitive and efficient sample preparation method to accurate measure MLT in crude plants extracts down to 100 mg F.W. by employing a high resolution Orbitrap MS and using stable isotope-labeled MLT (MLT-d4) as internal standard. Using this approach, we successfully minimize sample purification process into one-step LLE procedure to purify the crude extracts in rice, which greatly simplified the sample preparation procedure and improved the analytical throughput. The established method was successfully applied to investigate the dynamic distributions of MLT responses to Cd stress in shoot and root of rice. The developed method may facilitate to better understanding the physiological functions and regulatory mechanism of MLT in plant.

### AUTHOR CONTRIBUTIONS

TY designed and performed the experiments, analyzed the data, and wrote the manuscript; Y-HH performed the experiments and analyzed the data; LY performed the experiments; HS and RR revised the manuscript; Y-QF designed the experiments and revised the manuscript; and all authors approved the manuscript.

### ACKNOWLEDGMENTS

We thank Prof. Jiaqin Shi (Oilcrops Research Institute, Chinese Academy of Agricultural Sciences, China) for providing rape leaves, and thank Prof. Yuxian Zhu (College of Life Science, Wuhan University, China) for providing cotton leaves and fibers. This research was supported by the National Natural Science Foundation of China (31670373, 21475098).

## REFERENCES

fpls-08-00064 January 23, 2017 Time: 15:30 # 9

in Solanum lycopersicum L. Front. Plant Sci. 6:601. doi: 10.3389/fpls.2015. 00601

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Ye, Hao, Yu, Shi, Reiter and Feng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fpls-08-00064 January 23, 2017 Time: 15:30 # 10

# Comparative Transcriptomic Analyses of Differentially Expressed Genes in Transgenic Melatonin Biosynthesis Ovine HIOMT Gene in Switchgrass

### Edited by:

*Haitao Shi, Hainan University, China*

### Reviewed by:

*Chunyu Bai, Chinese Academy of Agricultural Sciences, China Jidi Xu, Northwest A&F University, China*

#### \*Correspondence:

*Yunwei Zhang<EMAIL_ADDRESS>Fuyu Yang<EMAIL_ADDRESS>† These authors have contributed equally to the work.*

#### Specialty section:

*This article was submitted to Plant Metabolism and Chemodiversity, a section of the journal Frontiers in Plant Science*

> Received: *12 September 2016* Accepted: *12 October 2016* Published: *08 November 2016*

### Citation:

*Yuan S, Guan C, Liu S, Huang Y, Tian D, Cui X, Zhang Y and Yang F (2016) Comparative Transcriptomic Analyses of Differentially Expressed Genes in Transgenic Melatonin Biosynthesis Ovine HIOMT Gene in Switchgrass. Front. Plant Sci. 7:1613. doi: 10.3389/fpls.2016.01613* Shan Yuan1 †, Cong Guan1 †, Sijia Liu1 , Yanhua Huang1, 2 , Danyang Tian1 , Xin Cui 1 , Yunwei Zhang1, 3, 4 \* and Fuyu Yang1, 5 \*

*1 College of Animal Science and Technology, China Agricultural University, Beijing, China, 2 College of Agriculture, China Agricultural University, Beijing, China, 3 Beijing Key Laboratory for Grassland Science, China Agricultural University, Beijing, China, 4 National Energy R&D Center for Biomass, Beijing, China, 5 Beijing Sure Academy of Biosciences, Beijing, China*

Melatonin serves pleiotropic functions in prompting plant growth and resistance to various stresses. The accurate biosynthetic pathway of melatonin remains elusive in plant species, while the *N*-acetyltransferase and *O*-methyltransferase were considered to be the last two key enzymes during its biosynthesis. To investigate the biosynthesis and metabolic pathway of melatonin in plants, the RNA-seq profile of overexpression of the ovine *HIOMT* was analyzed and compared with the previous transcriptome of transgenic *oAANAT* gene in switchgrass, a model plant for cellulosic ethanol production. A total of 946, 405, and 807 differentially expressed unigenes were observed in *AANAT* vs. control, *HIOMT* vs. control, and *AANAT* vs. *HIOMT*, respectively. Two hundred and seventy-five upregulated and 130 downregulated unigenes were detected in transgenic *oHIOMT* line comparing with control, including the significantly upregulated (F-box/kelch-repeat protein, zinc finger BED domain-containing protein-3) genes, which were potentially correlated with enhanced phenotypes of shoot, stem and root growth in transgenic *oHIOMT* switchgrass. Several stress resistant related genes (SPX domain-containing membrane protein, copper transporter 1, late blight resistance protein homolog R1A-6 OS etc.) were specifically and significantly upregulated in transgenic *oHIOMT* only, while metabolism-related genes (phenylalanine-4-hydroxylase, tyrosine decarboxylase 1, protein disulfide-isomerase and galactinol synthase 2 etc.) were significantly upregulated in transgenic *oAANAT* only. These results provide new sights into the biosynthetic and physiological functional networks of melatonin in plants.

Keywords: switchgrass, melatonin, transgene, oHIOMT, RNA-seq

## INTRODUCTION

Melatonin (N-acetyl-5-methoxytryptamine) was first discovered in vertebrates (Lerner et al., 1958), then detected in higher plants (Dubbels et al., 1995; Hattori et al., 1995), and now is widely accepted of its distribution in all kingdoms, from prokaryotes to eukaryotes, from animals to plants (Manchester et al., 2000; Hardeland and Poeggeler, 2003; Simopoulos et al., 2005; Hardeland, 2015), and functions as a direct scavenger of reactive oxygen species (ROS), a mediator hormone of circadian rhythms and an activator of antioxidant enzymes (Arnao and Hernández-Ruiz, 2015; Reiter et al., 2015; Bai et al., 2015, 2016; Gao et al., 2016; Shi et al., 2016a). In plants, melatonin regulated the seed germination, growth of roots, and shoots, and development of flowering (Hernández-Ruiz et al., 2005; Zhao et al., 2013; Arnao and Hernández-Ruiz, 2015). Moreover, the responses to abiotic stresses, including extreme temperature, drought, salt, radiation, make melatonin a concerned candidate as a natural stimulator for crop cultivation (Li et al., 2012; Zhang et al., 2014a; Fan et al., 2015; Shi et al., 2015a). However, the involved regulations of the functional gene expression and physiological mechanism of melatonin biosynthesis and metabolic pathways remains poorly understood in plants (Hardeland, 2015; Zhang et al., 2014b).

The classic pathway of melatonin biosynthesis comprises four steps beginning with tryptophan, firstly decarboxylation by tryptophan decarboxylase (TDC), then N-acetylation by arylalkylamine N-acetyltransferase (AANAT) in animals/serotonin N-acetyltransferase (SNAT) in plants, and the final O-methylation by N-aceylserotonin methyltransferase (ASMT) in animal/hydroxyindole-O-methyltransferase (HIOMT) in plants (Tan et al., 2016). The last two steps are presumed as rate-limiting which are catalyzed by AANAT/SNAT and HIOMT/ASMT (Morton and Forbes, 1988; Byeon et al., 2015, 2016). Furthermore, the overexpression of rice caffeic acid O-methyltransferase (COMT) exhibited upregulation of melatonin contents in transgenic rice plants, indicating that the N-acetylserotonin methyltransferase activity was required for melatonin biosynthesis (Byeon et al., 2015). Recently, genetic engineering modifications of coding genes for melatonin biosynthesis and metabolism enzymes were applied to alter the melatonin contents in transgenic rice, tomato, Arabidopsis thaliana, and Nicotiana sylvestris (Kang et al., 2010; Okazaki et al., 2010; Park et al., 2012; Zhang et al., 2012; Wang et al., 2014; Zuo et al., 2014).

Switchgrass (Panicum virgatum) is a Poaceae warm season C4 perennial grass native to the U.S., and is regarded as a model plant of cellulosic biofuel production for its desirable characteristics, such as large biomass and strong ability to thrive in marginal areas (McLaughlin and Kszos, 2005; Keshwani and Cheng, 2009; Brown et al., 2016). Molecular breeding has made much progress toward improving biomass yield, biofuel quality, and stress resistance in switchgrass, which is the model plant for cellulose ethanol production (Bouton, 2007; Fu et al., 2011; Shen et al., 2013; Poovaiah et al., 2014; Wuddineh et al., 2015). Recently, a series of transgenic overexpression of transcriptional factor (NAC, AP2/ERF, MYB4), sucrose synthesis gene (SUS) to prompt growth (Shen et al., 2013; Poovaiah et al., 2014; Yang et al., 2014; Wuddineh et al., 2015), the inhabitation of cinnamyl alcohol dehydrogenase (CAD) to reduce the lignin content (Fu et al., 2011), and the exploration of abiotic stress related genes and miRNA (Sun et al., 2012; Sharma et al., 2015; Liu et al., 2016) were produced in switchgrass. In addition, the activities of antioxidant and free radical scavenger for melatonin provide opportunities for prompting growth and development in plants (Arnao and Hernández-Ruiz, 2015; Reiter et al., 2015; Shi et al., 2016a). The exogenous applications of melatonin exhibited the enhanced seed germination in cucumber (Zhang et al., 2014a), lateral root formation in both Brassica juncea (Chen et al., 2009) and cucumber (Zhang et al., 2014c), and salt resistance in soybean plants (Wei et al., 2015) etc. Otherstudies revealed that melatonin also exercised some control over root architecture as observed in St. John's Wort, wild leaf mustard, sweet cherry root stocks, and lupin (Murch et al., 2001; Arnao and Hernández-Ruiz, 2007; Chen et al., 2009; Sarropoulou et al., 2012). Moreover, the endogenous modifications of related genes to gain the melatoninrich plants displayed cold resistance in rice (Kang et al., 2010), the drought-tolerant phenotypes of tomato (Wang et al., 2014) and A. thaliana (Zuo et al., 2014). Here, the HIOMT gene encoding the last enzyme in melatonin biosynthesis pathway was overexpressed in switchgrass, and the transcriptomic profile was analyzed in order to disentangle the melatonin biosynthesis pathway and also the potential functions of melatonin in plants.

### MATERIALS AND METHODS

### Plant Materials and Morphological Traits

Transgenic switchgrass (cultivar Alamo) lines expressing ovine HIOMT gene (ID: JF815374.1) were grown under 16 h light (26◦C, 120µmol/m2 /s) and 8 h dark (18◦C) conditions in greenhouse. Fully matured plants were chosen from each genotype for molecular characterization and transcriptome sequencing. Three replicates of stems for each transgenic line and in total six (three transgenic HIOMT (H) lines: H1, H2, H3, and three transgenic empty vector (EV): EV1, EV2, and EV3) were collected and frozen in liquid nitrogen and stored at 80◦C until analysis. The t-test was applied to compare the differences of the morphological traits between transgenic oHIOMT lines and control.

At the transgenic reproductive third (R3) stage (Hardin et al., 2013), the tiller number, plant height, stem node number, internode length, internode diameter, leaf blade length, leaf blade width, root number, root length, root diameter, and spike length were determined (**Figure 1**). The third internode was chosen for measuring internode diameter. Leaf blade length and leaf blade width of third internode were measured. Twelve replicates were randomly sampled for each transgenic line.

### RNA Isolation and Qualification

Total RNA was extracted from the sampling stems using Trizol method (Invitrogen, Carlsbad, CA, USA). RNA purity and integrity were assessed using the NanoPhotometer R spectrophotometer (IMPLEN, CA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA),

respectively. RNA concentration was measured using Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer.

### Transcriptome Preparation

The 1.5µg RNA per sample was prepared for the RNA-seq. Sequencing libraries were generated using NEBNext R UltraTM RNA Library Prep Kit for Illumina R (NEB, USA). The cDNA fragments of 150∼200 bp were selected from the library by purification with AMPure XP system (Beckman Coulter, Beverly, USA).

### Clustering

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia). After cluster generation, paired-end reads were generated by sequencing with the library preparations on an Illumina Hiseq platform.

### Validation of RNA-seq Data by Real-Time Quantitative

RNA-seq results were validated by Real-time quantitative PCR of 16 different genes using 7500 Real-Time PCR System (Applied Biosystems) (primer sequences please see Supplementary Table 1). Gene expression levels were calculated by the 2-11Ct method (Livak and Schmittgen, 2001). Each plate was repeated three times in independent runs for all reference and selected genes.

## Data Analysis

### Quality Control

After removing the adapter containing reads, ploy-N containing reads and low quality reads, clean data/clean reads were obtained from raw data. Meanwhile, the Q20, Q30 values, GC-contents, and sequence duplication level of the clean data were calculated (Supplementary Table 2). The RNA-seq raw data have been submitted to NCBI SRA dataset (http://www.ncbi.nlm.nih.gov/ bioproject), and the project accession number is PRJNA322585.

### Transcriptome Assembly and Gene Functional Annotation

De novo transcriptome assembly was constructed and accomplished for the incomplete genome of switchgrass using Trinity with min\_kmer\_cov set to 2 and all other parameters set default (Grabherr et al., 2011). Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology).

### Differential Expression Analysis

Gene expression levels were calculated by RSEM (Li and Dewey, 2011). To reveal the transcriptionally regulatory event occurring during the transgenic process, comparative transcriptomic analysis was performed among the pools of control and the two transgenic (oAANAT and oHIOMT) RNA samples. Only the genes with a P-value-adjusted (padj) < 0.05 were identified as being significantly changed by the Benjamini and Hochberg's approach.

### GO Enrichment Analysis and KEGG Pathway Enrichment Analysis

Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs. KEGG (Kanehisa et al., 2008) is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem, from molecular-level information, especially largescale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). The KOBAS software was performed to test the statistical enrichment of differential expression genes in KEGG pathways (Mao et al., 2005).

### RESULTS

### Promoted Vegetative and Reproductive Growth

To compare the effects of HIOMT gene during the biosynthesis processes, the phenotypic traits were measured. Transgenic HIOMT switchgrass exhibited significant promotion of growth comparing with EV (**Figure 1**; Supplementary Table 3). Average plant height and internode length were 35.6 and 52.9% higher in HIOMT (90.60 and 13.43 cm) than that of EV (66.81 and 8.78 cm, respectively, P < 0.05). Stem node number of HIOMT was 15.6% higher than that of EV (P > 0.05). There was no significant difference in internode diameter between the two groups (P > 0.05). Average leaf blade length was 20.7% longer in HIOMT (52.53 cm) than that of EV (43.53 cm, P < 0.05). Root number, root length and root diameter of HIOMT were 75.8, 18.3, and 39.4% larger than those of EV, respectively (P < 0.05). Spike length in HIOMT (19.5 cm) was nearly 4-fold of the EV (5.4 cm, P < 0.01).

### Differential Expression Profiling between Transgenic oHIOMT and EV Lines

To analyze the similarities and differences in the transgenic transcriptome, a hierarchical clustering was represented and indicated the significant differences in transcripts of all of the DEGs in the three replicates of the control and two transgenic groups (**Figure 2**). After calculating the expression level of each mapped unigene, a total of 1556 unigenes were detected that had levels of expression that were significantly different among the two transgenic and control libraries. The 946, 405, and 807 differentially expressed unigenes were observed in AANAT vs. EV, HIOMT vs. EV, and AANAT vs. HIOMT, respectively (**Figure 3**). A total of 183 genes were differentially expressed in both AANAT vs. EV and HIOMT vs. EV, but not significantly different in AANAT vs. HIOMT. Three unigenes overlapped with all three groups, including c52804\_g3 (hypothetical protein), c56995\_g3 (O-methyltransferase), c56995\_g4 (ribosomal protein). One hundred twenty nine unigenes were upregulated in both transgenic oAANAT and oHIOMT comparing with EV, while 55 unigenes were downregulated. There were 145 and 76 genes that significantly upregulated and downregulated in transgenic oHIOMT, respectively, but not differentially expressed in transgenic oAANAT. Comparatively, there were 607 and 153 genes that significantly upregulated and downregulated in transgenic oAANAT, while not differentially expressed in transgenic oHIOMT (Supplementary File 1). Twenty-five unigenes specifically and differentially expressed in HIOMT, including 13 upregulated (c64524\_g1, c52718\_g1: glycine-rich domain-containing protein 1-like; c43728\_g1: Bowman-Birk type wound-induced proteinase inhibitor WIP1-like; c62742\_g1: SPX domain-containing membrane protein) and 12 downregulated ones (c46269\_g1: NAC domain-containing protein 67-like; c61012\_g1: protein tyrosine kinase) comparing with AANAT and EV (Supplementary File 1).

### Functional Classification of the DEGs by Gene Ontology Analysis

Among all differentially expressed unigenes, 1071 were upregulated among the three libraries, and 513 were downregulated. However, only 511, 193, and 504 unigenes were functionally annotated with GO terms in AANAT vs. EV, HIOMT vs. EV, and AANAT vs. HIOMT, respectively, revealed by DEG analysis were functionally assigned to the relevant terms in three categories (Biological Process, Cellular Component, and Molecular Function) of the GO database (**Figure 4**). The GO terms "organonitrogen compound biosynthetic process," "oxidoreductase activity," "NAD(P)H oxidase activity," growth factor activity, and "hyperosmotic response" were highly enriched in the DEGs. At the same time, other terms that were related to the response to various other types of abiotic and biotic stresses, such as SOS response "phototropism," "cellular response

to external stimulus," and "response to external stimulus," were also highly enriched in the DEGs.

By comparing transgenic AANAT vs. EV, 2012 GO terms were DEGs, including 1687 upregulated and 1094 downregulated GO items. The GO terms of the "peptide metabolic process" and "amide biosynthetic process" in Biological Process and "structural molecule activity" in Molecular Function were significantly overrepresented (**Figure 4**). By comparing transgenic HIOMT vs. EV, 1314 DEGs, including 724 downregulated GO items and 1062 upregulated GO items, were functionally assigned to the relevant terms; "metabolic process" in Biological Process was significantly overrepresented, followed by "oxidation-reduction process" in Biological Process and "oxidoreductase activity" in Molecular Function. Notably, none of the GO terms was significantly enriched after multiple testing corrections (corrected P > 0.05) in HIOMT vs. EV. Considering the AANAT vs. HIOMT, 1904 DEGs, including 1493 upregulated and 1126 downregulated GO items, were significantly overrepresented of the "macromolecule biosynthetic process" and "cellular macromolecule biosynthetic process" in Biological Process, "cytoplasm," and "cytoplasmic part" in Cell Component, and "structural molecule activity" and "structural constituent of ribosome" in Molecular Function (**Figure 4**).

### KEGG Pathway Analysis of the Melatonin-Related Genes

The 50 pathways were identified between the HIOMT vs. EV from the KEGG database, with 36 upregulated and

17 downregulated pathways (**Figure 5**). The DEGs were significantly enriched in "homologous recombination," "mineral absorption," "isoquinoline alkaloid biosynthesis," "tropane, piperidine, and pyridine alkaloid biosynthesis," "nitrogen metabolism," "flavonoid biosynthesis," "steroid biosynthesis," "alpha-Linolenic acid metabolism," "tyrosine metabolism," and "beta-Alanine metabolism." The top 20 obviously enriched pathways are shown in **Figure 5**. Specifically, the photosynthesisantenna proteins pathway, MAPK signaling pathway and oxidative phosphorylation were significantly upregulated in KEGG pathway maps during the transgenic oAANAT line. While nitrogen metabolism, flavonoid biosynthesis, beta-Alanine metabolism, Brassinosteroid biosynthesis, Phenylalanine metabolism, and Ascorbate and aldarate metabolism were significantly upregulated in transgenic oHIOMT line comparing with EV (Supplementary File 2). Furthermore, a total of 166 pathways, with 141 upregulated and 66 downregulated pathways, represented the DEGs between AANAT and HIOMT groups, including "ribosome," "oxidative phosphorylation," "photosynthesis," "DNA replication" and "calcium signaling pathway."

### Validation of Gene Expression Profiles Using RT-qPCR

The 16 DEGs were selected to evaluate the different expressions for qRT-PCR. Total RNA samples extracted from switchgrass leaves at reproductive stages were used as templates. Histograms were generated by comparing the FPKM determined by transcriptome analysis and qRT-PCR. Expression quantities of the selected genes using RT-qPCR were consistent with the results obtained with RNA-Seq analysis, which means the RNA-seq data were credible. (R 2 = 0.861 for AANAT, R 2 = 0.933 for HIOMT, P < 0.01) (**Figure 6**).

### DISCUSSION

The catalytic reactions of the two key enzymes (AANAT and HIOMT) are the last two steps during the melatonin biosynthesis processes, and the overexpression of the two genes significantly enhanced the melatonin levels in rice (Kang et al., 2010) and tomato (Wang et al., 2014). In our study, the average melatonin contents of transgenic oAANAT and oHIOMT lines was 12 and 36% (P < 0.05) higher than that of EV control, respectively. Melatonin aids plants in terms of root growth, leaf morphology, chlorophyll preservation, and fruit development (Reiter et al., 2015). The enhanced phenotypic growth traits in transgenic oHIOMT switchgrass were in accordant with the increasing of the melatonin contents, and average shoot height, leaf length, root number, root length, and spike length of transgenic HIOMT lines surpassed those of EV control (**Figure 1**; Supplementary Table 3). The transcriptional level of the related genes was significantly upregulated, for instance, F-box/kelch-repeat protein, zinc finger BED domaincontaining protein 3, RAX2-like protein and so on (**Table 1**). Other researches demonstrated that increasing of melatonin significantly stimulate vegetative growth by both exogenous and endogenous applications in several plant species (Hernández-Ruiz et al., 2005; Wang et al., 2014).

Furthermore, the transgenic oHIOMT switchgrass exhibited earlier flowering than that of the EV, suggesting that the reproductive development was expedited by the endogenous addition of melatonin contents (**Figure 1**). The reduction and delay of flowering were reported under the exogenous melatonin treatment of Chenopodium rubrum (Koláˇr et al., 2003) and transgenic melatonin-rich rice (Byeon and Back, 2014), respectively. Melatonin is potentially involved in the regulation of flowering process by the similar pattern of the antioxidant ascorbic acid (Kotchoni et al., 2009) or the impediment of the floral transition from vegetative growth to reproductive growth by the repressor of gibberelic acid (GA) pathway (Shi et al., 2016b). However, other researches implied that melatonin was transiently induced with a peak level during flower development in rice (Park et al., 2013). In our study, one of the flowering regulator encoding gene APETALA2 significantly upregulated in transgenic oHIOMT line (log2FoldChange = 6.11 for c52062\_g5, P < 0.01, **Table 1**) comparing with the EV control, which acts on the regulation of establishment of the floral meristem (Chen, 2004). The accumulation of melatonin represented the resistance for internal and external oxidative stress during reproductive development (Park et al., 2013).

Melatonin also plays innate immune responses to the complex of biotic and abiotic stresses in various plants species, dicots and monocots (Lee et al., 2014; Shi et al., 2015b; Hardeland,

component (CC), and molecular function (MF).

factor. A high *q*-value is represented by blue, and a low *q*-value is represented by red.

2016). As summarized in recent publications, increases of melatonin were typically induced by diverse forms of stresses, including extreme temperatures, drought, salinity, and oxidant (Li et al., 2012; Zhang et al., 2014a; Fan et al., 2015). From our results, the overexpression of oHIOMT gene also drives a series of defense genes activated, for instance, c43076\_g1 gene (annotated to "response to stress") was significantly upregulated comparing with the EV (Supplementary File 1). In addition, the several amino acid and secondary metabolite related pathways were significantly upregulated, including flavonoid biosynthesis, tyrosine metabolism, beta-Alanine metabolism, glycine, serine, and threonine metabolism, brassinosteroid biosynthesis, phenylalanine metabolism, suggesting the activated growth and stress responses in transgenic oHIOMT lines

(Supplementary File 2). Other researches exhibited similar tolerance to drought in transgenic MzASMT A. thaliana (Zuo et al., 2014), and the resistance to herbicide-induced oxidative damages in transgenic melatonin-rich rice (Park et al., 2012). The transcriptional regulation of defense genes alleviate the oxidative damages driven by various stresses (Wei et al., 2016).

Notably, the differentially expressed unigenes in HIOMT vs. EV (405) was less than those of AANAT vs. EV (946), suggesting that the overexpression of oAANAT could affect more growth-related genes than that of oHIOMT. In contrast, the significantly enriched pathways of the DEGs in HIOMT were more than those of AANAT, indicating that the alterations of involved processes driving by the O-methyltransferase were more complex (**Figures 3**, **5**). The previous RNA-seq analysis of overexpression of oAANAT gene reported a large number of differentially expressed unigenes, which were majorly involved in various compound biosynthetic processes and organelle component (Yuan et al., 2016). Here, the synchronously upregulated transcriptional factors and growth related genes in transgenic oAANAT and oHIOMT implied the prompting roles of melatonin in plant growth and development comparing with EV. In addition, a number of differentially expressed genes specifically driving by transgenic oAANAT or oHIOMT potentially indicated the typical functions of each enzyme. For instance, several stress resistant related genes (SPX domaincontaining membrane protein, copper transporter 1, late blight resistance protein homolog R1A-6 OS etc.) were specifically and significantly upregulated in transgenic oHIOMT only, while metabolism-related genes (phenylalanine-4-hydroxylase, tyrosine decarboxylase 1, protein disulfide-isomerase, and galactinol synthase 2 etc.) were significantly upregulated in transgenic oAANAT only (**Table 1**). The previous transgenic tomato represented the branching phenotype in overexpressing of oAANAT and oHIOMT and drought tolerance in transgenic oHIOMT lines (Wang et al., 2014). These indicated the both similarity and differences of the functions of the two key enzymes in the melatonin biosynthesis and metabolic networks. Furthermore, several researches regard the AANAT/SNAT is the rate-limiting enzyme during the melatonin biosynthesis (Morton and Forbes, 1988), while the other recent studies consider that the HIOMT/ASMT is the functional one (Byeon et al., 2015, 2016). The previous observation of the nearly two times of leaf melatonin level in oHIOMT line of transgenic tomato than those in oAANAT lines, provided the proof that ASMT, the homologous to HIOMT, was possibly the ratelimiting enzyme in plants (Wang et al., 2014). Although the similar patterns of the more melatonin accumulation and metabolic pathways were involved in the overexpression of the oHIOMT lines than those from the oAANAT ones in switchgrass seem to support the HIOMT is the rate-limiting enzyme (**Figure 5**), the large amount of KEGG pathways related to animal diseases possibly indicated the limited knowledge of melatonin functions in plants through the animaloriented transgene. Therefore, further investigations from plantoriented transgenic donor will be expected to provide valuable information for rate-limiting enzyme during the melatonin biosynthetic process.

The multiple membrane receptors and signal transduction mechanism were clarified in animals (Reppert et al., 1994). However, no high-affinity melatonin receptor was identified in plants to date (Park, 2011). Therefore, the receptor-independent mechanism of melatonin in plants was proposed as the neurons (Jan et al., 2011). In our study, the c58996\_g3 unigene significantly downregulated (log2FoldChange = −1.768, padj < 0.01), which was annotated as G-protein coupled receptor signaling pathway by GO classification. The discrepancies between increasing melatonin contents and decreasing signaling receptor roughly sustain the hypothesis that melatonin might act via a receptor independent mechanism.

In conclusion, the transcriptomic data from transgenic oHIOMT gene switchgrass revealed amount of differentially expressed unigenes comparing with transgenic oAANAT and

*U, significantly up-regulated (corrected P* < *0.05); D, significantly down-regulated (corrected P* < *0.05); and F, not differentially expressed (corrected P* > *0.05). No. 1–8 and No. 11–15 genes indicated the specific up-regulation and down-regulation in H; No. 9 gene was upregulated in H but downregulated in A line, while No. 10 was oppositely downregulated in H but upregulated in A line; No. 16–34 and No. 35–43 showed the specific up-regulation and down-regulation in A line. NA or Inf presented if the adjusted readcount is zero.*

control lines, implying that the last step of catalytic reaction probably is the rate-limiting step during the melatonin biosynthesis in plants. The potential roles of melatonin in plants were indicated by the upregulation of a series of transcriptional factors and functional genes involving growth and resistance for various stresses. Moreover, additional experiments under adverse environmental stresses will confirm the melatonin physiological functions in plants. The definitions of the membrane receptors and signal transduction will largely drive the explanation of the melatonin biosynthesis processes and functional metabolic pathways.

### AUTHOR CONTRIBUTIONS

Conceived and designed the experiments: YZ and FY. Performed the experiments: SY, CG, SL, and YH. Analyzed the data: SY, CG, YH, XC, and DT. Wrote the paper: SY, CG, and YZ. All authors reviewed and approved the final manuscript.

### REFERENCES

### ACKNOWLEDGMENTS

Funding for this work was provided by the Ministry of Science and Technology, China (2014BAD23B03, 2012AA101801), National Natural Science Foundation of China (31272493), Natural Science Foundation of Beijing (6162016) and China Agricultural University (2014FG062).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01613/full#supplementary-material

bovine amniotic epithelial cells and recovery from spinal cord injury. J. Pineal Res. 32, 147–152. doi: 10.1111/jpi.12311

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2016 Yuan, Guan, Liu, Huang, Tian, Cui, Zhang and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# RNA-seq Analysis of Overexpressing Ovine AANAT Gene of Melatonin Biosynthesis in Switchgrass

Shan Yuan1 , Yanhua Huang1,2, Sijia Liu1 , Cong Guan1 , Xin Cui1 , Danyang Tian1 , Yunwei Zhang1,3,4 \* and Fuyu Yang1,5 \*

1 College of Animal Science and Technology, China Agricultural University, Beijing, China, 2 College of Agriculture, China Agricultural University, Beijing, China, 3 Beijing Key Laboratory for Grassland Science, China Agricultural University, Beijing, China, 4 National Energy R&D Center for Biomass, Beijing, China, 5 Beijing Sure Academy of Biosciences, Beijing, China

Melatonin serves important functions in the promotion of growth and anti-stress regulation by efficient radical scavenging and regulation of antioxidant enzyme activity in various plants. To investigate its regulatory roles and metabolism pathways, the transcriptomic profile of overexpressing the ovine arylalkylamine N-acetyltransferase (oAANAT) gene, encoding the penultimate enzyme in melatonin biosynthesis, was compared with empty vector control using RNA-seq in switchgrass, a model plant of cellulosic ethanol conversion. The 85.22 million high quality reads that were assembled into 135,684 unigenes were generated by Illumina sequencing for transgenic oAANAT switchgrass with an average sequence length of 716 bp. A total of 946 differentially expression genes in transgenic line comparing to control switchgrass, including 737 up-regulated and 209 down-regulated genes, were mainly enriched with two main functional patterns of melatonin identifying by gene ontology analysis: the growth regulator and stress tolerance. Furthermore, KEGG maps indicated that the biosynthetic pathways of secondary metabolite (phenylpropanoids, flavonoids, steroids, stilbenoid, diarylheptanoid, and gingerol) and signaling pathways (MAPK signaling pathway, estrogen signaling pathway) were involved in melatonin metabolism. This study substantially expands the transcriptome information for switchgrass and provides valuable clues for identifying candidate genes involved in melatonin biosynthesis and elucidating the mechanism of melatonin metabolism.

#### Keywords: switchgrass, melatonin, overexpression, oAANAT, RNA-seq

### INTRODUCTION

Since melatonin (N-acetyl-5-methoxytryptamine) was first discovered in plant in Dubbels et al. (1995) and Hattori et al. (1995), varying concentrations of melatonin have been found in many other plant species (Murch et al., 1997; Manchester et al., 2000; Simopoulos et al., 2005). The fundamental issues of melatonin biosynthetic pathways and physiological functions still need to be deciphered to utilize in plants (Tan et al., 2003, 2012; Arnao and Hernández-Ruiz, 2009; Arnao, 2014; Wei et al., 2016). Both exogenous melatonin treatments (Kolár et al., 2003; Arnao and Hernández-Ruiz, 2007; Sarropoulou et al., 2012; Zhang et al., 2013) and melatonin-rich transgenic plants (Kang et al., 2010; Park et al., 2012; Byeon and Back, 2014; Wang et al., 2014) have been carried on in order to determine its potential functional roles. The physiological role of

### Edited by:

Haitao Shi, Hainan University, China

#### Reviewed by:

Ruediger Hardeland, University of Göttingen, Germany Yang-Dong Guo, China Agricultural University, China

#### \*Correspondence:

Yunwei Zhang<EMAIL_ADDRESS>Fuyu Yang<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Plant Genetics and Genomics, a section of the journal Frontiers in Plant Science

Received: 15 July 2016 Accepted: 12 August 2016 Published: 31 August 2016

#### Citation:

Yuan S, Huang Y, Liu S, Guan C, Cui X, Tian D, Zhang Y and Yang F (2016) RNA-seq Analysis of Overexpressing Ovine AANAT Gene of Melatonin Biosynthesis in Switchgrass. Front. Plant Sci. 7:1289. doi: 10.3389/fpls.2016.01289

melatonin in plants involved growth regulation, scavenging reactive oxygen species and increases of antioxidant enzyme activities (Paredes et al., 2008; Arnao and Hernández-Ruiz, 2014; Hardeland, 2016). Through the encoding genes for the catalytic reactions during melatonin biosynthesis have been identified and cloned recently, its biosynthetic pathways still not clearly deciphered yet in plants (Kang et al., 2011, 2013; Byeon et al., 2013a; Hardeland, 2016). The classic pathways of melatonin biosynthesis consist of the four catalytic reactions from tryptophan: decarboxylation by tryptophan decarboxylase (TDC), hydroxylation by tryptamine 5-hydroxylase (T5H) to serotonin, N-acetylation by serotonin N-acetyltransferase (SNAT) and the final O-methylation to melatonin by N-acetylserotonin O-methyltransferase (ASMT). The penultimate step catalyzes the same reaction as the nonhomologous AANAT of vertebrate, which catalyzes conversion of serotonin into N-acetylserotonin (Byeon et al., 2013a). Overexpression of AANAT gene can promote the ability of biosynthesis of melatonin, and significantly improve the melatonin content in plants (Kang et al., 2010; Wang et al., 2014; Zhang et al., 2014).

Switchgrass (Panicum virgatum L.) is a perennial C4 grass native to North America, and is well-researched on germplasm collection, cultivation, genetic breeding as a model plant for cellulosic bioethanol production over the past several decades (Sanderson et al., 1996; Parrish and Fike, 2005). Aiming at the production demand of large biomass, the traditional breeding techniques have been challenged by its nature of outcrossing and polyploidy, as well as the general infertile environments for the biomass grass cultivation. Fortunately, modification of the functional genes that related to the growth and resistance could effectively assist the breeding process of switchgrass (Fu et al., 2011; Poovaiah et al., 2014; Baxter et al., 2015; Wuddineh et al., 2015). Several transgenic plants expressing the key genes for melatonin biosynthesis enzymes have been demonstrated to increase the resistance to environmental stresses (Kang et al., 2010; Park et al., 2012).

In this study, the transgenic melatonin-rich switchgrass (overexpressing of the oAANAT gene) and the transgenic control with empty vector (EV) were used to conduct a RNA-seq and analyze the effects of melatonin on gene expression. This research intends to provide differentially expressed genes (DEGs) in the melatonin-rich switchgrass and mechanism information on the biosynthesis, regulation, metabolism of melatonin for further investigations in plants.

### MATERIALS AND METHODS

### Plant Materials and Cultivation

The switchgrass (Panicum virgatum L. var. Alamo) plants expressing sheep AANAT gene (ID: 25120) were used for this study. The transgenic plants were grown under 16 h light (26◦C, 120 µmol/m2 /s) and 8 h dark (18◦C) conditions with watered every other day. Fully matured plants were chosen from each genotype for molecular characterization and transcriptome sequencing, which were all phenotypically identical with respect to overall size, tiller number and leaf pattern. Three replicates of relative high melatonin contents for each transgenic stem were frozen in liquid nitrogen and stored at 80◦C until further analysis.

### Characterization of Growth and Development of Transgenic Switchgrass

Tiller number, plant height, internode number, internode length, internode diameter, leaf blade length, leaf blade width, root number, root length, root diameter, and spike length were measured at the transgenic reproductive third (R3) stage (**Figure 1**), a stage with fully emerged spikelets and an emerged peduncle (Hardin et al., 2013). Internode 3 (I3) was used for measuring internode diameter. The leaves of I3 were used to measure leaf blade length and leaf blade width. Twelve replicates were measured for each transgenic line.

### RNA Isolation and Qualification

Total RNA was extracted from the stems using the TRIzol reagent method (Invitrogen, Carlsbad, CA, USA) and was treated with RNase-free DNase or 30 min at 37◦C to remove residual DNA. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer spectrophotometer (Implen, Westlake Village, CA, USA). RNA concentration was determined by Qubit RNA Assay Kit in Qubit 2.0 Flurometer. RNA integrity was assessed with the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

### Transcriptome Sample Preparation for Sequencing

The total amount of 1.5 µg RNA per sample was prepared for the RNA-seq. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA). The cDNA fragments of 150∼200 bp were preferentially selected from the library by purification with AMPure XP system (Beckman Coulter, Beverly, MA, USA).

### Clustering and Sequencing

The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia). After cluster generation, paired-end reads were generated by sequencing with the library preparations on an Illumina Hiseq platform.

### Validation of RNA-seq Data by Real-Time Quantitative

Real-time quantitative PCR validation of RNA-seq data for 12 random genes was performed in 20 µl of reaction mixture (Supplementary Table S1). The amount of the amplified DNA was monitored by fluorescence at the end of each cycle using 7500 Real-Time PCR System (Applied Biosystems). Each plate was repeated three times in independent runs for all reference and selected genes. Gene expression was evaluated by the 2−11Ct method (Livak and Schmittgen, 2001).

### Data Analysis

### Quality Control

Clean data (clean reads) were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. Q20, Q30, GC-contents and sequence duplication level of the clean data were calculated. Our raw data have been uploaded to NCBI1 , and the accession number is PRJNA322585.

### Transcriptome Assembly

The complete genome of switchgrass has not been released through the draft version was released in 2012 and updated to V3. De novo transcriptome assembly was accomplished using Trinity (Grabherr et al., 2011) with min\_kmer\_cov set to 2 by default and all other parameters set default.

### Gene Functional Annotation

Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology).

### Differential Expression Analysis

Gene expression levels were estimated by RSEM (Li and Dewey, 2011). Differential expression analysis of two conditions/groups was performed using the DESeq R package (1.10.1). Genes with an adjusted P-value < 0.05 found by DESeq were assigned as differentially expressed using the Benjamini and Hochberg's approach.

### GO Enrichment Analysis and KEGG Pathway Enrichment Analysis

GOseq R packages were implemented to define the DEGs based Wallenius non-central hyper-geometric distribution (Young et al., 2010). We used KOBAS (Mao et al., 2005) software to test the statistical enrichment of differential expression genes in KEGG pathways2 (Kanehisa et al., 2008).

## RESULTS

### Morphological Characterization of Transgenic Plants