Metabolic and transcriptional regulatory mechanisms underlying the anoxic adaptation of rice coleoptile

Rice is a unique crop plant since it can survive under anoxia condition by flooding and also germinate and grow up to coleoptile through its distinctive adaptations. Despite several decades of research on this topic, the current knowledge on the molecular machinery of rice under anoxia is very limited. Therefore, we unraveled the possible regulatory mechanisms by resorting to systems biology approach which combines the metabolic modeling and transcriptome analysis. Such integrative analysis highlight the critical role of MYB, bZIP, ERF and ZnF transcription factors in up-regulating the fermentation and sucrose metabolism genes to generate sufficient energy for cellular growth.


Introduction
Rice is one of the few plant species that can survive submerged conditions. In particular, it can successfully germinate by extending a coleoptile underwater in the complete absence of oxygen (Magneschi and Perata 2009). Most of the early biochemical studies reported that even under anoxia, rice is able to metabolize reserve carbohydrates via glycolysis and fermentation to † These authors contributed equally. generate sufficient cellular energy, i.e. ATP needed for plant growth (Mocquot et al. 1981;Alpi and Beevers 1983;Mohanty et al. 1993;Guglielminetti et al. 1995). Microarray data show that this is associated with a strong up-regulation of sucrose and starch metabolism, glycolysis and fermentation genes (Lasanthi-Kudahettige et al. 2007). This work also showed many other significantly up-and down-regulated genes associated with various metabolic and signalling pathways. These findings highlight that rice globally reconfigures its molecular machinery under anoxic stress by selectively synthesizing the necessary transcripts coding for enzymes that produce and conserve energy.
Despite past efforts to characterize the biochemical adaptations of rice under anoxia, current knowledge on how oxygen deficiency is sensed and on the regulatory cascade that fine-tunes the transcriptional and translational changes is still very limited. To date, only the induction of GA-response-free RAmy3D under anoxic conditions (Loreti et al. 2003a, b) is the notable trait unravelled at the molecular level. Even with the availability of abundant high-throughput data, such limitations still exist mainly due to the lack of appropriate systematic frameworks to analyse and derive a valid hypothesis from them. In this regard, constraint-based in silico metabolic modelling and analysis can be helpful. Not only do they allow predictions concerning the physiological behaviour and metabolic states of an organism exposed to various environmental or genetic changes but they also serve as a scaffold to contextualize multiple '-omics' data, thereby enabling biologically meaningful correlations to be identified (Lewis et al. 2012). As a result, several frameworks are now available to integrate metabolic models and high-throughput data such as transcriptomics (Becker and Palsson 2008;Colijn et al. 2009;Bordel et al. 2010), metabolomics (Mo et al. 2009;Selvarasu et al. 2012) and proteomics (Shlomi et al. 2008) within the context of systems biology. One such approach compares the flux levels simulated by the metabolic model with corresponding gene expression data resulting from environmental or genetic changes. This has, for example, highlighted key transcriptional mechanisms in Saccharomyces cerevisiae (Famili et al. 2003;Bordel et al. 2010;Osterlund et al. 2013). Therefore, following a similar method, here, we identified the transcriptionally regulated reactions in rice during its anoxic germination. Our approach used initial random flux sampling using the recently published rice central metabolic network (Lakshmanan et al. 2013), and subsequently compared the differences in flux levels with previously published gene expression data (Lasanthi-Kudahettige et al. 2007) between air and anoxia. Furthermore, to gain a deeper insight into the key regulatory mechanisms during anoxic adaptation, we also found the possible transcription factors (TFs) of transcriptionally regulated enzymes by analysing the distribution of putative cis-elements in their promoter regions.

Rice central metabolic network
We used the recently published rice central metabolic model (Lakshmanan et al. 2013) to characterize the metabolic states of germinating rice seeds under air and anoxia. The gene -protein reaction (GPR) mappings were updated for all reactions since the original model did not include the unique genes associated with enzymes having multiple isozymes across different cellular compartments. The updated model is provided in Supporting Information.

Microarray data
The raw transcriptome data generated by Lasanthi-Kudahettige et al. (2007) were first downloaded from the Gene Expression Omnibus (accession no. GSE6908) and normalized using the quantile method (Bolstad et al. 2003). Differentially expressed genes were then identified by a linear model (Wettenhall and Smyth 2004). The resulting P values were also corrected for multiple testing using Benjamini-Hochberg correction.

Random flux sampling
The artificial centring hit-and-run (ACHR) Monte Carlo sampling (Thiele et al. 2005;Mo et al. 2009) was utilized to sample uniformly the metabolic flux solution space under aerobic and anaerobic conditions with appropriate flux constraints. Under both conditions, we constrained sucrose uptake rate, oxygen uptake rate and cell growth rate with experimentally measured values (Lakshmanan et al. 2013). To make a fair comparison between both conditions, oxygen uptake and growth rates were normalized with respect to sucrose uptake rate. The solution space was sampled with 100 000 randomly distributed points for 10 000 iterations in each simulation. In this study, the COBRA toolbox (Schellenberger et al. 2011) was utilized to implement the random flux sampling.
The differences in flux samples between aerobic and anaerobic conditions were quantified using a Z-score approach as described previously (Mo et al. 2009). In this approach, two random flux vectors, v j , one from each sample, i.e. aerobic and anaerobic, were chosen and the difference is calculated as follows: this flux difference sample, the sample mean, m j , and standard deviation, s j , were computed to calculate the Z-score as follows: Finally, the absolute Z-scores were translated to P values using the normal cumulative distribution function and the reactions with P values ,0.05 were deemed as statistically different between aerobic and anaerobic conditions.

Identification of transcriptionally regulated enzymes
Using the P values of transcriptome data and flux sampling, we identified the reactions that are transcriptionally regulated (Bordel et al. 2010). Briefly, if the flux and gene expression (for both up-and down-regulated) significantly changes in the same direction, the corresponding enzyme is classified as 'transcriptionally regulated'.
On the other hand, if the values significantly change in the opposite direction, the enzyme is classified as 'metabolically regulated'. In the case of reactions with multiple isozymes, the gene expression was considered to be upor down-regulated based on the expression values of the majority of the transcripts. For example, if a gene has more up-regulated transcripts than down-regulated ones, then it is considered as up-regulated.

Motif detection and identification of putative TFs
Promoter sequences [21000, +200 nt relative to the transcription start site (TSS)] for transcriptionally up-and down-regulated genes related to anaerobic central metabolism were extracted from our in-house rice promoter sequence database. Known and novel promoter motifs were detected using the Dragon Motif Builder program with EM2 option (Huang et al. 2005

Random sampling reveals significant differences in rice central metabolism between aerobic and anaerobic conditions
We sampled the plausible metabolic states of rice coleoptile in air and under anoxia using ACHR Monte Carlo sampling to estimate the range of possible steady-state flux values through each of the reactions in the rice model (see Methods). The resulting probability distributions of individual reaction fluxes revealed significant differences in central metabolic pathways such as glycolysis, the tricarboxylic acid (TCA) cycle, the pentose phosphate pathway and oxidative phosphorylation between aerobic and anaerobic conditions [see Supporting Information]. Under anoxia, the TCA cycle, the pentose phosphate pathway and oxidative phosphorylation had only a small range of fluxes, with some even having zero flux due to the imposed capacity constraints, i.e. the absence of oxygen exchange [see Supporting Information]. As a result, the fluxes across various amino acid and lipid synthetic pathways were also severely restricted under anoxic conditions (data not shown). On the other hand, random sampling allowed us to observe the possibility of high fluxes through cytosolic glycolysis and fermentation, mainly to produce all the ATP required for anaerobic cell growth. Other energyproducing pathways such as oxidative phosphorylation and the TCA cycle are grossly impaired under such conditions (Mohanty et al. 1993;Lakshmanan et al. 2013).

Transcriptionally regulated reactions during anaerobic adaptation
We compared the differences in flux samples and gene expression data between air and anoxia, thereby identifying the reactions that are likely to be transcriptionally regulated (see Methods). Overall, among the 63 reactions in rice central metabolism, 37 and 5 exhibit transcriptional and metabolic regulation, respectively. The remaining 38 reactions could not be classified in any of these categories as they exhibited an insignificant change in either flux or gene expression. The complete lists of transcriptionally and metabolically regulated reactions are provided in Supporting Information.
All the reactions in the fermentation pathway including pyruvate decarboxylase (PDC), alcohol dehydrogenase (ADH) and aldehyde dehydrogenase (ALDH) were upregulated ( Fig. 1). In rice, both alcohol and aldehyde dehydrogenases have multiple isozymes and a few of them such as ADH1, ADH2 and ALDH2a showed a drastic increase in gene expression (Lasanthi-Kudahettige et al. 2007). Interestingly, the up-regulation of ALDH at both flux and transcriptional level is in very good agreement with earlier experimental reports on both tolerant and intolerant lines of rice under hypoxic conditions during submergence (Nakazono et al. 2000;Ismail et al. 2012;Miro et al. 2013), indicating that ALDH is likely to play a key role in detoxifying the excess acetaldehyde from PDC that cannot be metabolized via ADH. In this regard, Lasanthi-Kudahettige et al. (2007) hypothesized that the acetate resulting from ALDH can enter the TCA cycle to further fuel amino acid synthesis pathways. However, our flux analysis revealed that the acetate from ALDH primarily fuels lipid biosynthesis in plastids, as fatty acid synthesis needs sufficient carbon flux to synthesize the required saturated fatty acids under anoxia.
Unlike other central metabolic pathways, the glycolytic reactions did not show an overall up-or down-regulation in gene transcription and fluxes (Fig. 1). Out of a total of 25 reactions, only 11 showed transcriptional regulation. Among the remaining 14, 5 were oppositely correlated between flux and gene expression, indicating that these enzymes are most likely to be metabolically regulated. Plastidic pyruvate kinase, UDP-glucose pyrophosphorylase, PPi-dependent phosphofructokinase (PFP), and cytosolic and plastidic pyruvate phosphate dikinase (PPDK) are the five enzymes with metabolic regulation. Such observations are in good agreement with the earlier hypothesis by Plaxton (1996), who suggested that PK, PPDK, PFK and PFP are involved in the fine control of plant glycolysis and are most likely to be regulated based on the variation in substrate(s) and cofactor concentrations, pH and other metabolite effectors.

TFs associated with transcriptionally regulated genes
To identify the potential TFs involved in the transcriptional control of the transcriptionally regulated reactions, we performed promoter analysis of the corresponding genes. In total, promoter sequences for 26 transcriptionally up-regulated genes and 27 down-regulated genes were used for cis-element detection [see Supporting Information]. This analysis identified several highly enriched putative cis-elements associated with different potential TFs for both up-and down-regulated transcriptionally controlled genes (Tables 1 and 2). A high enrichment of putative cis-elements such as AT-hook/PE-1-like, GT element-like, pyrimidine-box-like, GARE-like and MYb-box-like associated with MYB TFs was observed among both up-and down-regulated genes. However, the high enrichment of these cis-elements in the upregulated genes signifies that MYB TFs especially play an important role in the transcriptional control of up-regulated genes under anoxia. Similarly, putative cis-elements such as AS-1/ocs-like and ABRE-like associated with bZIP TFs were also found in both up-and down-regulated genes but with a high percentage of occurrences only in up-regulated genes. Besides, a number of other putative cis-elements such as ERE-like/GCC-boxlike and zinc finger binding element-like associated with ERF and ZnF TFs were also more overrepresented in the up-regulated genes albeit being present in both up-and down-regulated genes. Furthermore, we noticed a moderate enrichment of other ERE like-elements such as JA response element-like elements among the upregulated genes (Table 1). Collectively, these results indicate that all these TFs, i.e. MYB, bZIP, ERF and ZnF, may work together in the transcriptional control of the upregulated genes in rice central metabolism. Interestingly, we also noticed the enrichment of MBF1C binding element and CRT/DRE-like elements associated with MBF1C and CBF/DREB TFs only in down-regulated genes. This highlights the possibility that this TF could be involved in transcriptional control of down-regulated genes under anoxia (Table 2). pyrophosphate; Q, ubiquinone; QH2, ubiquinol; R-5-P, ribose-5-phosphate; Ru-5-P, ribulose-5-phosphate; UDP-G, UDP-glucose; X-5-P,

Motif detection in negative sets
To ensure that our cis-element analysis did not identify an excessive number of false positives, we sought motifs from the promoters of a similar number of genes that are not anoxia specific using the same protocol (see Methods). For this purpose, we used two negative datasets: (i) randomly selected genes from the non-differentially expressed gene list under anoxia (Lasanthi-Kudahettige et al. 2007) (negative Set 1) and (ii) up-regulated genes associated with drought response in rice (Zhang et al. 2012) (negative Set 2). The list of genes from the negative sets and the results of this analysis are provided in Supplementary Information. Using the negative control Set 1, we identified a few common motifs such as MYB, bZIP and DOF, and their total enrichment was much lower compared with the motifs detected in transcriptionally up-or down-regulated genes of anoxia (Tables 1, 2 and Supplementary Information). Similarly, the motif analysis of negative Set 2 also revealed a significantly different cis-element enrichment pattern where most of the identified TFs such as MYB, bZIP, ERF, NAC and MYC have been experimentally confirmed to play a regulatory role in drought stress response (Shinozaki and Yamaguchi-Shinozaki 2007;Zhang et al. 2012). Collectively, these results clearly demonstrate that our cis-element method identifies reasonably precise motifs that are specific to anoxic stress.

Discussion
Our study analysed the plausible metabolic states of rice under aerobic and anaerobic conditions using random sampling and compared their differences. We also compared the changes in flux levels of individual enzymes in the central metabolism of rice grown in the two conditions with that of gene expression and identified the transcriptionally controlled reactions. Accordingly, the major contribution of the current work is the demonstration of the utility of random flux sampling and gene expression data for elucidating metabolic differences and the identification of corresponding putative TFs responsible for such changes under the stress conditions.
Monte Carlo flux sampling has been shown to be an effective tool to analyse and compare the plausible solution space of a metabolic network across different environmental/genetic conditions (Famili et al. 2003;Bordel et al. 2010;Osterlund et al. 2013). When compared with the commonly used constraint-based analysis methods such as flux balance analysis (FBA) and minimization of metabolic adjustment, which predict just a single flux distribution based on a particular objective function, random sampling allows us to account for uncertainty in the flux distributions as it uniformly samples the possible solution space in an unbiased manner and, thus, potentially minimizes the possibility of overestimating the actual flux differences between the two conditions. This can be clearly exemplified by comparing the results of previous work (Lakshmanan et al. 2013), which utilized FBA for simulating flux differences between air and anoxia, with our current results. Although both studies show that most of the central metabolic reactions possess different flux levels in aerobic and anaerobic conditions, flux sampling can provide confidence scores to the reactions, thus allowing us to consider the reactions with the most plausible fluxes. Furthermore, by comparing the significant flux differences with gene expression data, we determined the reactions that are most probably regulated at the transcriptional level. Most importantly, we revealed that although the glycolytic reactions are not transcriptionally regulated, genes coding for reactions upstream and downstream of glycolysis, i.e. sucrose metabolism and fermentation, respectively, are regulated transcriptionally. It will be important now to compare and contrast such findings with anoxia-intolerant higher plants such as Arabidopsis thaliana, wheat and maize to analyse whether the corresponding gene orthologues are commonly regulated at a similar level.
In recent years, in silico promoter analysis of differentially expressed genes from microarray data has improved our understanding of the transcriptional regulation of gene expression markedly. Several in silico analyses on the promoter cis-elements and their cognate TFs have been experimentally verified for predicted functions, highlighting its reliability (Meier et al. 2008;Yun et al. 2010 Stamm et al. 2012). Therefore, in this study, we analysed the cis-regulatory content of transcriptionally controlled reactions through such a method and unravelled the possible involvement of a number of TFs in transcriptional control of sucrose metabolism and fermentation under anoxia. Among all the putative cis-elements identified, the elements associated with MYB family TFs were found to be highly enriched in transcriptionally up-regulated genes. Coincident with these observations a number of MYB family genes are also up-regulated under anoxia (Table 3). MYB proteins represent important plant TFs and are found to be involved in various developmental and physiological processes including transcriptional activation, kinase activity, protein binding and transcription repressor activation under abiotic and biotic stresses (Dubos et al. 2010). In rice, a recent analysis highlighted that 98.70 % of total MYB proteins are fully involved in transcriptional activation (Katiyar et al. 2012). Furthermore, the role of MYB TFs that bind specifically to MYB box/GTelement during hypoxic and anoxic responses has already been reported through promoter analysis, both computationally and experimentally (Dolferus et al. 2003 2005, 2012). Specifically, the AtMYB2 in Arabidopsis has been shown to be a key regulator of the ADH1 promoter under low oxygen conditions (Hoeren et al. 1998). When this AtMYB2 was driven by a constitutive promoter, it was able to transactivate ADH1 expression not only in Arabidopsis but also in Nicotiana plumbaginifolia and Pisum sativum. Collectively, these results fully support the hypothesis that MYB TFs play an important role in the upregulation of sucrose metabolism and fermentation enzymes at the transcriptional level. As mentioned above, our promoter analysis also highlighted significant overrepresentation of several other motifs associated with bZIP, ERF and ZnF TFs in the transcriptionally up-regulated genes. In this regard, an increase in transcript levels of bZIP TF (AtbZIP50) in anoxia-exposed root cultures of Arabidopsis (Klok et al. 2002), the up-regulation of an ABRE-binding bZIP TF (OsABF1) in rice shoot and root under anoxic treatment (Amir Hossain et al. 2010) and the up-regulation of a number of bZIP TFs in rice anoxic coleoptile (Table 3) support the view that bZIP could orchestrate the transcription of up-regulated genes upon anoxic stress. The major role of ERF TF has been identified as a positive regulator of Sub1A expression in a flood-tolerant rice variety during hypoxic conditions caused by submergence (Xu et al. 2006 Table 3. List of anoxia-stressed up-regulated TFs with potential significance to the pattern of cis-element enrichment among the up-regulated genes.

Family
Locus_ID ( 'Nipponbare' that germinates and elongates under anoxia does not have this gene. Surprisingly, the presence of highly enriched putative GCC-box-like/ERE-like cis-elements in the promoters of all anoxia up-regulated genes in this variety of rice highlights the possiblity that these TFs might also control the transcription of relevant genes as it is in Arabidopsis under hypoxic conditions (Hinz et al. 2010;Licausi et al. 2010Licausi et al. , 2011. Unlike the ERF TFs, the role of ZnF TF in response to submergence/anoxia is not yet known. However, matching with our identification of ZnF binding site-like elements in transcriptionally up-regulated genes, the increase in expression of ZnF TF genes (Table 3) in response to hypoxia/anoxia in both rice and Arabidopsis allows us to speculate on its potential regulatory role in transcriptional control under oxygen stress (Loreti et al. 2005;Pandey and Kim 2012).
As mentioned earlier, we utilized the promoter sequences of 26 transcriptionally up-regulated and 27 down-regulated genes for the cis-element detection in this study. When compared with the normal promoter analysis method that considers the whole set of differentially expressed genes, the current study takes into account relatively few genes. Therefore, in order to argue that this small set of genes is sufficient to generate biologically meaningful results, we first compared the TFs identified in this study with previous work that considered the whole list of up-regulated (842) and down-regulated (1794) genes during anoxic adaptation ). Interestingly, the comparison result revealed that the TFs associated with up-and down-regulated genes are the same between two studies where the motif enrichment scores of certain TFs such as MYB, bZIP, ZnF and ERF are much higher in the current study, highlighting the possible role of these TFs in anoxic adaptation. To further confirm our findings, we also analysed the individual promoter sequences of four of the key transcriptionally regulated enzymes, i.e. SUS, PDC, ADH and ALDH, and identified the same TFs as fully responsible for their transcription (Fig. 2). Therefore, based on these motif analysis results and experimental evidence from the literature, we hypothesize a positive role of MYB together with bZIP, ERF and ZnF in the transcriptional control of sucrose metabolism and fermentation during germination and coleoptile elongation of rice under anoxia.
The current methods we used successfully identified several transcriptionally regulated reactions and their related TFs, some of which are experimentally confirmed. However, the overall results of this approach await experimental validation since the current prediction relies on several assumptions concerning model completeness, constraints used during simulation and the statistical cutoff values chosen for comparative analysis. In our work, we utilized the central model to simulate the differences in metabolic fluxes between air and anoxia. It should be noted that although this model predicts the overall Figure 2. Presence of common putative cis-elements in the key transcriptionally regulated genes. The presence of potential putative cis-elements and their cognate known TFs are shown in different strands of the promoter region (21000, +200 nt relative to TSS) of the key genes. Each putative cis-element/motif is represented by its consensus logo. TATA boxes are located between 25 and 30 nt upstream from the TSS. cellular phenotype quite accurately (Lakshmanan et al. 2013), it may not capture the global changes in cellular metabolism. This reservation is needed since the model does not take account of the secondary metabolic pathways. Moreover, we have utilized only a few external fluxes, i.e. sucrose uptake, O 2 uptake and growth yield, as the sole constraints to identify the internal flux distributions during flux sampling simulations. In such cases, the internal fluxes of some reactions are determined with statistically low confidence scores due to the possibility of multiple flux solutions. Such limitations can be overcome by the use of C 13 flux measurements as constraints to the internal reactions. Therefore, the list of transcriptionally regulated reactions and the TFs identified in this study need further confirmation by direct experiments.

Conclusions
In this study, we presented a combined in silico modelling and gene expression data analysis framework to identify transcriptionally regulated reactions in rice central metabolism. Random flux sampling simulations highlighted significant changes in flux levels across oxidative phosphorylation, the TCA cycle, the pentose phosphate pathway, glycolysis, sucrose metabolism and fermentation. The subsequent comparative analysis of changes in flux levels and gene expression between aerobic conditions identified 37 reactions from these pathways that are regulated at the transcriptional level. The motif enrichment analysis of transcriptionally regulated enzymes highlighted the potential involvement of TFs MYB, bZIP, ERF and ZnF in controlling the transcription of sucrose metabolism and fermentation genes under anaerobic conditions. Thus, the current study successfully demonstrated how condition-specific gene expression data can be exploited to narrow down the metabolic changes in rice under anaerobic adaptations and in silico analysis to help us unravel the relevant transcriptional mechanisms.

Supporting Information
The following Supporting Information is available in the online version of this article -File S1. Table. Rice central model with updated GPR associations.
File S2. Diagram. Sampling histogram of glycolysis, pentose phosphate pathway, TCA cycle, oxidative phosphorylation and sucrose metabolism under aerobic and anaerobic conditions. File S3. Table. List of rice central metabolic reactions and their flux differences with P values. Expression patterns of corresponding gene louses and type of regulation, i.e. transcriptional or metabolic, are also listed.
File S4. Table. List of genes used in negative sets. File S5.