Deﬁning Core Metabolic and Transcriptomic Responses to Oxygen Availability in Rice Embryos and Young Seedlings

Analysis reveals that there is limited overlap in the sets of transcripts that show signiﬁcant changes in abundance during anaerobiosis in different plant species. This may be due to the fact that a combination of primary effects, changes due to the presence or absence of oxygen, and secondary effects, responses to primary changes or tissue and developmental responses, are measured together and not differentiated from each other. In order to dissect out these responses, the effect of the presence or absence of oxygen was investigated using three different experimental designs using rice ( Oryza sativa ) as a model system. A total of 110 metabolites and 9,596 transcripts were found to change signiﬁcantly in response to oxygen availability in at least one experiment. However, only one-quarter of these showed complementary responses to oxygen in all three experiments, allowing the core response to oxygen availability to be deﬁned. A total of 10 metabolites and 1,136 genes could be deﬁned as aerobic responders (up-regulated in the presence of oxygen and down-regulated in its absence), and 13 metabolites and 730 genes could be deﬁned as anaerobic responders (up-regulated in the absence of oxygen and down-regulated in its presence). Deﬁning core sets of transcripts that were sensitive to oxygen provided insights into alterations in metabolism, speciﬁcally carbohydrate and lipid metabolism and the putative regulatory mechanisms that allow rice to grow under anaerobic conditions. Transcript abundance of a speciﬁc set of transcription factors was sensitive to oxygen availability during all of the different experiments conducted, putatively identifying primary regulators of gene expression under anaerobic conditions. Combined with the possibility of selective transcript degradation, these transcriptional processes are involved in the core response of rice to anaerobiosis.

Various species, including some bacteria, yeasts, and plants, have the ability to survive and grow in the absence of oxygen (Bunn and Poyton, 1996;Kwast et al., 1998;Bailey-Serres and Voesenek, 2008). Although the ability of plants to germinate and survive under anaerobic conditions varies greatly, one common response is the synthesis of a number of anaerobic proteins, such as alcohol dehydrogenase, which was identified almost 30 years ago (Sachs et al., 1980) and has been extensively studied at transcript, protein, activity, and regulatory levels (Dolferus et al., 2003;Bailey-Serres and Voesenek, 2008). Analysis of the effects of anaerobic conditions on various plant species reveals that widespread changes limiting energydemanding processes is a common response (Bailey-Serres and Voesenek, 2008). Overall, under conditions where oxygen is limited, an increase in glycolytic flux, modified Suc degradation, and a decrease in processes such as transport and various biosynthetic pathways occur. More specifically, a switch to the pyrophosphate (PPi)-linked enzymes, such as PPi-dependent phosphofructokinase and pyruvate orthophosphate dikinase, as well as the breakdown of Suc via Suc synthase conserves ATP, a response reported in a variety of studies Magneschi and Perata, 2009).
Despite differences in tolerance to anaerobic conditions, plant size, and life cycle, the use of transcriptomic approaches to study the response to anaerobic conditions in Arabidopsis (Arabidopsis thaliana; Liu et al., 2005;Loreti et al., 2005;Van Dongen et al., 2009), rice (Oryza sativa; Lasanthi-Kudahettige et al., 2007;Magneschi and Perata, 2009), and poplar (Populus 3 canescens; Kreuzwieser et al., 2009) reveals similar changes in primary carbon metabolism. However, more global analysis of these studies also reveals that the number of genes that displayed significant changes in abundance varied greatly between species; in Arabidopsis, only 1,266 transcripts were altered in abundance in response to anoxia (Klok et al., 2002;Liu et al., 2005), 3,134 probe sets were observed to change sig-nificantly in rice coleoptiles under the conditions analyzed (Lasanthi-Kudahettige et al., 2007), while as many as 5,000 transcripts were altered in roots of poplar (Kreuzwieser et al., 2009). Additionally, analysis of the changes between the anaerobiosis-tolerant plants, rice and poplar, from these two studies revealed that less than 5% (107) of the transcript changes that occurred in poplar were shared by the orthologous rice transcripts. Given that the species studied have different genome sizes and that varying experimental conditions were used, it is not clear whether or not the observations indicate substantial differences in the primary responses of these species to a lack of oxygen.
Rice plants can survive periods of anaerobiosis, and rice seeds also have the ability to germinate in an environment completely lacking in oxygen. The identification of the SUB1A locus provides an important breakthrough for understanding molecular responses of rice under submergence (Xu et al., 2006), but this differs from germination under anaerobic conditions, where there is no oxygen present to synthesize ethylene or many other hormones (Magneschi and Perata, 2009), and cultivars lacking the SUB1A gene successfully germinate under anaerobic conditions (Magneschi and Perata, 2009, and refs. therein). In efforts to increase tolerance to anaerobic conditions, a number of genes encoding components involved in glycolysis, fermentation, or sugar metabolism have been overexpressed in order to increase survival, with varying degrees of success (Christianson et al., 2009). An alternative approach has been the identification of regulatory factors that control expression of multiple genes involved in tolerance to anaerobic conditions. This type of approach has been undertaken to increase tolerance to a wide variety of stresses (Century et al., 2008), although as yet there are no commercial varieties produced (Arce et al., 2008). The expression of a NAC transcription factor from Arabidopsis has been shown to increase tolerance to low-oxygen treatments (Christianson et al., 2009).
In order to gain greater insight into the regulatory processes that modulate responses to anaerobic conditions, three experimental strategies were used to define oxygen-responsive genes in rice. It was hypothesized that many of the changes that occur are not primarily due to the effects of a lack of oxygen but instead are dependent on prior growth conditions or secondary responses. We have previously characterized changes in water content and metabolic activity in rice embryos during germination up to 48 h after imbibition (HAI) under aerobic and anaerobic conditions (Howell et al., 2006. We have used this system to investigate the temporal changes in metabolite and transcript abundance under three different experimental conditions: during germination under anaerobic conditions, highlighting the differences compared with germination under aerobic conditions , and we also switched young rice seedlings (24 HAI) grown under aerobic or anaerobic conditions to the other condition and determined changes in transcript and metabolite abundance (Fig.  1A). These three experimental conditions were used in combination to define core aerobic and anaerobic genes (Fig. 1A).

Changes in Transcript and Metabolite Abundance during Germination in Aerobic or Anaerobic Conditions
An analysis of changes in transcript abundance during germination between aerobic and anaerobic conditions revealed that there was little difference up to 3 HAI; as a result, transcriptomic data for aerobic and anaerobic germinated seeds at this time point were not differentiated by principal component analysis (PCA) analysis ( Fig. 1B; Supplemental Fig. S1). As approximately 2,000 genes had significantly altered transcript abundances 3 HAI under both growth conditions when compared with 1 HAI (Fig. 1B, 1A v 3A and 1N v 3N; Supplemental Tables S1 and S2) and greater than half of these changes were significant increases in abundance, this indicated that initial changes in transcript abundance occurring 3 HAI were independent of oxygen levels. In contrast, after 3 HAI, 2,628 (1,148 up-regulated, 1,480 down-regulated) and 4,892 (2,195 up-regulated, 2,697 down-regulated) transcripts were observed to have significantly different levels at 12 and 24 HAI, respectively (Fig. 1B,A v N), with the different samples forming distinct groups as visualized by PCA analysis (Supplemental Fig. S1). In total, there were 5,948 nonredundant probe sets that showed significant differences in transcript abundance between the two growth conditions up to 24 HAI, compared with approximately 14,000 or approximately 12,000 that were observed to change between 0 and 24 HAI under aerobic  or anaerobic growth conditions, respectively (Supplemental Tables S1 and S2; Supplemental Fig. S2). Differences between tissues grown for 12 h under aerobic and anaerobic conditions are likely to represent a primary response, given that at this time point developmental/tissue-specific differences are not apparent . However, at later time points (24 and 48 HAI), there are likely to be secondary changes reflecting differences in morphology, given that trace amounts of oxygen are required for root growth development .
Analysis of the functional categorization of the genes for the 5,948 transcripts that were higher in abundance in anaerobic germinated seeds (A v N) revealed that lipid metabolism and transcription factor categories were overrepresented based on z-scores where P , 0.01 (Fig. 1C, i; Supplemental Table S3). In contrast, functional categorization of the genes for transcripts that were lower in abundance under anaerobic conditions was enriched in transcripts encoding metabolism, carbohydrate and energy metabolism, and membrane transport functions, with transcripts encoding transcription factors and associated with translation and metabolic/genetic information processing being underrepresented (Fig. 1C, ii).
An analysis of changes in metabolites during germination under aerobic and anaerobic conditions revealed that the number of metabolites that changed was similar up to 12 HAI (Fig. 1B,3A v 12A and 3N v 12N). However, comparing metabolite profiles from aerobic and anaerobic growth conditions showed a number of significant changes, particularly after the 3-h time point, similar to the pattern seen for the transcripts. Specifically, more than 20% of all detected metabolites showed significant differences when 12-HAI samples were compared (Fig. 1B, 12A v 12N), and this increased to more than 40% at the 48-h time point (48A v 48N). Examples included several amino acids such as Gly and Tyr that had significantly (P , 0.05) higher levels under anaerobic germination, while the several metabolites associated with the tricarboxylic acid (TCA) cycle, such as isocitrate, 2-oxoglutarate, and citrate, had significantly (P , 0.05) lower levels under anaerobiosis (Supplemental Table S5).

Defining Common Changes in Response to Oxygen Availability
As there are developmental differences between aerobic and anaerobic germinated seeds , two additional experimental designs were used that monitored the response to switching from aerobic to anaerobic conditions and vice versa at 24 HAI (Fig. 1A). This allowed complementary responses in all three experimental conditions to be defined as aerobic or anaerobic responses. Thus, aerobic transcripts and metabolites were defined as (1) significantly down-regulated during anaerobic germination (A v N), (2) significantly down-regulated in response to switching to anaerobic conditions (A / N), and (3) significantly up-regulated in response to switching to aerobic conditions (N / A). Likewise, anaerobic transcripts and metabolites were defined as those that are (1) significantly up-regulated during anaerobic germination (A v N), (2) significantly upregulated in response to switching to anaerobic conditions (A / N), and (3) significantly down-regulated in response to switching to aerobic conditions (N / A). Analysis revealed that 10 detected metabolites (six of which could be identified) could be defined as aerobic while 13 (11 identified) could be defined as anaerobic (Fig. 2).
Applying this principle of seeking complementary responses between three experiments to the transcriptomic analysis defined 1,136 and 730 transcripts as aerobic and anaerobic, respectively (Fig. 3). Functional analysis of the aerobic set (1,136) revealed significant enrichment of metabolism and membrane transport functions and underrepresentation of transcription factors, translation, folding, sorting, degradation, and replication and repair functions (Fig. 3). For the anaerobic set, there was significant enrichment of carbohydrate and lipid metabolism and underrepresentation of translation functions (Fig. 3).

Linking Metabolite and Transcript Changes to Visualize Modification of Metabolism under Anaerobic Conditions
A comparison of the differences in transcripts after 24 HAI under aerobic or anaerobic conditions, or upon transferring seedlings from growth in anaerobic conditions to aerobic conditions, revealed differences in lipid metabolism, cell wall metabolism, secondary metabolism (particularly flavonoids and phenylpropanoids and phenolics), carbohydrate metabolism, and amino acid metabolism (Supplemental Fig. S3). A variety of metabolite and transcript changes were observed that relate to starch/Suc breakdown via glycolysis, the TCA cycle, g-aminobutyrate (GABA) shunt, and amino acid metabolism (Fig. 4). These included a transcript encoding cytosolic pyruvate orthophosphate dikinase (LOC_Os03g31750.1), which was induced greater than 100-fold, in contrast to the transcript level of the plastidic isoform (LOC_ Os05g33570.1), which was reduced greater than 26fold during anaerobic germination. Several other Figure 1. Experimental design and summary of changes in transcript and metabolite abundance in rice in response to oxygen. A, Overview of the experimental design. Rice seeds were imbibed in liquid medium under aerobic or anaerobic conditions. Three biological replicates for analysis of the transcriptome were taken at 0, 1, 3, 12, and 24 HAI under both growth conditions (the 0-h time point was common to both treatments) to give nine sample sets where genome-wide changes in transcript abundance were analyzed. Additionally, after 24 h of growth, samples were switched, where aerobically grown samples were switched to anaerobic conditions and samples taken 3 and 6 h after switching (27 and 30 h after imbibition, respectively). Anaerobically germinated seeds were switched to aerobic conditions, and samples were taken at the same time points. This yielded eight data sets where genome-wide changes in transcript abundance were analyzed. In gray, time points for metabolite profiling analysis are indicated. Three additional time points were included. Red indicates aerobic conditions, and blue indicates anaerobic conditions. B, Summary of the changes in transcript and metabolite abundance during the experiments outlined in A. The number of probe sets significantly changing in abundance (left axis, dark color bars) and significant changes in metabolite levels (right axis, pale color bars) are indicated. The time points compared are indicated on the y axis. nd, Not determined. C, Functional categorization of the encoded proteins of the 5,944 transcripts observed to be significantly (P , 0.01) changing in abundance over germination (A v N). Genome-wide distribution of the functional categories of proteins and those that are overrepresented (shown in red) or underrepresented (shown in blue; based on z-scores with P , 0.01) for transcripts up-regulated (i) or down-regulated (ii) during anaerobic germination is shown.
transcripts showed large increases in abundance during anoxia, including those encoding pyruvate decarboxylase (LOC_Os05g39320.1), hexokinase (LOC_Os05g09500.1), and L-allo-Thr aldolase (LOC_ Os04g43650.1). A number of transcripts encoding TCA cycle components were found to decrease, including citrate synthase (LOC_Os02g10070.1), aconitase (LOC_Os03g04410.1), isocitrate dehydrogenase (LOC_Os04g40310.1), 2-oxoglutarate dehydrogenase (LOC_Os04g32020.1 and LOC_Os07g49520.1), and succinyl-CoA ligase (LOC_Os07g38970.1 and LOC_ Os02g40830.1). These were accompanied by decreases in associated metabolites, citrate, isocitrate, and, to the greatest extent, 2-oxoglutarate (Fig. 4). In contrast, a 16-fold increase was observed for succinate, and fumarate increased more than 2-fold. The simplest explanation for succinate accumulation is a block in its catabolism via the respiratory chain component succinate dehydrogenase due to a lack of oxygen. However, the observed increase in succinate also may have been due to enhanced synthesis from GABA, via the well-described pathway of Glu decarboxylase, GABA transaminase, and succinic semialdehyde dehydrogenase (Fig. 4;Fait et al., 2008). It has already been proposed that the production of GABA is due to the acidification of the cytosol under anoxic conditions, which stimulates the activity of Glu decarboxylase (Ratcliffe, 1995). GABA can then be converted into succinate semialdehyde via transamination with pyruvate to produce Ala; the latter was also observed to increase 16-fold at 48 HAI (Figs. 2 and 4). Both transcript and metabolite profiles suggested major reprogramming of pathways feeding into and out of the aromatic amino acids, Phe and Tyr. Transcriptional down-regulation of two enzymes encoding phospho-2dehydro-3-deoxyheptonate aldolase (LOC_Os03g27230.1, 4.2-fold down-regulated; LOC_Os07g42960.1, 2-fold down-regulated), which controls entry of phosphoenolpyruvate and erythrose-4-phosphate into the shikimate pathway (Herrmann, 1995), together with a 5-fold down-regulation of the bifunctional enzyme responsible for shikimate synthesis, 3-dehydroquinate dehydratase/shikimate dehydrogenase (LOC_ Os01g27750.1), were accompanied by 10-fold lower levels of shikimate at 24 HAI. It appears that many pathways downstream of shikimate were downregulated at a transcript level (i.e. chorismate, Trp, Phe, and Tyr synthesis as well as those associated with flavonoids, phenylpropanoids, and lignins; Fig. 4; Supplemental Fig. S3). Despite the apparent downregulation of aromatic biosynthetic pathways, increases were seen in a variety of aromatic metabolites, including the amino acids Phe and Tyr and the phenylpropanoids coumarate (4-hydroxycinnamate) and sinapinate (3,5dimethoxy-4-hydroxycinnamate), which are precursors for lignin biosynthesis.

Using Germination and Switch Experiments to Explore the Oxygen Dependence of Specific Processes
Lipid metabolism is known to be sensitive to oxygen at three distinct steps: in b-oxidation of fatty acids at the acyl-CoA oxidase step (Graham, 2008), in desaturation of fatty acids, and in sterol synthesis (Buchanan et al., 2002). Several transcripts involved in lipid metabolism appear to be affected by oxygen in that they were up-regulated in anaerobically grown rice compared with aerobically grown rice (Supplemental Fig. S4, 24 A v 24 N). There was a significant enrichment of transcripts up-regulated under anaerobic germination and upon switching to anaerobic conditions (A to N; i.e. 3.47% and 3.61%, respectively; Supplemental Fig. S4; Supplemental Table S3). In total, 13 (4.08%) of the anaerobiosis-responsive subset encoded genes involved in lipid metabolism (Supplemental Fig. S4, in boldface). Notably, this included genes encoding acyl-CoA ligase (LOC_Os03g03790.1 and LOC_Os04g58710.1), triacylglycerol lipases, and several transcripts encoding proteins involved in sphingolipid metabolism, which were up-regulated under anaerobic conditions (Supplemental Fig. S4). Interestingly, it was found that the anaerobic subset also contained several other transcripts encoding proteins involved in lipid metabolism-related functions (e.g. those classified as "lipid metabolism/signal transduction," such as the transcript encoding an acylactivating enzyme (LOC_Os03g03790.1), which was up-regulated over 50-fold upon switching to anaerobic conditions (Supplemental Fig. S4). Finally, with respect to lipid metabolism, rice contains 28 genes encoding putative transcription factors that contain the lipid/sterol-binding StAR-related lipid transfer protein domains (START; Schrick et al., 2004;Supplemental Table S6). One (LOC_Os6g02590.1) was defined as a core aerobic gene, and one (LOC_Os07g08760.1) was defined as a core anaerobic gene, implicating a regulatory role for changes in lipid metabolism under aerobic and anaerobic conditions. Defining core sets of aerobic and anaerobic metabolites in rice. A, Metabolites found to have significant changes in abundance under the three different experimental conditions were compared to determine those metabolites that were consistently responsive to oxygen availability and could be defined as aerobic (i) or anaerobic (ii). B, A heat map showing the fold changes over germination and both switch conditions, sorted to highlight the aerobic (red boxes) and anaerobic (blue boxes) metabolites that could be reliably identified. Red color indicates a significant increase, and blue color indicates a significant decrease. The full list of fold changes for all metabolites detected and associated P values are shown in Supplemental Table S5. *, Metabolite tentatively identified by matching against the Q_MSRI_ID mass spectral and retention index library available at http://csbdb.mpimp-golm.mpg.de/csbdb/gmd/msri/gmd_msri.html. **, Metabolite tentatively identified by matching against the NIST05 mass spectral library (www.nist.gov). The core lists of aerobic and anaerobic genes defined above allow a broader analysis of the regulatory processes that may determine transcript abundance under aerobic or anaerobic conditions. Twelve subsets of the genes that were defined as aerobic (1,136) or anaerobic (730) were analyzed for the presence of common regulatory elements in their upstream promoter regions, based on their potential to regulate transcript abundance (transcription factors) or on their relative enrichment in the set of aerobic or anaerobic genes relative to the whole genome (Supplemental Table S7A). The subsets chosen were aerobic and anaerobic sets of genes encoding transcription factors, proteins targeted to mitochondria and plastids, proteins involved in membrane transport and carbohy-drate metabolism, and random sets of 50 aerobic or anaerobic genes. The MEME database was used to search for conserved elements in the sequence 1 kb upstream of the translation start site (Supplemental Table S7B). Only elements that occurred in 70% of any subset were used for further analyses (Supplemental Table S7C). These elements were taken and their enrichment in these subsets of genes was compared with the genome (Table I). Furthermore, 23 elements or variations of elements previously characterized to mediate responses to anaerobic conditions were tested for enrichment in the various subsets of genes (Supplemental Table S7C). Twelve unique elements were found to be significantly enriched in the promoter regions of anaerobic genes (Table I, highlighted in blue), two elements were found to be enriched in the promoter regions of aerobic genes (Table I, highlighted in pink, percentages in red), while four elements were . Defining core sets of aerobic and anaerobic genes in rice. A, Venn diagram outlining how the aerobic core set of genes was defined. The number of genes in each group and overlapping groups are shown. For the group of 1,136 genes defined as aerobic, functional categorization analysis was carried out to show which categories were significantly enriched (red) or depleted (blue; based on z-scores with P , 0.01). B, Venn diagram outlining how the anaerobic core set of genes was defined. The number of genes in each group and overlapping groups are shown. For this group of 730 genes defined as anaerobic, functional categorization analysis was carried out to show which categories were significantly enriched (red) or depleted (blue; based on z-scores with P , 0.01).
significantly underrepresented in aerobic genes (Table  I, highlighted in pink, percentages in blue). Fifteen elements were found to be enriched in both anaerobic and aerobic genes (Table I, highlighted in gray; Supplemental Table S7C).
Under the three experimental conditions, approximately 40% to 50% of transcripts were down-regulated in abundance in response to oxygen availability (Fig. 1B). In order to assess the role of transcript degradation, potential regulatory elements were ex- . Parallel display of transcripts and metabolites for starch-Suc metabolism, glycolysis, the TCA cycle, GABA shunt, mitochondrial respiratory chain, and amino acid metabolism. Fold changes in transcripts and metabolites were log transformed and displayed using the MapMan tool. Metabolite changes are indicated in the circles next to the corresponding metabolite name (gray boxes). Metabolite and transcript levels in rice embryos germinated under anaerobic conditions versus aerobic conditions at 24 HAI are compared. For both metabolites and transcripts, changes are represented by shading, where the color saturates at a log 2 fold change value of 3 (i.e. an 8-fold change). Transcripts and metabolites defined as core aerobic or anaerobic are boxed in red or blue, respectively. Note that blue/pink squares are still drawn around some genes/metabolites that appear to be unchanging; these are transcripts/metabolites for which the significant change in abundance occurred at a comparison different from the one visualized (e.g. at 12A v 12N). amined by searching for common sequences within the 3# untranslated regions (UTRs) of genes in the aerobic and anaerobic sets using MEME. It is important to note that the analysis of the 3# UTRs of the aerobic and anaerobic sets was limited by the small number of genes in rice for which this annotation exists. Nevertheless, we did find five enriched elements within the 3# UTRs of the aerobic set (P , 0.01) but none significantly enriched in the anaerobic set. Table I. Comparison of the presence of the putative motifs between the genome and the aerobic and anaerobic subsets The Institute for Genomic Research full genome sequence (Seqs) information, which included the 1-kb upstream regions (66,710) and 3# UTR sequences (3,027), was used for comparison with the genome. The number of sequences (1 kb upstream or 3# UTRs) that contained each motif is shown as a percentage of sequences in that subset. A z-score analysis was carried out (Supplemental Table S7C), and the putative motifs significantly overrepresented/underrepresented (at P , 0.01) are indicated by the red/blue coloring, respectively. The source column indicates the source of the motif (i.e. MEME output from aerobic and/or anaerobic subsets or other studies with known motifs). Shading in blue (anaerobic) or pink (aerobic) indicates that it was significantly enriched or depleted in the respective subsets. Gray shading indicates significant enrichment in both the aerobic and anaerobic responsive sets. Coloring of the motifs indicates overlapping or reverse complement sequences.

Identifying Transcription Factors Changing in Response to Oxygen Availability
From a complete list of 3,098 rice transcription factors (described in "Materials and Methods"), it was found that 2,009 transcripts encoding for transcription factors were detected in at least one sample from the germination and switch samples (Supplemental Table S8). Of these, 437 were found to change significantly in abundance in the germination comparison, and 734 displayed significant differences within the switch samples. Four main subsets were selected for analysis: genes that were significantly upregulated during anaerobic germination (224), genes that were down-regulated in anaerobic germination (212), genes that were present in the aerobic set (63), and genes that were present in the anaerobic set (71; Supplemental Table S8). The distribution of the various transcription factor families within each subset was compared with the genome (2,009) to determine which families were significantly overrepresented or underrepresented in each subset (Supplemental Table  S8). Only four transcription families were significantly overrepresented or underrepresented compared with the genome in one or more subset(s) (Fig. 5). The bZIP family was enriched in the subset of transcripts upregulated under anaerobic germination and in the core anaerobic subset (Fig. 5). The enrichment of the bZIP family in both of these subsets is in agreement with the previous finding that this family is induced in tomato (Solanum lycopersicum) plants subjected to anaerobic conditions (Sell and Hehl, 2004). In contrast, in Arabidopsis, only three bZIP family members were induced under hypoxic conditions, from the 64 transcription factors/kinases that were differentially expressed (Liu et al., 2005). In contrast, the AUX/IAA, bHLH, and ZIM families were all significantly enriched in the subset of genes significantly down-regulated under anaerobic conditions (Fig. 5). Interestingly, it can be seen that bHLH and ZIM transcription factors are significantly overrepresented in the group of genes down-regulated under anaerobic germination but are not overrepresented in the oxygen-responsive set (Fig.  5). This suggests that these transcription factors may be involved at a later stage, such as regulation of secondary responses, or may have a more general role (i.e. not specifically responsive to changes in oxygen availability). In contrast, the AUX/IAA family was also enriched in the oxygen-responsive sets (Fig. 5), indicating a role for these transcription factors in regulating genes during aerobic germination and perhaps in response to oxygen availability.

Aerobic and Anaerobic Transcription Factors Display Contrasting Expression Patterns in Rice Tissues and in Response to Various Stresses
In an attempt to understand how regulation of responses to anaerobic conditions compares with other stresses or rice development and growth, we investigated the expression of the 63 transcription factors defined as aerobic (Fig. 6A) and the 71 transcription factors defined as anaerobic (Fig. 6B) in all publicly available Affymetrix rice microarray data (68 different conditions). Analysis of the expression pattern of the 63 aerobic transcription factors reveals that young root tissue, stigma, and anther have relatively high expression of approximately 15 of these genes ( Fig. 6A; Supplemental Fig. S5A; Supplemental Table   Figure 5. An analysis of the transcription factor families present in each subset. Column 1, up-regulated under germination A v N; column 2, down-regulated under germination A v N; column 3, defined as aerobic; column 4, defined as anaerobic. Only the significantly changing transcription factor families are shown (based on z-score with P , 0.01), with overrepresentation indicated by red asterisks (no families showed significant underrepresentation in any of the subsets). The percentage of transcripts from each subset that were in each transcription factor family is shown, and the significance scores of all transcription factor families are shown in Supplemental Table S8.  . Analysis of aerobic and anaerobic responsive transcription factor genes across publicly available rice arrays. All 63 aerobic and 71 anaerobic responsive genes encoding transcription factors were hierarchically clustered across the germination, switch, and other public Affymetrix arrays (68 conditions/sample types in total). Of these, the 63 aerobic (A) and 71 anaerobic (B) responsive genes are shown across 35 of the array conditions/samples. The full clusters are shown in Supplemental Figure  S4, A and B. S9). This may correspond to energy-demanding processes in these tissues. Expression in response to drought, salt, and cold stress of seedlings was also observed. Expression of these transcription factors was high in salt stress of crown and growing point, and salt-sensitive rice varieties have decreased expression compared with salt-tolerant rice varieties in salt stress experiments ( Fig. 6A; Supplemental Fig. S5A; Supplemental Table S9). Leaf material also appears to have high levels of expression for about 10 aerobic genes that may decrease slightly with cytokinin treatment ( Fig. 6A; Supplemental Fig. S5A; Supplemental Table S9). Analysis of changes in transcript abundance of the 63 aerobic transcription factors in published studies using quantitative reverse transcription-PCR reveal that one NAC transcription factor (LOC_ Os02g13710.1; Fig. 6A, no. 16; Supplemental Table  S9) has been shown to increase under salt stress (Fang et al., 2008), while three IAA response transcription factors (LOC_Os05g08570.1, LOC_Os03g43400.1, and LOC_Os05g14180.1; Fig. 6A, nos. 9, 14, and 63; Supplemental Table S9) decrease under salt stress, and one (LOC_Os03g0742900.1, OsIAA13; Fig. 6A, no. 34; Supplemental Table S9) increases under drought stress (Song et al., 2009). Several Myb transcription factors, including those in the aerobic set, have also been associated with stress response in rice drought stress (Yanhui et al., 2006). Thus, many of the aerobic transcription factors overlap with responses to stress.
For the 71 anaerobic transcription factors, several differences in expression were observed compared with the public array data (  Table S9) has previously been identified to encode a heat shock transcription factor (Yamanouchi et al., 2002), and another (LOC_ Os03g26210.1) increased under dehydration stress (Nijhawan et al., 2008). This group also contains a number of myb-type transcription factors previously associated with stress responses (Yanhui et al., 2006).

DISCUSSION
Using three different experimental conditions, 10 metabolites could be defined as aerobic across all three conditions, out of a total of 39 that displayed a significant difference in the presence of oxygen in any one condition, and 13 were defined as anaerobic in all three conditions, out of a total of 82 that were up-regulated in response to anaerobiosis in any one condition. For the transcriptome, over 9,000 nonredundant transcripts were observed to change in at least one condition. However, only 1,866 of these transcripts showed complementary changes in all three conditions. This reveals that although rearrangement of the transcriptome underlies tolerance to anaerobic conditions in rice, the majority of these changes depend on the prior growth conditions and/or the stage of development tested. The identification of the common responses in all three experimental conditions provides an insight into the core set of genes that respond to oxygen, aerobic (1,136) and anaerobic (730), and provide a basis for detailed analysis of the oxygen regulatory mechanisms.
A previous study analyzed the differences in transcript abundance in rice grown under aerobic and anaerobic conditions (Lasanthi-Kudahettige et al., 2007). Care needs to be taken when making comparisons with this study, as several factors differed, namely the tissue examined (coleoptile versus embryo or young seedling), the means by which differential gene expression was defined, since this study only used two biological replicates, and the period of anaerobic conditions (4 d versus either 24 h or switching growth conditions for up to 6 h). Of the 730 anaerobiosis-responsive genes defined in this study, 681 were present in the coleoptile study (Lasanthi-Kudahettige et al., 2007), and 151 had significant increases in transcript abundance in the anaerobically grown coleoptile. Of the 1,136 oxygen-responsive genes defined in this study, 987 were present in the coleoptile study, and 345 displayed a significant decrease in transcript abundance in the anaerobically grown coleoptile (Supplemental Table S2). Given the large difference in the time period of anaerobic growth between this and the previous rice study, it may indicate that some of the changes observed in this study represent primary changes that may not be observed after a longer period of anaerobiosis, indicating the need for both long-and short-term studies to gain a full insights into the processes that occur during anaerobiosis.
Metabolite and transcript profiles revealed a large number of changes in central carbon and nitrogen metabolism under oxygen-deprived conditions, especially during anaerobic germination. It is generally accepted that anoxia-tolerant plant species such as rice offset losses in aerobic ATP production by increasing fluxes through glycolytic and fermentative pathways, a phenomenon known as the Pasteur effect (Bailey-Serres and Voesenek, 2008). Here, we confirm previous reports that, in rice, this adjustment involves transcriptional up-regulation of glycolytic and fermentative enzymes with a shift toward ATP-conserving, PPidependent pathways of Suc catabolism (Bailey-Serres and Voesenek, 2008;Huang et al., 2008;Magneschi and Perata, 2009). A number of the metabolite responses we observed also support previous reports in oxygendeprived rice seedlings. These include increases in the amino acids Ala, Val, Gly, Leu, Arg, Tyr, Phe, Pro, and 4-aminobutyrate (Bertani et al., 1981;Reggiani et al., 1988), the organic acids succinate and lactate (Menegus et al., 1988(Menegus et al., , 1989, and the hexose phosphate Fru-6-P (Menegus et al., 1991) as well as decreases in the amino acids Glu, Gln, Asp, and Asn (Bertani et al., 1981;Reggiani et al., 1988), the organic acid 2-oxoglutarate (Reggiani et al., 1988), and the sugars Glc and Fru (Menegus et al., 1991). The polyamine putrescine has been previously shown to accumulate approximately 2.5-fold in 3-d-old aerobically grown rice seedlings transferred to anoxia for 2 d (Reggiani et al., 1989). Interestingly, we observed a similar approximately 2-fold increase in putrescine when 24-h aerobically grown rice seedlings were transferred to anaerobic conditions; conversely, we saw a transient 4-fold decrease in putrescine in anaerobically germinated rice seedlings compared with aerobically germinated seedlings at 24 HAI, indicating that preacclimation of seedlings in aerobic conditions is a prerequisite for anoxia-induced putrescine accumulation.
Our nontargeted gas chromatography-mass spectrometry (GC-MS) approach revealed a wide variety of oxygen-dependent metabolic responses that, to our knowledge, have not previously been reported in rice. In addition to the previously reported aerobic rice metabolite 2-oxoglutarate (Reggiani et al., 1988), the core aerobic-induced metabolite set included Ara, trehalose, citrate, shikimate, and gluconate, which have not been previously reported as oxygen-dependent metabolites in rice. Decreases in citrate and 2-oxoglutarate under oxygen-deprived conditions are consistent with a decrease of carbon flow into the TCA cycle due to a decrease in the NAD + /NADH ratio in the mitochondrial matrix, causing inhibition of the pyruvate dehydrogenase complex and other redox-sensitive TCA cycle enzymes. Citrate has been previously reported to decrease sharply in the oxygen-deprived roots of the relatively anoxia-sensitive species tomato (cv Bison; Dubinina, 1961), while both citrate and 2-oxoglutarate have been shown to decrease in Arabidopsis cell suspensions treated with the respiratory inhibitor rotenone (Garmier et al., 2008). However, decreases in citrate and 2-oxoglutarate do not appear to be general responses of plants to oxygen deprivation or respiratory perturbation, with citrate being quite stable and 2-oxoglutarate being increased in the roots of poplar during a 1-week period of flooding (Kreuzwieser et al., 2009) and citrate showing a moderate increase and 2-oxoglutarate showing a strong increase in the oxygen-deprived roots of pumpkin (Cucurbita pepo 'Mozoleevskaya'; Dubinina, 1961).
In our core anaerobiosis-induced metabolite set, only 4-aminobutyrate, succinate, and Phe have previously been reported to accumulate in anaerobic rice tissues (Bertani et al., 1981;Menegus et al., 1988Menegus et al., , 1989Reggiani et al., 1988), while Orn, coumarate, urate, glycolate, fumarate, nicotinate, glycerol-3-phosphate, and myoinositol have not. Glycolate accumulation has previously been reported in anaerobic pumpkin and tomato roots (Dubinina, 1961). In addition, urate and glycerol-3-phosphate have been reported to accumulate, fumarate and nicotinate were stable, and glycolate, Orn, coumarate, and myoinositol were not recorded in flooded poplar roots (Kreuzwieser et al., 2009). The accumulations of urate and glycolate in anaerobic tissues are readily explained by the oxygen requirement of their respective catabolic enzymes, urate oxidase (EC 1.7.3.3) and glycolate oxidase (EC 1.1.3.15), as suggested for anaerobic glycolate accumulation by Dubinina (1961). Similarly, major catabolic pathways of succinate and glycerol-3-phosphate are linked to oxygen availability via the mitochondrial respiratory chain-associated enzymes succinate dehydrogenase and glycerol-3-phosphate dehydrogenase (Shen et al., 2003(Shen et al., , 2006Rasmusson et al., 2008), respectively. In rice, coumarate may be formed from Tyr via Tyr ammonia lyase activity (Kishor, 1989), thus avoiding an oxygen-dependent step required for its synthesis from Phe. In lignin biosynthesis, coumarate is presumed to be metabolized to caffeate (3,4-dihydroxycinnamate) in an oxygen-dependent cytochrome P450-catalyzed reaction (Jaiswal et al., 2006;www. gramene.org). Hence, its accumulation may be explained by a block in this reaction. The accumulations of the remaining core anaerobic metabolites nicotinate, Phe, fumarate, Orn, and myoinositol are more difficult to explain, not having any obvious metabolic links with oxygen, and make interesting targets for further experimentation.
In addition to changes in transcripts and metabolites involved in carbon and nitrogen metabolism, 28 transcripts encoding proteins involved in lipid metabolism were present in the core anaerobic set. These included 13 transcripts specifically categorized as lipid metabolism and 15 transcripts encoding proteins associated with lipid metabolism (e.g. signaling/lipid metabolism). To our knowledge, a role for lipid metabolism in response to anaerobic conditions in plants has not been suggested before. These changes in lipid metabolism have the potential to play a central role in sensing and signaling anaerobic conditions. An absence of oxygen will result in an inability to synthesize polyunsaturated fatty acids due to the oxygen-requiring desaturase step, thus affecting membrane composition and fluidity. This has the potential to affect signaling, as has been previously described for stress signaling pathways in plants (Penfield, 2008). Second, changes in lipid metabolism will result in an inhibition of sterol synthesis, again affecting membrane fluidity and also signaling via sterol-binding transcription factors (Schrick et al., 2004). The enrichment of lipid metabolism function in the anaerobic set does correspond with observations in yeast, where sterol metabolism has been shown to be up-regulated and to play an important role in energy production and signaling under hypoxic conditions (Zitomer and Lowry, 1992;Kwast et al., 1998).
The changes observed in transcript abundance are ultimately mediated by transcription factors binding CAREs in the upstream regions of targeted genes. Twelve CAREs were found to be overrepresented in the promoter regions of anaerobic genes, whereas only two were found to be overrepresented in aerobic genes and four were underrepresented (Table I). Many CAREs enriched in the anaerobic set have been previously associated with hypoxia (refs. in Table I) and shown to be functional in individual genes. However, in addition to CAREs enriched in anaerobic genes, 15 CAREs were found to be enriched in both the anaerobic and aerobic sets. This may appear paradoxical at first; however, it has been proposed that a downregulation of genes that encode proteins that require oxygen to function may be an energy-saving mechanism under anaerobic conditions, but this requires repression of transcription (Bailey-Serres and Voesenek, 2008;Magneschi and Perata, 2009). Thus, CAREs present in both groups may act as positive regulators of anaerobic genes and negative regulators of aerobic genes under anaerobic conditions. Such a situation has been demonstrated for the Hap1 transcription factor in yeast, as it acts as a positive regulator of aerobic genes in the presence of oxygen (heme) and a negative regulator of ergosterol biosynthetic genes in the absence of oxygen (heme; Hickman and Winston, 2007). Another example of factors that can act as either positive or negative regulators are mammalian nuclear hormone receptors, where binding of the ligand switches from a repressor to an activator (Xu et al., 1999).
In terms of transcription factors that may bind the CAREs identified above, the bZIP family of transcription factors were overrepresented, both in germination under anaerobic conditions and in the anaerobic core set. bZIP transcription factors are associated with various stress responses in plants (Jakoby et al., 2002;Jain et al., 2008). Notably, some are membrane bound, in the endoplasmic reticulum, and are released by proteolysis, under various stimuli (Chen et al., 2008). In mammalian systems, a similar system occurs where transcription factors are released under conditions of low cholesterol (Hoppe et al., 2001), while in yeast, a similar system has been shown to be oxygen responsive (Martin et al., 2007). Thus, even though the 71 anaerobic transcription factors are candidate factors for the regulation of the response to anaerobic conditions, posttranslational mechanisms such as the activation of transcription factors may also occur, a process well described for bZIP transcription factors (Schutze et al., 2008). In addition, from the 71 transcription factors defined as core anaerobic genes, three have been annotated as having an insertion in the Tos17 Insertion Mutant Database library (Miyao et al., 2007). The phenotype for all three is noted as low fertility. Previously, it has been documented that anaerobic metabolism may occur during pollen development (Tadege et al., 1999). This is consistent with the role of these transcription factors as central regulators of anaerobic metabolism in rice. In terms of posttranscriptional regulation, it was observed that the transcripts of genes defined as aerobic contained five motifs enriched in the 3# UTR, one previously defined as a destabilization element (Narsai et al., 2007). Also, aerobic genes were significantly underrepresented in one element defined as a stabilization element (Narsai et al., 2007). In contrast, the anaerobic genes analyzed contained no enriched elements in the 3# UTR. One explanation for this could be that in rice, the transcripts defined as aerobic transcripts may be actively degraded under anaerobic conditions, a process that has been previously described in yeast (Dagsgaard et al., 2001). Notably, posttranscriptional regulation of ribosome loading has also been considered to be an important regulatory mechanism under hypoxic conditions (Branco-Price et al., 2005, indicating that under anaerobic conditions specific regulation can occur at the posttranscriptional and posttranslational levels.

Rice Growth
Dehulled, sterilized rice seeds (Oryza sativa 'Amaroo') were grown under aerobic or anaerobic conditions in the dark at 30°C as described previously . Embryos were rapidly dissected from the endosperm and snap frozen in liquid nitrogen.

RNA Isolation
Total RNA was isolated from rice embryos as described previously (Howell et al., 2006). Three independent RNA preparations were performed for each developmental stage/growth condition, and the concentration of RNA was determined spectrophotometrically.

Microarray Analyses
Affymetrix GeneChip Rice Genome Arrays were used for the transcriptomic analysis as described previously ). Data quality was assessed using GCOS 1.4 (Affymetrix) before CEL files were imported into Avadis 4.3 (Strand Genomics) for further analysis. Raw intensity data were initially normalized using the MAS5 algorithm allowing probe set identifiers that were called present to be determined. Only those probe sets that were called present in at least two of three replicates in at least one time point were included for further analysis. Ambiguous probe sets and bacterial controls were also removed, resulting in a final data set of 29,087 probe sets across the germination and switch array experiments.
Following filtration of the data set, all probe intensities were then analyzed using the GC (content) robust multiarray average (GC-RMA) algorithm, and three-dimensional principal component analysis was carried out to visualize the global differences between the arrays. To perform differential expression analysis, the GC-RMA normalized data were then log transformed and false discovery rate correction (Benjamini and Hochberg, 1995) was applied at the level of P , 0.01. This allowed the number of transcripts significantly changing to be calculated, which were then visualized on a heat map. The differential expression analysis was carried out using Avadis 4.3 (Strand Genomics), while the heat maps, PCAs, and hierarchical clustering were all carried out using Partek Genomics suite software, version 6.3. MapMan (Thimm et al., 2004;Usadel et al., 2005) analyses were performed using a reduced set of unique probe sets (17,722). The mapping file used was generated as described by Howell et al. (2009). All GC-RMA normalized array data and present/absent calls are given in Supplemental Table S1.

Generation of Transcription Factors and Organelle Lists
The transcription factor list was generated using three main sources, DRTF (Gao et al., 2006), RiceTFDB (Riano-Pachon et al., 2007), and Caldana et al. (2007). These lists were compiled, and all of the 3,098 unique transcription factors were matched to the 29,087 genes to generate a list of 2,009 transcription factors. To examine the transcripts in the aerobic and anaerobic responsive sets that encoded mitochondrial, chloroplast, and peroxisomal proteins, we matched these identifiers to the localizations assigned by Howell et al. (2009), and any remaining identifiers were matched exactly as described previously ).

Functional Annotation and Statistical Analysis
These were performed exactly as described by Howell et al. (2009).

Public Rice Microarray Data Analysis and Comparison
In order to examine transcript abundance changes across different tissues, under different conditions, and compare these with the transcript abundance profiles generated from this study, rice array data were retrieved from the Gene Expression Omnibus within the National Center for Biotechnology Information database as described previously ) with the addition of data derived from cytokinin treatment of roots and leaves (GSE6719; Hirose et al., 2007) and infection with Striga hermonthica (GSE10373; Swarbrick et al., 2008).

Promoter Motif Analysis
Following expression analysis, the 12 groups of transcripts selected for promoter analysis (Supplemental Table S7A) were analyzed for the presence of overrepresented sequence elements in the 1-kb upstream regions. In order to study these coexpressed transcripts more closely, all 1-kb upstream regions of the 29,087 transcripts were retrieved and the upstream regions of the 12 subgroups were extracted and examined for putative cis-acting elements. Programs designed to detect sequence elements generally have limits of approximately 50 to 60 input sequences; thus, these smaller lists were selected. The MEME Web server (Bailey et al., 2006) was used under default settings with the length of the motif set to 6 to 8 bp and the number of motifs to find set to five (instead of the default of three). The output from this analysis is shown in Supplemental Table S7B. The lists of motifs from Supplemental Table S7B were then filtered only to include motifs present in more than 70% of all input sequences (Supplemental Table S7C), and the presence of these motifs was then examined in the whole genome and the 12 subsets.

3# UTR Sequence Analysis
The full genome 3# UTR and 5# UTR sequences available from The Institute for Genomic Research were downloaded and filtered to retain only the 3# UTRs. However, this only included a total of 3,027 UTRs available for the "whole genome." Taking this small number into consideration, it was not feasible to look at the 12 subsets, as these lists were too small. Thus, for the 3# UTR, all of the genes in the aerobic and anaerobic responsive subsets were used for the analysis of overrepresented sequence elements in the 3# UTRs. The settings were set to search for five motifs that are 6 to 8 bp (default) in each of the subsets, and the outputs are shown at the bottom of Supplemental Table  S7B. It is important to note that setting the output to five motifs can result in false present calls for motifs that are not significant when the input list is small; therefore, only the significantly enriched motifs (present in 60%-70% of all input sequences) were included for further analysis (Supplemental Table  S7D). In addition to these putative predicted motifs, other motifs known to be associated with RNA stability/instability Ohme-Takagi et al., 1993;Narsai et al., 2007) were examined for their presence in these data sets compared with the whole genome (Supplemental Table S7C).

Metabolomic Analysis
Metabolite extraction, derivatization, and GC-MS analysis were all carried out as outlined by Howell et al. (2009). Briefly, metabolites were extracted with a hot aqueous methanol solution containing ribitol as an internal standard. Aliquots of extracts were dried under vacuum and derivatized by methoximation and trimethylsilylation. An n-alkane internal retention index standard mixture was added to each sample prior to GC-MS analysis on an Agilent GC/MSD system. Raw data files were automatically processed and matched to an in-house mass spectral and retention index library using in-house MetabolomeExpress software (A.H. Millar and A. Carroll, unpublished data), which generated a [metabolite 3 sample] data matrix without missing values that contained metabolite peak area values that had been normalized to dry tissue mass and internal standard (ribitol) peak area. Automatic statistical analysis of processed data was carried out by calculating, for each set of biological replicates, the mean normalized signal intensity for each metabolite and then, for each metabolite, dividing the mean signal in treated sample sets by the mean signal in control sample sets to calculate fold difference and then testing the statistical significance (P , 0.05) of this difference by Welch's t test. Heat map statistical results tables were exported from MetabolomeExpress and manipulated in Microsoft Excel 2008. Principal component analysis of metabolite profiles was carried with Partek Genomics suite software, version 6.3, using the above [metabolite 3 sample] data matrix after filtering with MetabolomeExpress to remove redundant metabolite signals, internal standards, and known analytical artifacts.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Principle component analysis of changes in transcripts and metabolites under the three experimental conditions used in this study.
Supplemental Figure S2. Hierarchical clustering of the changes in transcript abundance.
Supplemental Figure S3. MapMan overview of changes in transcript abundance for genes encoding components involved in metabolism.
Supplemental Figure S4. Changes in transcripts encoding components of lipid metabolism and lipid metabolism-related functions in the core anaerobic set.
Supplemental Figure S5. Hierarchical clustering analysis of the aerobic and anaerobic transcription factor genes across publicly available rice arrays.
Supplemental Table S1. All 29,087 expressed genes, GC-RMA normalized values, and these values normalized to maximum.
Supplemental Table S2. The number of differentially expressed genes during germination and following the switch conditions.
Supplemental Table S4. Lists of all detected metabolites under the germination and switch conditions with averaged raw metabolite abundance data and SE.
Supplemental Table S5. All 166 metabolites that were found to significantly change in abundance under the germination and/or switch conditions Supplemental Table S6. The 28 identified START domain-containing genes (Schrick et al., 2004).
Supplemental Table S7. Promoter analysis Supplemental Table S8. Z-score analysis of all transcription factor families, including those shown in Figure 5.
Supplemental Table S9. Public data were compared with the rice oxygenand anaerobiosis-responsive transcription factor sets.