Transcriptome analysis of hormone-induced gene expression in Brachypodium distachyon

Brachypodium distachyon is a new model plant closely related to wheat and other cereals. In this study, we performed a comprehensive analysis of hormone-regulated genes in Brachypodium distachyon using RNA sequencing technology. Brachypodium distachyon seedlings were treated with eight phytohormones (auxin, cytokinine, brassinosteroid, gibberelline, abscisic acid, ethylene, jasmonate and salicylic acid) and two inhibitors, Brz220 (brassinosteroid biosynthesis inhibitor) and prohexadione (gibberelline biosynthesis inhibitor). The expressions of 1807 genes were regulated in a phytohormone-dependent manner. We compared the data with the phytohormone responses that have reported in rice. Transcriptional responses to hormones are conserved between Bracypodium and rice. Transcriptional regulation by brassinosteroid, gibberellin and ethylene was relatively weaker than those by other hormones. This is consistent with the data obtained from comprehensive analysis of hormone responses reported in Arabidopsis. Brachypodium and Arabidopsis also shared some common transcriptional responses to phytohormones. Alternatively, unique transcriptional responses to phytohormones were observed in Brachypodium. For example, the expressions of ACC synthase genes were up-regulated by auxin treatment in rice and Arabidopsis, but no orthologous ACC synthase gene was up-regulated in Brachypodium. Our results provide information useful to understand the diversity and similarity of hormone-regulated transcriptional responses between eudicots and monocots.

Human diets greatly depend on grasses, which have applications as energy crops and sustainable energy sources. Rice (Oryza sativa) and maize (Zea mays) are commonly used as model plants for monocots, as their genome sequence and extensive genetic resources are available 1 . Draper et al. proposed Brachypodium (Brachypodium distachyon) as a novel model plant in 2001 because it has advantages such as a small genome size, short life cycle, small plant size, and the ability to self-pollinate 2 . Furthermore, Brachypodium is more closely related to wheat [Triticum aestivum] and rye [Secale cereale] than rice and maize are, which makes it convenient for analysis of different grass groups using model systems 1 . The whole genome sequence of Brachypodium inbred line Bd21 was released in 2010 as the first genome sequence available for the Pooideae subfamily 3 . Recently, a comprehensive collection of full-length cDNAs was reported by Mochida et al. 4 . Therefore, fundamental resources for Brachypodium molecular biology are being rapidly developed. Recently, Priest et al. reported an analysis of global gene expression in Brachypodium in response to abiotic stress 5 . Large-scale gene expression data are essential resources for model systems. For . In this study, we found that the transcriptional regulation of BR-regulated genes in BR-treated det2 plants was negatively correlated with BR-biosynthesis inhibitor treatment using the AtCAST database (http://atpbsmd.yokohama-cu.ac.jp/atcast/html/genelist/A0036_A0081.html) 9 . These findings suggest that the application of hormone biosynthesis inhibitors (Brz220 for BR and Phx for GA) to generate hormone-deficient conditions is a practical means of identifying BR-or GA-regulated genes in Brachypodium.

Results
Determination of hormone treatment conditions and identification of hormone-regulated genes. We selected phytohormone-responsive marker genes in Brachypodium using the following procedure. Twenty-eight phytohormone-responsive genes from rice 13 and Arabidopsis thaliana 11,15 were selected and their sequences were subjected to a BLAST search for the most similar sequences in Brachypodium. Top-hit Brachypodium genes were used to examine transcriptional responses in Brachypodium (Table 1). We checked the transcriptional responses with RT-PCR. In total, 19 of 28 genes showed reproducible gene expression responses to one of the phytohormones or the inhibitors in more than two independent experiments. These 19 genes were selected as phytohormone-responsive marker genes in Brachypodium (Table 1). Next, hormone treatment conditions in Brachypodium were determined using these marker genes. We used 1, 3, and 10 μ M for IAA, tZ and GA; 10, 30, and 100 μ M for Phx, ABA, MJ, ACC, Brz220 and SA; and 1, 10, 100, and 1000 μ M for BL. The concentrations of phytohormones and inhibitors used in the comprehensive transcriptome analysis were determined accordingly. The optimum concentration of each phytohormone (Table 2) was determined as that at which the marker gene showed the greatest change in response to treatment compared to mock treatment. After the determination of treatment conditions, the hormone treatment and sampling were repeated twice independently and the hormonal responses were confirmed against the transcriptional responses of the marker genes. Total RNAs extracted in two independent experiments were mixed and used in RNA-seq analysis.
About 90% of obtained reads from the RNA-seq analysis were successfully mapped to the Brachypodium genome (Table 3). Differentially expressed genes (DEGs) in response to each treatment (FDR < 0.05) were analyzed using edgeR 16 in a comparison with all other treatments ( Table 4). Hundreds of genes were differentially expressed in response to ABA and MJ. Dozens of DEGs were detected for the IAA, tZ and SA treatments. We designated these analyses as being high stringency; the hormone-responsive genes are listed in Data S1-S8. In contrast, only four genes were differentially expressed in the Brz220 treatment. No DEG was detected in the mock 1-h, mock 3-h, Phx, Phx+ GA, Brz220+ BL and ethylene treatments, indicating that the gene expression responses to GA, BR and ethylene were weaker than those to the other hormones. Therefore, we conducted additional statistical analyses to identify the genes responsive to these hormones. For these analyses, DEGs were compared between the hormone treatment vs. two Scientific RepoRts | 5:14476 | DOi: 10.1038/srep14476 control treatments (Table 5). Two hundred and twenty-four DEGs were detected in response to GA, 474 DEGs in response to BR, and 52 DEGs in response to ethylene. We designated these analyses as low stringency; the DEGs are listed in Table 4, Data S9-S16. Comparison with rice and Arabidopsis. First, we compared the phytohormone-regulated genes in Brachypodium with those in rice. A comparison of our experimental conditions with those used in previous rice studies is presented in Table S1. Garg et al. 13 reported the responses of rice seedlings to six hormones: auxin, CK, SA, ethylene, ABA, and JA. We made a list of predicted orthologs of the Brachypodium DEGs in rice and evaluated whether they are regulated in the same manner (Table 6). Since Garg et al. 13 reported phytohormone-regulated DEGs in rice without statistical information for each hormone treatment, we simply checked whether the orthologous genes from Brachypodium were included in the list of Garg et al. 13 or not. We predicted 408 orthologs to the Brachypodium DEGs in rice. In total, 326 orthologous genes (80%) were regulated by phytohormones in the same manner. The transcriptional response to CK and ABA in Brachypodium was consistent with that in rice: 100% of the orthologs of the Brachypodium DEGs that were regulated by CK (tZ) were also regulated in rice by CK (BAP), and 98% of the orthologous genes were regulated by ABA in rice. Brachypodium and rice shared many auxin-and JA-responsive orthologous genes: 67% of the orthologs were also regulated by IAA in rice. Garg et al. 13     Bradi2g46420, Bradi2g31820, Bradi3g55410, Bradi4g35960, Bradi3g54610, and Bradi2g11120) and GH3 family (Bradi2g50840, Bradi1g22830, and Bradi2g52000) were found to be auxin-responsive in Brachypodium. In total, 42% of the orthologous genes were JA-responsive in rice. JAZ1-like genes in Brachypodium (Bradi5g24410), OsJAZ1, a MYC2-like gene in Brachypodium (Bradi3g34200), and OsMYC2 were regulated in both plants. Transcriptional responses to SA and ethylene were less conserved. Only 26% and 31% of the orthologous genes were regulated by SA and ethylene in rice. A transcription factor gene, OsWRKY45, which acts in SA-mediated pathogen resistance 17 , and its orthologous gene in Brachypodium (Bradi2g44270), were found to be regulated by SA in both plants. Ethylene receptor-like genes (Bradi5g00700 and LOC_Os04g08740) were found to be regulated in both plants. Sato et al. 14 performed a transcriptional analysis using auxin, CK, GA, BR, ABA, and JA and presented a coexpression database. Since they did not show phytohormone-regulated genes, we did not directly compare our hormone-induced gene expression results with those of Sato et al. 14 .
A gene ontology (GO) term enrichment (GOE) analysis was performed to evaluate the DEGs for Brachypodium. In this study, the GO terms for Brachypodium were redefined using the latest GO annotation for Arabidopsis thaliana (see the Methods for details) because the GO terms inferred from InterProScan (http://geneontology.org/page/go-annotation-standard-operating-procedures) lack useful GO terms for the study of phytohormone responses (e.g., the "Response to cytokinin stimulus"; Table  S2). Therefore in the following sections, a detailed analysis was performed with the help of GO terms predicted from Arabidopsis (Tables 7 and 8). In this study, we analyzed the GOE for Arabidopsis using the latest GO terms (R package ath1121501cdf version 2.14.0) and microarray data provided by Goda et al. 11 (Table S3) to compare the results of GOE analyses with those in Brachypodium.
Auxin response. One hundred and twenty-five genes were identified as auxin-responsive genes in the high-stringency analysis (Data S1). The genes with the term "response to auxin stimulus" were found to be enriched in the GOE analysis of auxin-responsive genes in Brachypodium ( Table 7). As we described in the comparison with rice, many genes in the Aux/IAA family and GH3 family were detected as auxin-responsive genes in Brachypodium. These results are also consistent with the auxin response in Arabidopsis ( Fig. 1, Fig. S1). In contrast to the auxin response in Arabidopsis, Brachypodium genes related to ethylene functions were not regulated in the same manner ( Fig. S2), although their expression levels were not low (present at sufficiently high levels for analysis) (log 2 cpm (number of counts per million reads) > 1). CK response. Twenty-three genes were identified as cytokinin-responsive genes in the high-stringency analysis (Data S2). The genes with the term "response to cytokinin stimulus" were enriched in the GOE analysis (Table 7). Genes in the type A ARR family (Bradi2g61000, Bradi4g43090, Bradi5g25860, Bradi3g45930 and Bradi5g11350) were detected as cytokinin-responsive genes (Fig. 2). In contrast, members of the type B ARR gene family were not detected as cytokinin-responsive genes. These results are consistent with the cytokinin response in Arabidopsis 18 . SA response. Eighty-one genes were identified as SA-responsive genes in the high-stringency analysis (Data S3). Genes with the term "defense response to fungus" were enriched in GOE analysis (Table 7). A homolog of the PR-1 gene (Bradi1g57590) was up-regulated by SA treatment (Fig. 3).
ABA response. Four hundred and forty-five genes were identified as ABA-responsive genes in the high-stringency analysis (Data S4). Genes with the term "response to abscisic acid stimulus" were enriched in the GOE analysis (Table 7). Genes homologous to those encoding transcription factors involved in transcriptional regulation in response to ABA (i.e., ABF, NAC, Myb, bHLH) were detected as DEGs (  JA response. Three hundred and eighty-three genes were identified as JA-responsive genes in the high-stringency analysis (Data S5). Genes with the term "response to wounding" were enriched in the GOE analysis (Table 7). Genes involved in the biosynthesis of jasmonic acid (LOX3 homolog (Bradi5g11590), (Bradi3g37650) and OPCL1 homolog (Bradi1g76280) were up-regulated by MJ treatment (Fig. 4). The JAZ genes in Arabidopsis encode JA receptors 19 . Two members of the JAZ family (Bradi3g23190, Bradi1g58490) were also up-regulated in Brachypodium.

GA response. No gene was detected as a DEG in the high-stringency analysis (Data S6). This indi-
cates that the GA-inducible gene expression response is relatively weaker than those to other hormones, which is consistent with the GA response in Arabidopsis 11 . When the effects on gene expression of the GA biosysthesis inhibitor Phx treatment were compared to those of the mock 3-h and Phx+ GA treatments (low-stringency analysis), 224 genes were detected as DEGs (Data S14). Of these genes, those with the term "response to nitrate" were the most enriched in the GOE analysis (Table 8). In Arabidopsis, many GA-biosynthesis genes are also GA-responsive, since there is a negative feedback regulation in the expression of GA-biosynthesis genes 11 . The read counts of many genes involved in gibberellin biosynthesis (homologs of GA2oxs, GA3oxs

BR response.
Only four genes were detected as DEGs in the Brz220 treatment compared with all other treatments (high-stringency analysis, Data S7). This indicates that the BR-inducible gene expression response is relatively weak compared with those to other hormones, which is consistent with the BR response in Arabidopsis 11 . When the effects on gene expression of the Brz220 treatment were compared with those of the mock 3-h and Brz220+ BL treatments (low-stringency analysis), 474 genes were detected as DEGs (Data S15). Of these genes, those with the term "RNA elongation" were enriched in the GOE analysis (Table 8). Genes with the term "photosynthesis" were also enriched in the GOE analysis. P450 genes involved in the biosynthesis of brassinosteroids (Bradi1g15030, homolog of BR6ox2 and Bradi5g12990, homolog of DWF4) were up-regulated by Brz220 compared to the mock 3-h and Brz220+ BL treatments (Fig. 5), although the fold change of Bradi5g12990 was less than twice. This is consistent with the BR response in Arabidopsis; BR6ox2 and DWF4 were up-regulated in the det2 mutant of Arabidopsis compared to the wild-type (Fig. 5).
Ethylene response. No gene was detected as a DEG in response to ethylene in the high-stringency analysis (Data S8). To confirm response of seedlings to ACC treatment under our experimental conditions, plants were treated with ACC for a further week. Seedlings showed growth inhibition when grown in the presence of 10 μ M ACC (Fig. S3), indicating that plants respond to exogenous ACC. Therefore, gene expression in response to ACC treatment was compared to that in response to the mock 1-h and mock 3-h treatments (low-stringency analysis). Fifty-two genes were detected as DEGs (Data S16). Of these DEGs, those with the term "RNA elongation" were enriched in GOE analysis (Table 8). A gene homologous to EIN4, which is involved in sensing ethylene (Bradi5g00700) was up-regulated, and a gene homologous to ACS6, which is a member of the ethylene synthesis gene family, (Bradi5g19100) was down-regulated by ethylene treatment compared to the mock 1-h and mock 3-h treatments.

Discussion
The comprehensive transcriptional response of Brachypodium to phytohormones was analyzed in this study. To obtain clear responses from Brachypodium, we optimized the experimental design using the expressions of selected Brachypodium marker genes to validate hormone response. The concentrations of phytohormones used in this study were similar (1-10-fold) to those used in the analysis of the response to hormones in rice and Arabidopsis (Table S1), with the exception of the concentration of BL (100-fold compared to that used for Arabidopsis). We successfully identified BR-and GA-regulated genes (e.g., Bradi1g15030, homolog of BR6ox2 as a BR-responsive gene and Bradi3g18690, homolog of TCH4 as a GA-responsive gene).  Gene IDs log 2 ratio stat cpm In rice, 80% of the orthologs of the Brachypodium DEGs were also regulated by CK, ABA, auxin, JA, SA, and ethylene. These data show that major transcriptional responses to these phytohormones are shared in monocots, although fewer orthologous genes were commonly regulated by SA and ethylene. It was reported that the SA level in rice was the highest among the tested plants 20 . This may be one reason for the different responses to SA at similar concentrations. The difference in ethylene-regulated genes between rice and Brachypodium was also reported by Pacheco-Villalobos et al. 21 . Sato et al. 14 also reported the responses of rice seedlings to six hormones: auxin, CK, GA, BR, ABA, and JA. However, it was difficult to compare our data with theirs because they reported their data as part of a large transcriptome dataset and no analysis of hormone-responsive gene expression was performed using the data.
The transcriptional responses to phytohormones and the mechanisms involved in transcriptional regulation in Arabidopsis are well characterized. In addition, manually curated GO terms for Arabidopsis is very helpful to analyze large transcriptome dataset. Therefore, we compared the transcriptional response to phytohormones in Brachypodium to that in Arabidopsis thaliana in detail. A comprehensive analysis of phytohormone-regulated genes of Brachypodium revealed common transcriptional responses between Brachypodium and Arabidopsis. To compare the accumulation of GO terms in Brachypodium to Arabidopsis, we assigned new GO terms to the Brachypodium genes using the Arabidopsis GO terms. Some results of the GOE analyses in Arabidopsis and Brachypodium were common following treatment with the same phytohormones. The GO terms "response to auxin stimulus" had the lowest p-values in auxin-regulated DEGs. Similarly, "cellular response to cytokinin stimulus" in CK-regulated DEGs, "defense response to fungus" in SA-regulated DEGs, "response to abscisic acid stimulus" in ABA-regulated DEGs and "response to wounding" in JA-regulated DEGs had the lowest p-values in Brachypodium (Tables 7,  8). These GO terms were also enriched in the hormone-responsive genes in Arabidopsis (Table S3). The GO term "response to nitrate" was enriched in DEGs in the Phx treatment (Table 7), which was consistent with the enriched GO of the GA-biosynthesis mutant ga1-5 when compared to the wild type (Table S5). The GO term "photosynthesis" was enriched in DEGs of Brz-treated Brachypodium (Table 7), which was consistent with the enriched GO terms of the BR-biosynthesis mutant det2 when compared to wild-type Arabidopsis. Consistent with the findings in BR-treated Arabidopsis det2, genes involved in BR-synthesis were detected as DEGs in Brz-treated Brachypodium (Fig. 5). Fewer were detected as DEGs in response to ACC than to other hormones (Tables 4, 5). This was also consistent with transcriptional regulation in ethylene-treated Arabidopsis 11,15 . These results suggest that many molecular mechanisms characteristic of the responses to phytohormones in Arabidopsis are also present in Brachypodium.
In contrast, there were some differences in the hormone responses of Brachypodium and Arabidopsis. In Arabidopsis, genes involved in ethylene synthesis and signaling, ACO1, ACO3, ACO4/EFE, ACS4, ACS8,

Figure 2. Phylogenic tree of ARR genes and their transcriptional responses to CK in Brachypodium and
Arabidopsis. Brachypodium ARR family members were retrieved using the BLAST software. Similarities in protein sequences (percentage identity) are shown as a phylogenetic tree. Brachypodium genes are shown in red and Arabidopsis genes in black. log 2 ratio represents the gene expression ratio between the read count in tZ treatments divided by the average read count from all other treatments. Stat represents the p-value of the microarray experiment in Arabidopsis treated with t-zeatin for 1 h 11 or FDR of RNA-seq data when the count from tZ treatment was compared with all other treatments. log 2 ratio is hatched in red if genes are upregulated. cpm represents the average of log 2 -scaled read counts per million reads from all experiments. and ERF1, are responsive to auxin treatment. Only one gene from the ACO gene family (Bradi2g41840) was detected as a DEG in auxin-treated Brachypodium (Fig. S2). Genes involved in ethylene signaling, such as ERF2, 8,11, ERS1, EBF2, ETR2, CTR1 and EIN2, are regulated by ACC treatment in Arabidopsis 15 . However, in Brachypodium, only one gene homologous to EIN4 (Bradi5g00700) was detected as a DEG. These differences in auxin-and ethylene-responsive genes suggest divergent auxin-ethylene relationships between eudicots and monocots. Pacheco-Villalobos et al. 21 reported that the cross-talk between ethylene and auxin in Brachypodium differs from that in Arabidopsis. Further analysis of these differentially expressed genes will improve our understanding of the different functions of these hormones and the related genes in Brachypodium and Arabidopsis.
We encountered some difficulties in the analysis of RNA-seq data. We prepared 13-30 million reads because Liu et al. 22 reported that the increase in number of DEGs drops after 10-15 million reads. However, the transcriptional responses of some genes might not have been detected because of their low read counts. For example, none of the genes involved in GA-biosynthesis was detected as Phx-or GA-responsive in Brachypodium. However, it was reported that the expressions of the GA-biosynthetic genes AtGA2ox1, AtGA2ox2, GA3ox1 and GA20ox1 changed in Arabidopsis GA-deficient mutants following endogenous GA treatment [23][24][25][26] . The read counts of these GA-biosynthetic genes in Brachypodium were low (Table S4). In Arabidopsis and rice, few genes respond to ethylene 11,13 . In this study, no gene in the expansin family was detected as ethylene-responsive in Brachypodium, while AtEXP8 and AtEXP18 in Arabidopsis 27 and OsEXP1, OsEXP2 in rice 13 are responsive. The expression levels of 27 expansin genes were low (cpm < 1, Table S6). The growth of Brachypodium was reduced by ACC treatment (Fig. S3); therefore, some expansin genes may be regulated by ethylene to control cell elongation. To examine the transcriptional regulation of these poorly expressed genes, further analysis using additional RNA-seq reads with different tissues or different growth stages, and/or complementary methods (such as microarray or quantitative PCR) are required. Gene IDs log 2 ratio stat cpm  Table 2). The phytohormones were dissolved in DMSO before addition to the liquid cultures. DMSO (0.1% (v/v)) was used for mock treatments. For the "Inhibitor and GA" and "Inhibitor and BR" treatments, GA and BL were added to liquid cultures 2 h after treatment with the GA inhibitor, prohexadione-calcium (Wako), or BR-inhibitor, brassinazole220 29 , respectively. All growth and treatment processes were performed at 22 °C under continuous light. Phytohormone treatments were conducted twice independently as biological replicates, and the hormone responses were confirmed using the transcriptional responses of marker genes (Table 1). Transcriptional responses of marker genes were checked with RT-PCR using gene specific primers (Table S7) prior to the RNA-seq analysis.
RNA extraction. RNA was extracted from each sample using the RNeasy Mini Kit (QIAGEN). The concentration and quality of the RNA samples were determined using a Nanodrop 2000 instrument (Thermo Scientific). RNA extracted from two biological replicates was mixed and used for the RNA-seq analysis.
Strand-specific RNA-seq analysis. Strand 16 , as follows. Parameters and function usages in edgeR were determined to maximize the specificity for detecting hormone marker genes (Table 1) as DEGs. Read counts were normalized with the function "calc-NormFactors". The TMM normalization method in edgeR can normalize different amounts of RNA-Seq data without increasing the false-positive rate of detecting DEGs 33 . Dispersions among genes and among replications were estimated by the functions "estimateGLMCommonDisp", "estimateGLMTrendedDisp" and "estimateGLMTagwiseDisp" and p-values were calculated using the likelihood-ratio test with the generalized linear models. Calculated p-values were adjusted using the false discovery rate (FDR) of Benjamini and Hochberg's approach 34 . DEGs were defined as FDR < 0.05 and up-or down-regulated genes were identified in the DEGs as those having a log 2 ratio of cpm > 1 or that of cpm < − 1, respectively.
Phytohormone-regulated genes in Arabidopsis. Microarray data from Arabidopsis thaliana treated with phytohormones were retrieved from the AtGenExpress project 11 . Signals in the microarray data were calculated using R package affy (1.42.3) with the MAS5 method 35 . DEGs were identified as p-value < 0.1 by Student's t-test and up-and down-regulated genes were defined in the DEGs as log 2 ratios of signals > 1 and < − 1, respectively. GOE analyses were performed with the DEGs. Gene log 2 ratio stat cpm Figure 5. P450 genes involved in BR biosynthesis and transcriptional responses to Brz in Brachypodium and Arabidopsis. Brachypodium BR biosynthetic genes were retrieved using the BLAST software. Similarities in protein sequences (percentage identity) are shown as a phylogenetic tree. Brachypodium genes are shown in red and Arabidopsis genes in black. log 2 ratio represents the gene expression ratio between read counts from the Brz treatment divided by the average read counts from mock 1-h and Brz+ BL treatments. Stat represents the p-value of the microarray experiment, Arabidopsis det2 compared with wild-type 11 for Arabidopsis genes or FDR of RNA-seq data when the counts from the Brz treatment were compared with mock 1-h and Brz+ BL treatments for Brachypodium genes. log 2 ratio is hatched in red if genes are upregulated and in blue if down-regulated. cpm represents log 2 scaled read counts per million reads of RNAseq data from mock 1-h, Brz and Brz+ BL treatments.