|Home | About | Journals | Submit | Contact Us | Français|
The transcriptional regulatory network involved in low temperature response leading to acclimation has been established in Arabidopsis. In japonica rice, which can only withstand transient exposure to milder cold stress (10°C), an oxidative-mediated network has been proposed to play a key role in configuring early responses and short-term defenses. The components, hierarchical organization and physiological consequences of this network were further dissected by a systems-level approach.
Regulatory clusters responding directly to oxidative signals were prominent during the initial 6 to 12 hours at 10°C. Early events mirrored a typical oxidative response based on striking similarities of the transcriptome to disease, elicitor and wounding induced processes. Targets of oxidative-mediated mechanisms are likely regulated by several classes of bZIP factors acting on as1/ocs/TGA-like element enriched clusters, ERF factors acting on GCC-box/JAre-like element enriched clusters and R2R3-MYB factors acting on MYB2-like element enriched clusters.
Temporal induction of several H2O2-induced bZIP, ERF and MYB genes coincided with the transient H2O2 spikes within the initial 6 to 12 hours. Oxidative-independent responses involve DREB/CBF, RAP2 and RAV1 factors acting on DRE/CRT/rav1-like enriched clusters and bZIP factors acting on ABRE-like enriched clusters. Oxidative-mediated clusters were activated earlier than ABA-mediated clusters.
Genome-wide, physiological and whole-plant level analyses established a holistic view of chilling stress response mechanism of japonica rice. Early response regulatory network triggered by oxidative signals is critical for prolonged survival under sub-optimal temperature. Integration of stress and developmental responses leads to modulated growth and vigor maintenance contributing to a delay of plastic injuries.
Plants in temperate environments are able to withstand sub-freezing conditions by cold acclimation, whereas those that are adapted to tropical or sub-tropical environments can only maximally endure milder chilling temperatures . Arabidopsis has been a particularly well scrutinized model system for low temperature stress response mechanisms because of its ability to cold acclimate [2-4]. Cold acclimation is a manifestation of the coordinated function of a large number of genes under the control of few transcriptional regulators that respond directly to low temperature signals through abscisic acid (ABA) independent signal transduction [5,6].
The established resources for functional genomics in Arabidopsis facilitated a systematic approach for assembling the components of regulatory networks involved in cold acclimation as well as those involved in similar types of abiotic stresses . The most well documented mechanism involves a class of ERF transcription factors (TF) known as the DREB/CBF, which interact with the DRE/CRT cis-elements in the promoters of their downstream target genes to execute a highly coordinated transcriptional response to low temperature signals [8-10]. The DREB/CBF regulon is a large cluster consisting of more than a hundred genes with associated sub-regulons controlled by RAP2.1 and RAP2.6 [8,9,11]. Acting upstream is the MYC-type DREB/CBF activator ICE1, which is regulated by phosphorylation . Discovery of other regulons that function during cold acclimation is a continuing story. Recently, a DREB/CBF-independent regulon controlled by a zinc finger TF (ZAT12) was discovered [3,8]. ZAT12 acts in parallel to the CBF/DREB pathway although each mechanism appears to activate overlapping sets of COR genes. Additional regulatory clusters controlled by RAV1 TF have also been inferred based on transcriptome profiles of overexpression lines [3,9].
Rice is another important model for dissecting stress response regulatory networks because of its 'finished' genome sequence and rapidly developing functional genomics resources. Unlike the freeze-tolerant Arabidopsis, most rice cultivars are irreversibly injured by prolonged exposure to temperatures above 10°C (chilling), particularly at early seedling stage. However, many japonica cultivars are moderately adapted to temperate climates and they generally exhibit lesser degrees of chilling sensitivity than most tropical indica cultivars. Previous efforts to profile the low temperature transcriptome of japonica rice was based on the responses to 4°C , which is below the typical seedling stage LT50 of 13°C for even the most chilling tolerant japonica cultivars [14-16]. The scope of transcriptional changes that occur when japonica rice is exposed to chilling (10°C) has been investigated only recently in an effort to establish a more agronomically meaningful picture of response mechanisms [14,17,18]. In our earlier analysis of about 6,000 genes, we recognized peculiar trends suggesting that transient oxidative stress during the initial 24 hours triggered an 'early response' regulatory network that seems to play an important role in short-term adaptive responses . To dissect the various components of this network and other aspects that are independent of oxidative signals, we performed two parallel genome-wide expression profiling experiments comparing the responses of japonica rice to chilling (10°C) and cold-independent mimic of oxidative stress.
With the systems level approach employed in this study, we established a broader picture of the hierarchical organization and compositional complexity of upregulated transcriptional network supported by direct physiological evidence that oxidative signal caused by transient increases in intracellular H2O2 is a primary trigger for network activation. Potential agronomic significance of oxidative-mediated network was established by integrating different layers of information derived from genome-wide, physiological and whole-plant level analyses. The results of this study indicate that expression of early response regulatory network is a critical determinant of the ability of japonica rice for prolonged survival under sub-optimal temperature conditions. This is achieved by integrating defense and growth related responses leading to the ability to delay the progression of irreversible stress injuries by maintenance or enhancement of vigor.
Two parallel transcript profiling experiments were performed using the NSF-45K oligonucleotide microarray . In the first experiment (GSE8767), gene expression in cv. Nipponbare was monitored during a 96-hour duration at 10°C covering both early (0.5, 2, 4, 6, 12, 18, 24 hours) and late (36, 48, 96 hours) responses. In the second experiment (GSE10062), gene expression in response to applied H2O2 (4 mM) after 1, 3 and 6 hours of treatment at optimum temperature (28°C) was investigated also in cv. Nipponbare. The overall assumption was that gene activation by oxidative signals and those triggered by other physiological perturbations can be inferred by subtracting the composition of one dataset from the other.
Chilling stress has massive effects on the transcriptome, which is comprised of 8,668 genes that exhibited temporal changes during a 96-hour duration at 10°C. The upregulated components were identified by querying the whole microarray dataset at an arbitrary threshold of log2 ≥ 1.8 (p < 0.05), revealing a total of 2,604 unique genes that were induced within a duration of at least two consecutive time points. Nearly 60% (1,516) of these genes were also induced by 4 mM H2O2 (p < 0.05) in at least one time point (Figure (Figure1).1). This overlap was assumed to represent the components of chilling stress response transcriptome that were activated by an oxidative-mediated regulatory mechanism.
The downregulated dataset does not overlap with the upregulated dataset. It includes all the genes (total = 6,064) with expression values of log2 ≤ -1.8 (p < 0.05) in at least one time point and without any value equal to or greater than the upregulated threshold (log2 1.8) in any one time point. Of these, 25% (1,516) were also downregulated by 4 mM H2O2 (Figure (Figure1).1). These genes were assumed to represent the components of chilling stress transcriptome that were negatively regulated by an oxidative-mediated mechanism. The composition and temporal profiles of the downregulated group were quite complex. A summary of the distribution of these genes according to broad functional categories is presented (see Additional file 1).
Downregulation appears to have major impact on cellular processes related to growth, development and morphogenesis, cell division, metabolism and transport mechanisms. Due to the complexity of the temporal patterns exhibited by this group of genes, the assembly of underlying regulatory network and their linkages with the upregulated network will be described in a separate report based on the analysis of transcription factor overexpression lines.
Transcription factors (TF) link stress signals to the downstream genes that execute cellular defense processes . Of the 2,604 chilling upregulated genes, about 6% (148) are TFs (Figure (Figure1,1, Table Table1),1), which were identified based on the classification of the Database of Rice Transcription Factors [21,22]. These genes represent more than 7% of the total TFs encoded by the rice genome and include members of AP2/ERF, bZIP, MYB, WRKY, bHLH and NAC families [3,6,23-26]. TFs exhibited complex expression patterns characterized by waves of induction initiating at different time periods (Figures (Figures2,2, ,3,3, ,4;4; see Additional files 2, 3, 4). For instance, activation of the 'early rapid response' group (phase-1) occurred during the initial 6 hours, while those that belong to the 'early slow response' group (phase-2) occurred between 6 and 24 hours. Few TFs exhibited 'late response' profiles (phase-3) with no significant induction until after 24 hours. Majority of chilling upregulated TFs were activated at phase-1 and phase-2, indicating that the critical transcriptional changes occurred during the initial 24 hours.
Induction waves were either short-lived or long-lived events lasting between 2 to 60 hours. Some genes are activated multiple times based on interrupted patterns (Figures (Figures2,2, ,3,3, ,4).4). Multi-phasic profiles reflect the temporal hierarchy of regulatory networks operating at 10°C and the complexity of signals that trigger the activity of regulatory clusters. Apparently, phase-1 TFs are closest to the original signal and are likely to function as primary switches that trigger subsequent cascade of gene activation. Reactivation of some phase-1 TFs at later time points suggests that multiple signals are acting on the same TF. Phase-2 and phase-3 TFs are likely controlled by other phase-1 TFs that respond directly to the primary signal. Similar hierarchical organization of stress related TFs has been shown earlier. For example, DREB/CBF controls two other types of ERF genes (RAP2.1, RAP2.6) involved in the activation of DRE/CRT element-containing genes as well as an R2R3-type MYB that control other potential sub-regulons related to CBF/DREB network [3,9,27].
Members of the bZIP family of TFs have important roles in ABA response and regulation of oxidative and pathogen defense responses [5,28]. A total of 11 genes from this family were upregulated by chilling (Figure (Figure2),2), representing five sub-classes based on recent phylogenetic classification [24,29]. The 'early response' (phase-1 and phase-2) genes belong to Group-D associated with oxidative stress and defenses against pathogens (Os06g41100: TGA10), Group-H associated with photomorphogenesis (Os02g10860: HY5, Os12g06520: PosF21), Group-I associated with GA response (Os08g43090: RF2b-group, Os04g10260: bZIP domain protein), and Group-S associated with stress response and sucrose signaling (Os08g38020: AtbZIP148-like) and a novel gene that appears to be specific to rice (Os0410260). The ABA-associated bZIP genes belonging to Group-A (Os02g52780: ABI5, Os01g64000: AREB) were not induced until after 24 hours (phase-3). Several Group-D (Os05g37170: TGA6) and Group-S (Os02g09830, Os02g03960: ocs-binding proteins) genes were also induced at phase-3. Majority of chilling induced bZIP genes were also induced by H2O2 (Figure (Figure2).2). Many (Os06g41100, Os05g37170, Os02g09830, Os02g03960) are known to interact with as1/ocs/TGA-like elements involved in H2O2, auxin, disease, salicylic acid and jasmonic acid regulated gene expression [30-33].
ABA induces H2O2 synthesis in guard cells via NADPH oxidase [34,35]. This was indicated by impaired stomatal closure in double mutants for NADPH oxidase AtrbohD and AtrbohF genes, which was reversed by exogenous H2O2 . ABA-related bZIP genes of rice (Os02g52780, Os01g64000) were also induced by H2O2 but their late induction implies direct consequences of ABA rather than oxidative signaling. We observed that rice seedlings had visible symptoms of leaf rolling not earlier than 24 hours after exposure to 10°C (data not shown), which was coincident with the temporal expression of ABA-related bZIP genes at phase-3. This pattern indicates that oxidative-mediated gene expression occurred much earlier (within the initial 6 hours) than those mediated by ABA (after 24 hours), and that ROS (H2O2) was the primary signal that integrates early responses (Figure (Figure5).5). This is consistent with our earlier observations that many chilling induced genes that were also induced by ABA were 'late response' genes [, unpublished data].
About 18% (31) of AP2/ERF genes of rice were induced by chilling (Table (Table1;1; Figure Figure3).3). Most are ERF-type belonging to nine phylogenetic groups [3,5,20,37]. Genes involved in indole alkaloid metabolism (ORCA-type) and responses to diseases, wounding and elicitors such as ethylene, salicylic acid and jasmonic acid comprise one of the larger classes which include Group-VI/subfamily B-5 (Os06g06540), Group-VIII/subfamily B-1 (Os04g52090, Os06g47590, Os05g41780, Os06g47590, Os01g58420, Os04g57340), and Group-IX/subfamily B-3 (Os02g43790, Os07g22730, Os08g44960, Os09g39810). Many are known to interact with GCC-box and jasmonic acid response element (JAre) [38-41]. Genes in this category were expressed in all three phases but only few were responsive to H2O2 (Figure (Figure33).
Another large category of chilling-induced ERF represent a class of genes homologous to the low temperature, dehydration and salinity stress regulators in Arabidopsis belonging to Group-II/subfamily A-5 (RAP2.1/2.9/2.10 group: Os03g15660, Os06g11940), Group-III/subfamily A-5 (DREB/CBF group: Os08g43210, Os02g45450, Os01g73770, Os09g35010, Os06g03670, Os09g35020) and Group-VII/subfamily B-2 (Os03g08500, Os03g08490, Os12g40960). Members of these groups interact with DRE/CRT and related cis-elements [6,42,43]. Only few of these genes were induced by H2O2, suggesting that this class of ERF is not a direct target of oxidative signaling although cross-talk is possible as indicated by few DREB/CBF genes that responded to the H2O2 treatment (Figure (Figure3).3). Specific DREB/CBF genes from Group-III were activated at phase-1 (Os02g45450: DREB1A/CBF3, Os08g43210/Os09g35010: DREB1B/CBF1, Os01g73770: DREB1C/CBF2), phase-2 (Os06g03670: DREB1A/CBF3) and phase-3 (Os09g35020: DREB1D/CBF4), indicating that different members respond to specific signals at various stages of stress.
The third category of chilling induced ERF genes are associated with growth regulation and ABA signaling from Group-I/subfamily A-6 (Os04g44670, Os06g11860, Os03g09170), Group-IV/subfamily A-2 (Os05g28350) and Group-X/subfamily B-4 (Os09g28440) [44,45]. Genes in this category were activated either at phase-1 or phase-2 and many were also induced by H2O2 (Figure (Figure3).3). Early induction appears to be a consequence of oxidative signaling rather than ABA, given that ABA response does not seem to occur until after 24 hours after chilling exposure.
Several non-ERF-type TFs were also induced by chilling, including the RAV1-type ABI3/VP1 homologs (Os01g04750, Os01g04800), genes encoding baby boom-like proteins (Os09g25600, Os02g40070) and a novel ERF3-like gene (Os08g41030) (Figure (Figure3).3). RAV1-type TFs function in low temperature response and circadian regulation [5,9], while the baby boom-like proteins are involved in development . Rice seedlings at three-leaf (V3) stage, which is about 12 days after germination, were used in the microarray experiments. Activation of RAV1-type and baby boom-type transcription factors reflects the cross-talk that typically occurs between low temperature response and developmental processes and this is probably important to integrate signals necessary for normal seedling growth and immediate defenses to stress.
Homologs of ERF genes that function maximally during cold acclimation at 4°C in Arabidopsis were activated at milder (10°C) temperature . The regulatory clusters and physiological consequences are likely to be different under non-acclimating (10°C) and acclimating (4°C) temperature regimes based on the observation that CBF/DREB induction at 10°C tend to be relatively short-lived and less robust . Nevertheless, the transcriptional network at 10°C provides an agronomically meaningful view of response mechanisms since the LT50 for most japonica cultivars is about 13°C [14,16,46]. ERFs involved in disease response via oxidative, jasmonic acid and salicylic acid signaling are also involved in chilling stress response. However, they appear to act on a different set of target genes from those regulated by DREB/CBF.
MYB domain proteins represent the largest class of TFs involved in cell cycle regulation, cell fate determination and responses to environmental stresses [47-49]. MYB is also the largest group of chilling induced TFs with 37 members (Figure (Figure4)4) belonging to three sub-families . The two largest are the R2R3 subfamily (Os02g41510, Os01g50110, Os01g18240, Os09g23620, Os05g04820, Os01g45090, Os09g36730, Os09g26170, Os10g33810, Os10g35660, Os04g43680, Os11g45740, Os12g37690, Os09g36730, Os05g46610) and R1-MYB/MYB-like subfamily (Os01g09760, Os10g30719, Os01g09280, Os05g07010, Os06g07640, Os06g19980, Os05g10690, Os02g56030, Os02g10060, Os04g40420, Os04g41830, Os08g04840, Os01g03660, Os02g45670). Many R2R3 genes are known to function in cold, salinity and ABA response mechanisms [51-53]. A total of 24 MYB genes were activated at phase-1, 5 at phase-2 and 8 at phase-3 (Figure (Figure4).4). Phase-1 MYB genes tend to be more responsive to exogenous H2O2 than those in phase-2 and phase-3, indicating that 'early response' genes are possible direct targets of oxidative signals including the OsMYB4 (Os02g41510) that we previously reported .
WRKY TFs function as positive or negative regulators of defenses against pathogens and herbivores via salicylic acid and jasmonic acid signaling pathways [54-58]. Their precise role in abiotic stress response regulatory network is not fully understood [59,60]. A total of 17 members of this family were induced by chilling (see Additional file 2). Most are 'early response' genes, 10 of which were induced at phase-1 and 5 at phase-2. Phase-1 genes include several known regulators of salicylic acid and jasmonic acid mediated disease response (Os01g14440: OsWRKY1v2, Os09g24070: OsWRKY62, Os03g55164: OsWRKY4, Os05g39720: OsWRKY53, Os11g45850: OsWRKY40, Os05g04640: OsWRKY5), and few have recently been shown to function in ABA-mediated gene expression (Os01g61080: OsWRK24, Os11g29870: OsWRKY72). OsWRKY24 is particularly interesting because it was recently shown to act as repressor of ABA-mediated induction of HVA22 [60,61]. Upregulation of this gene at phase-1 is quite consistent with our current data showing that ABA-regulated expression of ABI5 and AREB genes occurred largely at phase-3. OsWRKY24 appears to act as repressor of ABA-inducible genes during the early stages of chilling stress and it is likely involved in the integration of responses to oxidative and ABA signals.
Phase-2 WRKY genes are mostly disease-related (Os05g27730; OsWRKY53; Os05g50610: OsWRKY8v2, Os01g51690: OsWRKY26) and herbivory-related (Os03g58420: OsWRKY6). Downregulation of OsWRKY53 by ABA was reported recently . OsWRKY53 expression lasted only up to phase-2 which is probably due to increased ABA levels at phase-3. Three of the phase-2 WRKY genes were induced by H2O2 and they are probably secondary targets of oxidative signaling through other TFs that are induced at phase-1. Their activation by H2O2 may also be mediated by salicylic acid or jasmonic acid . It is possible that some of the phase-1 and phase-2 WRKY genes are involved in repression of gene expression to fine-tune the interaction of early (oxidative) and late (ABA) signaling mechanisms. Enrichment of W-box-like motifs among downregulated genes are consistent with this hypothesis (data not shown).
TFs from the bHLH family are involved in plant development, morphogenesis, circadian regulation and stress response . MYC-type factors are involved in ABA-regulated gene expression [12,26,37]. About 19% (29) of rice bHLH genes were induced by chilling, 18 of which were at phase-1, 5 at phase-2 and 6 at phase-3. Half of 'early response' bHLH genes were responsive to exogenous H2O2 and their profiles were characterized by multiple waves reflecting possible roles in combinatorial control mechanisms (see Additional file 3).
NAC (NAM, ATAF1/2, CUC2) genes encode a group of plant-specific TFs involved in growth and hormone signaling, apical meristem and floral development, senescence, disease response via salicylic acid and jasmonic acid and abiotic stress responses via ABA [23,64,65]. A total of 23 NAC genes were induced by chilling, including 21 'early response' (14 at phase-1 and 7 at phase-2) and 2 'late response' (Table (Table1;1; see Additional file 4). Only few of the NAC genes (mostly phase-1 genes) were also induced by H2O2.
Cellular concentrations of ROS including H2O2 are tightly regulated by balancing synthesis and degradation. Exposure to sub-optimal conditions often disturbs such balance [66,67]. At a critical level within specific cellular microdomains, H2O2 acts as secondary messenger that integrates a variety of responses associated with biotic and abiotic signals and growth and development [68-71]. We have shown in our previous report that exogenous H2O2 mimics to a certain degree the effect of chilling on gene expression by activating many chilling-associated ROS scavengers and other stress related genes when plants were treated with 4 mM H2O2 at ambient temperature . We hypothesized that a large component of chilling stress transcriptome is a consequence of ROS (H2O2) acting as primary signals for the activation of 'early response' regulatory networks. To address this hypothesis further, we profiled the temporal fluctuation of intracellular H2O2 in Nipponbare at 10°C and 28°C (Figure (Figure55).
Based on the changes in [chilling]/[control] ratios across time, H2O2 levels were particularly high during the initial 12 hours. H2O2 profile was characterized by transient waves indicating the interplay of stress-induced synthesis and degradation. H2O2 spikes differ in peak heights and occurred in three phases, first during the initial 2 hours followed by two more spikes between 4 and 6 hours and between 8 and 9 hours. After about 12 hours relative concentrations in stressed and control plants were stabilized as a consequence of the activation of ROS scavenging systems. Stabilization of H2O2 spikes was coincident with high expression of many ROS scavenging genes after 12 hours (see Additional file 4).
Of the 148 chilling upregulated TFs, 62 were also upregulated by exogenous H2O2 (Table (Table1).1). Activities at phase-1 and phase-2 mirror the temporal profile of cellular H2O2 accumulation. For instance, activation of many 'early response' TFs that were also responsive to H2O2 coincided with the occurrence of H2O2 spikes within the initial 12 hours. These temporal coincidences provide an indirect but strong support that H2O2 is involved as a primary signal for the activation of 'early response' regulatory network. H2O2 and other ROS are known to cause conformational changes of the DNA-binding domain of certain TFs [72-74]. We previously proposed that a constitutively expressed but inactive TF (e.g., MYB) could be activated by changes in redox status at 10°C and this possibly leads to rapid induction of its target TFs (e.g., TGA10) occurring at phase-1 . Our current data support this hypothesis and provide a basis for possible concentration dependence of H2O2-inducible gene expression .
Genes encoding non-transcription factor (NTF) proteins (2,456) were further classified based on expression profiles to establish the composition of various co-expressed groups (regulatory clusters) controlled by the different classes of TFs listed in Table Table1.1. Genes were assigned to smaller groups by K-mean clustering to determine the closest similarities in expression . A total of 18 clusters were established with members tightly distributed around the group mean (Figures (Figures6,6, ,7;7; see Additional file 5). Based on consensus expression profiles, genes can be assigned to a three-phase induction pattern similar to that observed among the various TF classes. Clusters C78, C79, C83A, C86, C86A, C88A, C90A, C91, C94, C97, C97A and C99A comprised the 'early rapid response' group showing activation at phase-1 (Figure (Figure6).6). Clusters C80, C81, C83, C87A and C100A comprised the 'early slow response' group showing activation at phase-2. Only C99 belongs to the 'late response' group (Figure (Figure7).7). Genes that were activated at phase-1, phase-2 and phase-3 accounted for 71%, 22% and 7% of the total NTF genes in the dataset, respectively.
Genes that were also induced by H2O2 were more predominant in phase-1 clusters than in phase-2 and phase-3 clusters. Of the twelve phase-1 clusters, only C78 (39%) and C88A (16%) had low proportion of H2O2 induced genes. Between 50 to 79% of the genes in ten other phase-1 clusters (C79, C83A, C86, C86A, C90A, C91, C94, C97, C97A, C99A) were induced by H2O2 (Figure (Figure6).6). In contrast, the proportion of H2O2-induced genes among phase-2 and phase-3 clusters ranged only between 12 to 49% per cluster (Figure (Figure7).7). Apparently, genes that were rapidly induced during exposure to chilling (initial 6 hours) are more likely to respond to exogenous H2O2 than those that were induced at a slower pace (see Additional file 6). These results further support the hypothesis that oxidative burst during the early stages of chilling stress was a primary signal for early gene induction events, and that 'early response' transcriptome is controlled largely by TFs responding directly to primary oxidative signals. This trend is consistent with the induction of large number of genes associated with disease response, wounding, and ROS scavenging at phase-1, many of which are associated with various signals (ethylene, ABA, auxin, salicylic acid, jasmonic acid) linked directly or indirectly to H2O2-mediated processes [38,62,69,75].
Promoter elements shared in common by the majority of genes within a cluster were identified to gain further insights on how gene induction is coordinated . Promoter regions (-1,000 to +200) of more than 70% of genes in the total dataset were determined by aligning full-length cDNA with corresponding genomic loci. Delineated promoter regions were used for ab initio detection of putative cis-elements by the Dragon Motif Builder algorithms. Motifs with occurrence of 50% or higher in a given cluster were matched with known elements in the TRANSFAC, PLACE and AGRIS databases. Motifs were assigned to specific classes based on known association to specific classes of TFs (Tables (Tables2,2, ,3,3, ,4).4). The most significantly enriched motifs are related to several classes of elements associated with MYB, bZIP and ERF factors. Motifs associated with TF family were enriched at various levels in 18, 13 and 10 clusters, respectively (Figures (Figures6,6, ,7;7; Tables Tables2,2, ,3,3, ,4).4). WRKY associated motifs were significantly enriched only in 2 clusters. However, some of the possible WRKY target motifs (W-box) were quite similar with the as1/ocs/TGA-like motifs associated with bZIP factors and this may have reduced the occurrence of WRKY-associated cis-elements in the total gene set [31,75]. Motifs associated with bHLH and NAC had low occurrence at the 50% threshold. Other classes of motifs were also detected but could not be assigned to specific TFs (see Additional file 7).
The predominance of MYB, bZIP and ERF associated motifs were also quite apparent when the enrichment analysis was performed on the entire dataset (2,456 genes) as a single group (without clustering) and compared to a subset of known genes not involved in stress response (background control). This result indicated that enrichment was not due to random occurrence of frequently occurring sequences in the japonica rice genome. To assess the biological significance of ab initio prediction, a select group of putative cis-elements from the major clusters were tested for binding with nuclear proteins in vitro by electrophoretic mobility shift assay (EMSA). The rav1-like (TCT(a/c)AACA), as1/ocs/TGA-like (AATTTGAT, TAATTTGA), pyrimidine box-like (AAAGAAAAA) and MYB2-like (TAGTTTTT) motifs enriched in C81, C83A, C79, C97A and C94, respectively, showed band shifts in the presence of pooled nuclear proteins from chilling stressed (10°C) seedlings (Figure (Figure8).8). These results indicate that the sequences are bona fide binding sites of low temperature-induced nuclear proteins.
Motifs of the as1/ocs/TGA-like and ABRE-like classes were the most significantly enriched among the possible bZIP-target sequences [13,75,76]. Occurrence of the as1/ocs/TGA-like class was particularly prominent (score = 100 to 350) in C79, C83, C83A, C86, C86A, C90A, C91, C94, C97, C99A and C100A (Figures (Figures6,6, ,7;7; Table Table2).2). Only two of these clusters (C83, C100A) were activated at phase-2 and the rest were activated at phase-1, indicating that as1/ocs/TGA-like motifs are critical features of rapid response genes. Clusters enriched with as1/ocs/TGA-like elements also appeared to be more responsive to exogenous H2O2 based on high proportion (53 to 79%) of H2O2-induced genes among those clusters. This is an indication of direct relationships between rapid induction, as1/ocs/TGA-like element and responsiveness to H2O2, which we also observed in our previous analysis of a smaller subset of upregulated genes . The as1/ocs/TGA elements are key components of regulatory modules involved in auxin, salicylic acid and jasmonic acid mediated gene expression in response to wounding, oxidative stress and pathogen attack, providing direct support that 'early response' genes are direct consequences of primary oxidative signals [28,30,31,33,68,77]. This is also in agreement with the timing of H2O2 spikes within the initial 6 to 12 hours at 10°C (Figure (Figure55).
Activation of as1/ocs/TGA-like element-enriched clusters can be associated with several phase-1 (Os08g43090: RF2b-like protein, Os06g41100: TGA10, Os08g38020: AtbZIP148-like protein) and a phase-2 (Os04g10260: bZIP domain protein) bZIP genes that were also induced by H2O2 (Figure (Figure2).2). In tobacco, TGA10 has been shown to directly bind to an as1/ocs/TGA-like sequence and its transient expression enhanced the expression of downstream pathogenesis-related (PR) target genes via auxin, salicylic acid and jasmonic acid signaling pathways . Consistent with this, our recent results showed that as1/ocs/TGA-like element-containing genes were among those that were activated in a chilling or H2O2-independent manner in transgenic rice overexpressing TGA10 (data not shown). In addition, three other H2O2-induced bZIP genes that were activated at phase-3 (Os02g09830/Os02g03960: ocs-binding factor, Os05g37170: TGA6) are also likely to be involved in the regulation of as1/ocs/TGA-like element-enriched clusters. These genes could be the direct targets of bZIP genes induced at phase-1, thus they are possible components of ROS-bZIP sub-regulons. The ROS-bZIP genes belong to the groups involved in defenses to pathogens, photomorphogenesis, GA response and sucrose signaling, and this suggests important roles of their downstream regulon in the integration of stress and developmental responses [24,29].
ABRE-like motifs typically found among ABA-regulated genes were particularly enriched (scores = 100 to 400) in C88A, C91, C94 and C99A (Figure (Figure6;6; Table Table2)2) [5,13,75]. These clusters are likely due at least in part to two ABRE-associated TFs induced at phase-3 (Os02g52780: AB15, Os01g64000: ABF). In addition to ABRE-like motifs, C91, C94, and C99A were also enriched with as1/ocs/TGA-like elements and all three clusters had large proportions of H2O2 induced genes (67%, 70%, 64%, respectively). These clusters were activated in two phases (early and late), which are likely due to early and late activities of bZIP-TGA and bZIP-ABF genes, respectively (Figure (Figure2).2). WRKY factors appear to be involved in modulation of ABA-mediated expression, particularly in C91 which is enriched with W-box-like motifs. Some of the as1/ocs/TGA-like motifs in C94 and C99A were also quite similar to W-box-like elements, which could be functioning as WRKY target elements . WRKY genes could be regulatory fine-tuners by acting as repressors of ABA-induced genes at phase-1 and phase-2 .
Enrichment of ERF associated motifs in several clusters was consistent with the activities of several ERF genes belonging to groups I, II, III, IV, VI, VII, VIII, IX and X . ERF associated motifs were divided into two classes. The first includes DRE/CRT-like and rav1-like motifs, which are typically found among abiotic stress induced genes controlled by DREB/CBF, RAV1 and RAP2 [3,5,6,27]. This class of motifs occurred in 10 clusters but their enrichment was relatively subtle compared to other non-ERF-type motif classes (Figures (Figures6,6, ,7;7; Table Table33).
Clusters C78, C81, C86, C97A and C90A were most enriched with DRE/CRT-like and rav1-like elements (scores = 61 to 130), while C79, C80, C83A, C86 and C99 had intermediate (scores = 50 to 60) enrichment (Figures (Figures6,6, ,7).7). Clusters C78, C79, C80, C83A, C86, C90A and C97A were activated at phase-1 and this is likely due to several Group-II (Os06g11940: RAP2.1/2.6/2.10 group), Group-III (Os02g45450: DREB1A/CBF3, Os08g43210/Os09g35010: DREB1B/CBF1, Os01g73770: DREB1C/CBF2), Group-VII (Os12g40960: AP2 protein, Os03g08500: ERF2-like) and RAV1-like (Os01g04750: ABI3/VP1) genes that were also activated at phase-1 (Figure (Figure2).2). Clusters C80 and C81 were activated at phase-2, which are likely due to a different set of Group-II (Os03g15660: RAP2.1/2.6/2.10 group), Group-III (Os06g03670: DREB1A/CBF3) and Group-VII (Os03g08490: ERF2-like) genes distinct from those acting on phase-1 genes. Two genes from Group-III (Os09g35020: DREB1D/CBF3) and RAV1-family (Os01g04800: ABI3/VP1) appear to be involved in 'late response' regulation in C99. Many DREB/CBF-regulated genes are also activated through ABA-mediated pathway . Thus, the Group-I (Os04g44670: AtERF053, Os06g11860: AtERF058, Os03g09170: AtERF058), Group-Group-IV (Os05g28350: AtERF052/ABI4) and Group-X (Os09g28440) genes associated with ABA response may also be involved in the activation of DRE/CRT/rav1-like element-enriched clusters in response to ABA.
The second group of ERF associated elements includes the GCC-box-like and jasmonic acid response element-like (JAre) motifs (Figures (Figures6,6, ,7;7; Table Table3)3) involved in disease and wounding responses and regulation of indole alkaloid metabolism [37-39]. These motif classes were particularly enriched in C81, C83, C88A, C91, C94 and C100A. Clusters C88A, C91 and C94 were activated at phase-1 and this can be attributed to a number of group-VIII (Os04g52090: AtERF4, Os03g08500: AtERF2, Os06g47590/Os04g57340: AtERF3) and group-IX (Os08g44960: AtERF15, Os09g39810: PTI5/AtERF14) genes. Another gene from Group-VIII (Os05g41780: AtERF4) likely contributed to the activation of C81, C83 and C100A (Figure (Figure33).
ERF genes are also likely components of oxidative-mediated regulatory mechanism. Genes induced by H2O2 belong to Group-I (Os06g11860/Os03g09170: AtERF058, Os04g44670: AtERF053), Group-II (Os06g11940: RAP2.1/2.9/2.10 group), Group-III (Os01g73770: DREB1C/CBF2, Os09g35020: DREB1D/CBF4), Group-IV (Os05g28350: AtERF052/ABI4), Group-VII (Os12g40960: AP2 protein, Os03g08490: ERF2) and Group-X (Os09g28440: B-3 group). Possible targets of these TFs are in C94, C91 and C100A. These clusters were particularly enriched with H2O2-induced genes and had significant enrichment of GCC-box/JAre-like and/or DRE/CRT/rav1-like motifs (Figures (Figures6,6, ,77).
Potential MYB-target sequences were ubiquitous in almost all clusters implying that MYB genes are functioning either singly or in synergy with other TFs (Figures (Figures6,6, ,7;7; Table Table4).4). MYB-related motifs are quite diverse, but can be classified into two classes. MYB2-box-like motifs were significantly over-represented (scores = 50 to 250) in C79, C86, C87A, C88A, C90A, C91, C94, C97A and C100A (Figures (Figures6,6, ,7;7; Table Table4).4). Their high occurrence among chilling induced genes appears to be linked to the activities of R2R3-type MYB, which have been shown to directly interact with MYB2-box-like elements in the promoters of osmotic, drought and ABA induced genes [50,52,79-81]. Several OsMYB2 (Os01g18240, Os05g04820) and OsMYB4 (Os04g43680, Os02g41510, Os10g33810) homologs are likely the major players in the activation of C79, C86, C87A, C88A, C90A, C91, C94, C97A and C100A, perhaps in response to ABA, H2O2 or both. The second group of MYB-target elements includes several species of GA response element-like (GARE) and pyrimidine-box-like motifs, which are typically found among gibberellic acid (GA) regulated genes . These elements were particularly enriched (scores = 60 to 500) in C78, C79, C80, C81, C83, C83A, C86A, C87A, C88A, C90A, C91, C94, C97, C97A, C99, C99A and C100A. Recent studies have shown that R1-MYB particularly the SHAQKYF-type activates their GA-responsive targets through GARE and pyrimidine box-like elements .
R1-MYB genes are generally involved in seedling development and growth and they are known to regulate genes encoding hydrolytic enzymes involved in carbohydrate breakdown (amylases, glucanases, carboxypeptidases) for energy maintenance during early seedling growth [82,83]. Consistent with this, a number of R1-MYB genes were induced at 10°C including few SHAQKYF-type (Os06g07640, Os08g04840). Activities of these and other related R1-MYB genes are likely responsible for activating the GARE or pyrimidine box-like element-enriched clusters. Seedlings at the early stages of growth (V3) were used in the microarray studies, during which the endosperm was only partially degenerated. Growth at this stage requires mechanisms to integrate stress and developmental signals. The GARE or pyrimidine-box-like element-enriched clusters may be important for the integration of early growth and stress-related energy demands, which have been suggested as the key aspects of seedling vigor expression .
Broad distribution of MYB-related elements among chilling upregulated genes has important regulatory and physiological implications. MYB-target elements often occur in combination with one or more other types of element (Figures (Figures6,6, ,7),7), which implies that MYB factors are key components of combinatorial control mechanisms in conjunction with other major classes of chilling induced TFs (bZIP, ERF, WRKY) either as activators or repressors [49,50]. MYB-regulatory clusters may be involved in both stress and growth related processes. Within the context of seedling cold tolerance, MYB factors may be involved in coordinating the expression of genes required to modulate growth related processes at sub-optimal temperature.
Our current data indicate that regulatory modules involving bZIP Group-D, Group-I and Group-S (as1/ocs/TGA-like clusters), bZIP Group-A (ABRE-like clusters), DREB/CBF, RAP2 and RAV1 (DRE/CRT/rav1-like clusters), ERF (GCC box/JAre-like clusters), R1-MYB (GARE/pyrimidine-box-like clusters), R2R3-MYB (MYB2-box-like clusters), and WRKY (W-box-like clusters) TFs define the overall chilling stress response regulome of japonica rice. Based on these information, hypothetical models depicting the integration of various network components were established reflecting the complexity and hierarchical interaction among the various TFs acting on a given subset of chilling induced genes (Figure (Figure9).9). Based on the cis-element enrichment trends shown in Figures Figures66 and and7,7, some clusters involve relatively simple activation mechanism with one dominant class of regulatory genes. Examples are C83A, C91, C97A and C99A, which were primarily controlled by as1/ocs/TGA (bZIP), ABRE (bZIP), GARE/pyrimidine box (R2R3-MYB) and as1/ocs/TGA (bZIP) regulatory modules, respectively. Other clusters involved complex mechanisms with two or more dominant classes of TFs that may either be activated independent of each other or functioning as components of combinatorial control mechanisms. Examples are C79, C86, C90A and C100A, which involved as1/ocs/TGA (bZIP) + GARE/pyrimidine box (R2R3-MYB), as1/ocs/TGA (bZIP) + Myb2-box (R2R3-MYB), as1/ocs/TGA (bZIP) + GARE/pyrimidine box (R2R3-MYB) + W-box (WRKY), and as1/ocs/TGA (bZIP) + GCC-box/JAre (ERF), respectively. Many examples of the combinatorial control are well known particularly those involved in ABA responses [26,77,85,86]. At the center of the chilling stress regulatory network is ROS (H2O2) which acts as primary signal triggering the activities of 'early response' expression clusters. For instance, as1/ocs/TGA-like motifs are dominant in C83A and this cluster was also enriched with H2O2-responsive genes (65%). Therefore, an oxidative signal-mediated bZIP-TGA regulon appears to be a major component of this cluster. Similarly, C97A appeared to be regulated by R1-MYB. This cluster only had moderate enrichment of H2O2-responsive genes (Figure (Figure6),6), thus an oxidative-independent regulon is likely to be part of this cluster.
Many of the putative ERF-regulated clusters represent examples of complex regulatory mechanisms (Figure (Figure9).9). Clusters C81 and C86 for instance are enriched with DRE/CRT-like or rav1-like motifs and one or more other class of elements. Both clusters have low proportion of H2O2-responsive genes (Figures (Figures6,6, ,7),7), thus DREB/CBF regulon is likely to be a component of these clusters. Genes belonging to these clusters are likely to be regulated by combinatorial mechanism involving other TFs such as MYB for C81 and bZIP-TGA and MYB for C86 (Figure (Figure9).9). Enrichment of DRE/CRT-like and rav1-like elements was less pronounced compared to other types of elements (Figures (Figures6,6, ,7,7, ,9).9). It is possible that milder cold stress activates only a fraction of the entire regulon that would otherwise be induced under more severe cold conditions such as those that induce CA in temperate plants. It has been suggested that differences among plant species with respect to DREB/CBF regulon could be due to limited distribution of DRE/CRT-like elements in the genomes of chilling sensitive species and this might be true in rice .
Unlike bZIP, ERF and MYB-related elements, motifs associated with WRKY, NAC and bHLH TFs had lower enrichment among the chilling induced genes (Figures (Figures6,6, ,7;7; see Additional file 7). W-box-like enrichment in C90A appears to be associated with the function of WRKY factors as repressor of ABA-inducible expression. Preliminary analysis of chilling downregulated genes indicate that putative WRKY, MYC and NAC-target elements were among the enriched motifs (data not shown). Thus, most WRKY, NAC and MYC genes are likely functioning primarily as repressors of chilling stress response gene expression [54,59,60,87].
Of the 2,456 chilling upregulated NTF genes, 93% (2,295) could be assigned to one or more molecular or biochemical process based on integrated analysis of Gene Ontology  and protein domain i.e., INTERPRO  databases (see Additional file 5). Nearly 30% (699) were described with various types of 'cellular response' or 'stress' related keywords including defense (GO:0006952), stress (GO:0006950), water (GO:0009725), wounding (GO:0009611), oxidative (GO:0006979), cold (GO:0050826, GO:0042309, GO:0050825) and various stimuli including biotic (GO:0009607), abiotic (GO:0009628), external (GO:0009605) and hormone (GO:0009725). Three functional categories reflecting the major processes involved in chilling stress response mechanism were established.
Processing and transduction of environmental signals to the nucleus constitute the initial events in the execution of adaptive responses [4,7,62]. About 20% (494) of the upregulated NTF genes are associated with signaling and response regulatory processes such as protein phosphorylation (GO:0004674, GO:0004713), receptor recognition (GO:0004872, GO:0007154, GO:0004871), Ca2+ signaling (GO:0005509, GO:0005544), small GTPase coupled signaling (GO:0007186, GO:0007264), and generation of second messengers (GO:0046488, GO:0007165, GO:0004435) (Figure (Figure10;10; see Additional file 5).
A significant increase in the rate of mRNA synthesis was implied by genes required for the assembly of basal transcription machinery including those involved in transcription initiation and elongation, tRNA synthesis, and transcript processing (GO:0003700, GO:0003723, GO:0003712, GO:0016439, GO:0003677, GO:0003899, GO:0030528, GO:0003702). Enhanced activities of these genes reflect an increased rate of de novo mRNA synthesis required for activation of response-related transcriptome. Nearly 60% of genes associated with signaling and response regulation were also induced by H2O2, further supporting the central role of oxidative signaling in chilling stress response mechanism. Although the genes in this category are widely distributed among 18 clusters, phase-1 clusters tend to be more enriched than phase-2 and phase-3 clusters, further indicating that critical signaling events occur during the initial 6 hours.
About 30% (769) of total NTF genes are associated with cellular defense and rescue processes. Genes involved in biogenesis or repair of damaged cellular components comprised about 14% of total NTF genes (Figure (Figure10;10; see Additional file 5). Dominant classes were chaperones, heat shock proteins and proteases involved in rescue or turnover of damaged proteins (GO:0006508, GO:0006511, GO:0016567, GO:0000151, GO:0030693, GO:0004185), proteins involved in membrane and cell wall biogenesis (GO:0008654, GO:0045300, GO:0006869, GO:0042546, GO:0042545) and DNA repair (GO:0006284, GO:0006281, GO:0003684) .
Genes associated with oxidative stress including cellular redox regulators (GO:0045454, GO:0004602, GO:0004601, GO:0006801, GO:0006979; see Additional file 8) and pathogen defense and apoptosis-related proteins (GO:0009607, GO:0009605, GO:0009405, GO:0004568, GO:0006952, GO:0006915, GO:0042742; see Additional file 9) represent another major class of defense and rescue associated genes with about 11% of total NTF genes. Cellular redox regulatory processes include the components of glutathione system, respiratory burst proteins, and ROS-scavengers such as peroxidase, superoxide dismutase, catalase, ascorbate peroxidase, thioredoxins, and metallothionein-like proteins . Genes associated with defenses against pathogens include a large number of NBS-LRR-type disease resistance proteins and various pathogenesis-related (PR) proteins such as chitinases, glucanases, thaumatin-like, germin-like, cupin-like, and dirigent proteins, viral, bacterial, and elicitor-induced proteins, and proteins involved in phenylpropanoid, salicylic acid, jasmonic acid and phytoallexin metabolism [90-94]. This trend shows striking similarities between chilling-induced defenses and those associated with responses to pathogens and wounding, all of which are linked to H2O2 signaling [75,95].
About 2% of genes associated with cellular defense and rescue are involved in detoxification and efflux processes (GO:0006855, GO:0015904, GO:0006904, GO:0006887, GO:0006814) (Figure (Figure10;10; see Additional file 10). Various classes of multi-drug resistance, MATE efflux and other types of proteins involved in the extrusion or sub-cellular sequestration of toxic compounds were in this group [96,97]. Many of these genes are also involved in defenses against pathogens. Based on the expression of osmotins and aquaporins and other genes involved in galactinol and trehalose biosynthesis, maintenance of root and leaf water balance and cellular osmotic adjustment (GO:0005992, GO:0009719; see Additional file 5) are important components of defense mechanism [98-100]. Genes associated with thermal regulation such as homoiothermic/freezing-associated proteins (GO:0042309, GO:0050825, GO:0050826) were also induced with more than 2% of total NTF genes (see Additional file 5).
A small group of novel stress-induced genes were also assumed to be part of the defense and rescue category, with nearly 1% of total NTF genes. These genes are known only by virtue of their responsiveness to developmental and environmental stress factors such as cold, drought, salinity, wounding, hypoxia, ABA and senescence (see Additional file 5). Included in this group are several homologs of stress-related genes belonging to the early responsive to dehydration (ERD) and LEA/dehydrin families.
Nearly 60% of genes associated with cellular defense and rescue mechanism were induced by both chilling and H2O2. Many involved in redox regulation, pathogen defense, apoptosis and proteolysis (Figure (Figure10).10). This trend implies that defense processes are essentially consequences of chilling-induced oxidative stress and this is consistent with the fact that abiotic and biotic stresses induce a common set of genes through an oxidative-mediated pathway. Cellular defense and rescue related genes are widely distributed in all of 18 clusters. However, phase-1 clusters tend to be more enriched with this functional category than phase-2 and phase-3 clusters (see Additional file 6) implying that short-term adaptation depends on timely activation of defense and rescue mechanisms.
Maintenance of cellular homeostasis is critical for sustained growth and survival under sub-optimal temperature. Cellular physiology must be adjusted to accommodate stress-related demands and metabolic requirements. Means to compensate or replenish metabolic intermediates that have been diverted towards costly defense and rescue processes are necessary for physiological sustenance. The largest category of NTF genes represents diverse molecular and biochemical functions suggesting the nature of physiological adjustment and sustenance processes. This category accounts for more than 40% (1,032) of the total NTF genes (see Additional files 5, 11).
Increase in de novo protein synthesis was indicated by genes involved in ribosome assembly and translation (GO:0006412, GO:0003735, GO:0008135, GO:0006446) and protein folding and modification (GO:0006464, GO:0019538, GO:0006457). Genes involved with cellular energetics (GO:0006118, GO:0006096, GO:0006094, GO:0015979, GO:0006091), transport mechanism and facilitation (GO:0006810, GO:0005215, GO:0006865, GO:0015986, GO:0016469, GO:0008643, GO:0046907), carbohydrate, lipid, amino acid and secondary metabolic processes (GO:0008610, GO:0006633, GO:0008299, GO:0006629, GO:0009058, GO:0005975, GO:0019748, GO:0006807), cellular biogenesis (GO:0005875, GO:0005856) and growth (GO:0006260, GO:0016575, GO:0006333, GO:0006334) were highly represented in this category. Genes associated with physiological adjustment and sustenance are widely distributed among 18 expression clusters but more highly represented in the slower clusters (phase-2, phase-3). This trend suggests that physiological adjustment and sustenance processes occurred largely during the middle to later stages of stress (see Additional file 6) and reflects the biochemical costs of regulatory and defense related activities during the earlier stages of stress. It also suggests that 'late response' genes are important for recovery process.
Stress tolerance potential is governed by quantitative trait loci (QTL) and TFs are likely to be major components of QTL. To investigate possible association between stress tolerance QTL and the major activators of oxidative-mediated chilling stress response, genomic locations of chilling upregulated TFs were determined in relation to genomic boundaries of relevant QTL that have been anchored to the Nipponbare genome sequence . The ratio of QTL-associated to non-QTL-associated genes in the total annotated gene set in the microarray (whole genome) was determined by Fisher exact test and compared with the corresponding ratios within the upregulated gene subset. Significant enrichment was defined by a higher ratio in the upregulated gene subset (relative to genome-wide ratio) at p < 0.25. A total of 42 chilling induced TFs were located within the genomic boundaries 46 QTL for cold tolerance (COLDTL), seedling vigor (SDLVIG), osmotic adjustment capacity (OSADJCAP), salt sensitivity (SALTSN), chlorophyll content (CHLCN), leaf rolling (LFRL), and speed of germination (GERMSP) (see Additional file 12). SDLVG, CHLCN, COLDTL and OSADJCAP had the most number of co-localized TFs with 24, 16, 12, and 10 genes, respectively.
We also examined the differences in expression of representative TFs between Nipponbare and INIAP12 rice cultivars (Figure (Figure11).11). These cultivars represent contrasting sensitivity to chilling based on seedling survival test (Figure (Figure12).12). Chilling-induced expression of TFs in Nipponbare occurred earlier and was more robust than INIAP12, consistent with results from previous studies [14,17]. Differential expression of TFs positively correlates with the responses of japonica and indica cultivars to chilling. About 4 weeks after continuous exposure to 10°C, 50% of indica seedlings were either dead or with severe injury symptoms. Mild injuries occurred in japonica seedlings only after about 6 weeks. Based on these results, stress response transcriptome is tightly associated with the expression of vigor.
We studied the regulatory and physiological implications of chilling stress (10°C) transcriptome of japonica rice based on integrative analysis of whole-genome expression profiles and promoter architectures. Our analysis has taken into consideration the hypothesis that oxidative signaling is central in integrating various components of the regulatory network. We addressed this hypothesis by dissecting the commonalities between a direct response to cold stimulus and responses to low temperature-independent mimic of oxidative signaling. From the systems-level approach employed in this study, several important themes emerged. First, oxidative signaling by H2O2 is at the center of the regulatory network, particularly in relation to the execution of 'early response' mechanisms (initial 24 hours). H2O2 has long been recognized as a major trigger of stress response signaling and its important role in interfacing abiotic and biotic stress responses with growth and development has been established [67,68,101,102]. Until now, direct links of the primary oxidative signals to chilling-associated transcriptional changes has not been demonstrated at the genome-wide scale. Our current results contribute important information that fills this knowledge gap.
About 60% of chilling stress induced genes are triggered by oxidative signals either as primary or secondary targets. Much of the oxidative-mediated changes in gene expression occurred during the initial 24 hours and the composition of the transcriptome has striking similarities to disease response mechanisms in terms of the activities of genes involved in redox regulation, responses to organic toxins and wounding, phenylpropanoid and indole alkaloid metabolism, and jasmonic acid, salicylic acid and ethylene signaling [39,62,73]. H2O2 is likely to be the primary signal closest to the original stimulus (chilling) providing the initial trigger for a cascade of molecular processes leading to an intricately interconnected gene expression circuitry [101,102]. The timing and scope of oxidative mediated regulatory clusters justify their major contribution to immediate defenses, which are critical for a species like rice that can only endure transient exposure to milder cold stress.
Second, our current data revealed a number of oxidative-mediated early response expression clusters associated with specific classes of regulatory sequences such as as1/ocs/TGA-like, GCC-box/JAre-like and Myb2-box-like cis-elements. Several classes of TFs are likely to act on these clusters including bZIP Group-D (TGA), Group-S (ocs), and Group-I (GA-associated), ERF Group-I (subfamily A-6), Group-IV (subfamily A-2), Group-VI (subfamily B-5), Group-VII (subfamily B-2), Group-X (subfamily B-4), and R2R3-MYB [24,29,37,60]. Known interactions of these TFs with such elements and their involvement in various oxidative mediated processes justify their presumed roles as regulators of early response transcriptome. Particularly interesting candidates for regulon engineering are the as1/ocs/TGA-like element-regulated clusters comprising the largest group of genes associated with oxidative defenses.
Third, we have identified regulatory clusters that appear to be independent of oxidative signals. Two of the well known regulons that function during cold acclimation at 4°C in Arabidopsis were among those and they are the ABRE-like (bZIP-ABF) and CRT/DRE/rav1-like (CBF/DREB/RAP2/RAV1) element-enriched clusters. The third cluster (GARE/pyrimidine-box-like element-enriched) which involves an R1-MYB appears to be a component of combinatorial control mechanism. ABRE-enriched clusters appear to be downstream to the oxidative-mediated clusters and this is consistent with the timing of ABA response, i.e., after 24 hours. Recent report estimated that 400 protein-coding genes of rice are regulated via an ABRE requiring mechanism . This number is close to the estimated number of genes that belong to the chilling induced ABRE-enriched clusters, which is much less than the estimated number of genes regulated by an oxidative mediated mechanism. This data provide additional support that oxidative-mediated mechanism is the primary route for activating early response genes.
The cold acclimation associated regulons controlled by DREB/CBF, RAP2 and RAV1 are induced maximally at 4°C in Arabidopsis [13,42,43]. Apparently, orthologous TFs in rice are induced at 10°C, which probably induce a different set of target genes distinct from those that are activated in cold acclimating species. An interesting trend that we observed was that chilling upregulated clusters in rice have only moderate enrichment of DRE/CRT/rav1-like elements compared to other types of elements. The DRE/CRT/rav1-like elements associated clusters are also enriched with one or more other types of elements. A possible explanation is that at milder cold stress, DREB/CBF/RAP2/RAV1-regulated genes account for a relatively smaller fraction of the transcriptome compared to the oxidative mediated components. It has also been suggested that the genomes of chilling insensitive (e.g., Arabidopsis) and sensitive (e.g., rice) species might differ in terms of the distribution of DRE/CRT/rav1-like elements .
Fourth, our data showed an important role of combinatorial control in chilling stress response gene expression. Combinatorial control has long been recognized and numerous examples are constantly being reported in relation to stress response [17,26,77,104]. Interaction of several classes of TFs facilitates fine-tuned regulation. Examples of these are the possible roles of WRKY TFs as repressors of ABA-mediated gene expression, which allows fine-tuned regulation of oxidative mediated transcriptome in relation to ABA-mediated responses. Interaction of several TF classes also facilitates integration of various signals associated with stress, growth and development. GARE/pyrimidine-box-like elements associated with R1-MYB TFs were the most widely distributed elements among the chilling upregulated clusters. This class of elements occurred in various combinations with other classes. R1-MYB TFs are primarily involved in gibberellic acid regulated growth processes [79,82,83]. R1-MYB factors appear to coordinate and fine-tune gene expression required for growth modulation and survival at sub-optimal temperature, which appear to be a key aspect of vigor enhancement and/or maintenance under stress conditions.
Finally, our results provide a global picture of biochemical and physiological consequences of the regulatory networks induced at 10°C. The obvious theme is the striking similarities of chilling-induced processes to those associated with defenses against pathogens and responses to wounding. The predominance of genes involved in disease response signaling was indicated by large number of LRR, NBS-LRR, NB-ARC types of receptor proteins, genes associated with hypersensitive response, wounding, herbivory, redox regulation, cellular detoxification, programmed cell death, and terpenoid metabolism most of which are linked to jasmonic acid, salicylic acid and ethylene signaling. Activities of these genes provide further evidence that oxidative signal was the origin of most of the early responses.
The agronomic significance of the current results is quite interesting. The data provide a global picture of how stress, growth and developmental responses are integrated by the interaction of at least four major classes of TF. The interesting relationship between TF expression and QTL for stress tolerance and seedling vigor might be the key for further understanding of the precise mechanisms that make japonica rice withstand chilling better than indica rice. Genes related to stress response, growth, development and energy partitioning are often enriched within the boundaries of QTL associated with vigor and yield. An emerging theme is that linkage blocks with certain allelic combinations of stress-related genes involved in responses to diseases, dehydration, temperature extremes, and light stress are important components of plant vigor and heterosis [18,105-108]. Modulated growth, physiological sustenance and maintenance of vigor appear to be the downstream consequences of the transcriptional networks induced by chilling, providing a means to delay the occurrence of irreversible injuries. This hypothesis is consistent with the fact that many japonica cultivars are able to survive chilling for a much longer period of exposure than most indica cultivars.
Rice cultivars with contrasting sensitivities to chilling, i.e., Nipponbare and CT6748 (less sensitive) and INIAP12 and IR36 (more sensitive) were used in this study. Seedlings were grown to three-leaf (V3) stage at optimum temperature (28°C) before exposure to chilling condition (10°C) in a Pervical E30BHO growth chamber (Percival Scientific, Perry, IA). These experiments were performed based on previously described methods [14,17]. For the H2O2 experiment, seedlings were grown to V3 stage in standard Yoshida hydroponic solution at 28°C. Experimental plants were transferred to fresh medium with 4 mM H2O2 while the control plants were maintained in medium without H2O2. Tissues were harvested after 1, 3 and 6 hours of treatment and another 6 hours of recovery in fresh medium without H2O2. Total RNA was isolated with Trizol reagent (Invitrogen, Carlsbad, CA) from leaves after 0.5, 2, 4, 6, 12, 18, 24, 36, 48 and 96 hours of exposure to chilling and after each sampling time in the H2O2 experiment. Seedling survival analysis was performed as previously described .
RNA samples from control and experimental Nipponbare seedlings were processed with the MessageAmp II aRNA amplification and T7 polymerase in vitro transcription kits producing at least 20 μg of aminoallyl-dUTP-labeled aRNA. The cDNA samples were synthesized with oligo-dT primer and ArrayScript reverse transcriptase. All procedures were according to manufacturer's instructions (Ambion, Austin, TX). Equal amounts of Cy5- and Cy3-labeled samples were combined in 1× hybridization buffer composed of 50% formamide, 5× SSC, 0.1% SDS and 10 mM DTT. Chemical labeling with Cy3 (control) and Cy5 (treatment) dyes was performed according to manufacturer's protocol (GE Healthcare-Amersham, Piscataway, NJ).
The NSF rice oligonucleotide microarray version 3, which contains roughly 45,000 probes representing all 40,000 predicted genes of japonica rice was used in all experiments. Pre-hybridization treatments were according to manufacturer's instructions . Hybridization was performed for 18 hours at 42°C in Hybex Humidified Thermal Blocks (Scigene, Sunnyvale, CA). Stringency washing was performed at 42°C in 2× SSC, 0.1% SDS and 0.1× SSC, 0.1% SDS and at room temperature in 0.1× SSC for 10 minutes each step. Expression data was acquired with the Axon 4000A scanner and processed with the GenePix Pro 5.1 (Axon-Molecular Devices, Sunnyvale, CA). Two independent biological replicates were performed for each control by treatment comparison. Dye-swap experiments were initially performed to assess dye by sample interaction, which was found to be negligible. Gene expression data was normalized globally prior to high level analysis. Microarray dataset can be accessed by GSE8767 and GSE10062 at the Gene Expression Omnibus . Values represent the log2 ratio of background subtracted intensity values (Cy5/Cy3). High level statistical analyses (T-test, hierarchical and K-mean clustering) were performed with the MeV Bioinformatics Tools . Protein-coding genes with log2 fold ≥ 1.8 (p < 0.05) in at least two consecutive time points for the chilling stress experiment and at least one time point for the H2O2 experiment were first identified from the total dataset. Genes that passed these criteria comprised the total gene set for high-level analysis. TFs were identified and grouped according to family according to the Database of Rice Transcription Factors . Members of each family were hierarchically clustered while the NTFs were divided into smaller groups by K-mean clustering. NTFs were further classified by functional categories using the RiceCyc pathway  and Interpro protein domain  databases. Co-location and enrichment of candidate genes within the boundaries of relevant QTL was performed with Fisher exact test based on current information in Gramene QTL database .
Promoter motif enrichment analysis was performed on the entire set of upregulated NTFs (2,456) and individually for each K-mean cluster . Sequences of bona fide promoter regions (-1,000 to +200) were extracted from the Nipponbare genome sequence by locating the experimentally validated transcription start site (TSS) by alignment with FLcDNA . The Dragon Motif Builder algorithm with EM2 option was used to detect over-represented motifs . Thirty motifs (8 to 10 nt) were detected each run with a threshold value of 0.875. Promoters of randomly selected genes not associated with stress response mechanisms were used for background subtraction of random motif occurrence. Significant motifs were selected based on a threshold occurrence of 50%. Motif classes were identified by significant matches with TRANSFAC [113,114], PLACE  and AGRIS [116,117] databases.
Comparative analysis of selected chilling upregulated genes between Nipponbare and INIAP12 was performed by real-time PCR as previously described . PCR reactions were performed with the Verso cDNA synthesis and Absolute QPCR Sybr Green mix (AbGene, Rochester, NY) using the myCycler real-time PCR system (Biorad, Hercules, CA).
Total H2O2 content of rice leaves was determined by the Amplex Red (10-acetyl-3,7-dihydroxyphenoxazine) assay (Invitrogen, Carlsbad, CA) according to established procedures . Leaf discs were homogenized in chilled 5% TCA. Crude extracts were centrifuged at 12,000×g for 10 minutes and further purified with Dowex anion exchange resin (AG1X100). Samples were buffered to pH 7.4 with 0.25M NaH2PO4. Concentration was measured every 30 minutes by the absorbance at 560 nm with three replicates.
Nuclear protein extraction was performed based on established procedures [118,119]. Proteins were extracted from leaf tissues (10 g) after exposure to chilling (10°C) for 6, 12 and 24 hours and quantified by Bradford method. EMSA was performed with the Gel Shift Assay System (Promega, Madison, WI) following the manufacturer's protocol. Double stranded oligonucleotide probes were synthesized and end-labeled with [γ-32P]-dATP (3,000Ci/mmol) by T4 polynucleotide kinase. Unlabelled probes were removed using the Bio-Spin 6 (Bio-Rad, Hercules, CA). Binding reactions (10 μg nuclear extract +1 μl 32P-labeled probe) were carried out for 30 minutes at room temperature. For competition experiments, excess unlabeled probes were added in approximately eight-fold molar ratio relative to labeled probes. Binding reactions were analyzed by electrophoresis in a 4% non-denaturing polyacrylamide gel.
KYY and MRP performed the microarray experiments including all statistical and bioinformatic analyses. RM and RB performed the gene ontology (GO) and QTL enrichment analyses. BM, EW and VBB performed the ab initio promoter analysis. VH contributed to the microarray and RT-PCR analyses and performed the leaf H2O2 quantitation experiments. FX performed the genotypic comparisons by RT-PCR and EMSA experiments. BGDR was responsible for the overall concept and experimental designs, data integration, analysis and interpretation, and manuscript preparation. All authors approved the final manuscript.
Dominant functional categories in the downregulated group of genes. This data shows the most highly enriched broad functional categories of downregulated genes classified according to gene ontology.
Expression matrix of chilling upregulated WRKY transcription factors. Heat map showing the temporal expression profiles of WRKY transcription factors under chilling stress. Gene designations were based on putative Arabidopsis orthologs according to the most recent genome annotation.
Expression matrix of chilling upregulated bHLH transcription factors. Heat map showing the temporal expression profiles of bHLH transcription factors under chilling stress. Gene designations were based on putative Arabidopsis orthologs according to the most recent genome annotation.
Expression matrix of chilling upregulated NAC transcription factors. Heat map showing the temporal expression profiles of NAC transcription factors under chilling stress. Gene designations were based on putative Arabidopsis orthologs according to the most recent genome annotation.
Functional categories and grouping by cluster of chilling upregulated NTF genes. List of all NTF genes included in the analysis of the chilling upregulated transcriptome grouped according to putative functions.
Correlation between timing of induction and function of chilling responsive genes. Distribution of functional categories in relation to activation timing of co-expressed gene clusters. (A) Signaling and response regulation; (B) cellular defense and rescue; (C) physiological adjustment and sustenance mechanisms. Gray: rapid phase-1 clusters; Red: phase-2 and phase-3 clusters.
Relative enrichment of other types of cis-elements detected in the promoters of chilling upregulated genes. List of cis-elements associated with other types of transcription factors not included in Tables Tables2,2, ,33 and and44.
Components of the chilling stress transcriptome with possible roles in oxidative stress and redox regulation. List of upregulated genes with possible roles in oxidative stress and redox regulation classified according to gene ontology.
Components of the chilling stress transcriptome with possible roles in response to diseases, elicitors and apoptosis. List of upregulated genes with possible roles in response to diseases, elicitors and apoptosis classified according to gene ontology.
Components of the chilling stress transcriptome with possible roles in cellular detoxification and efflux. List of upregulated genes with possible roles in cellular detoxification and efflux classified according to gene ontology.
Possible components of physiological adjustment and sustenance mechanisms. Functional categories relevant to physiological adjustment and sustenance processes classified according to gene ontology.
Association between chilling-induced transcription factors and stress-associated QTL. Genomic location of chilling upregulated transcription factors relative to the boundaries of known QTL of rice associated with seedling vigor and stress response.
This project was supported by USDA-CSREES-NRI Grant 2006-35604-1669, Maine Agricultural and Forestry Experiment Station (3042) and University of Maine-MEIF Program. M.R. Park was supported by a fellowship grant from the Korea Research Foundation (KRF-2006-352-F00002) and BioGreen 21 Program-Rural Development Administration (20080401034024), Republic of Korea. We thank Drs. Mitch McGrath and Yulin Jia for their critical review of the manuscript. We also thank the anonymous reviewers for their insightful reviews and suggestions.