Plants of japonica rice cv. Dongjin (‘DJ’) and indica rice ‘IR64’ (‘IR64’) were grown in a walk-in chamber (Koencon, Hanam, South Korea) under conditions of 30°C [200 μmol m^-2 s^-1 (day)]/22°C (night) and a 12-h photoperiod for 10 days in plastic boxes containing 100 g of soil used in growing rice (Punong, Kyung-Ju, Korea) (Kumar et al., 2017). The effects of cold stress (exposure at 4°C) on the light intensity 110 μmol m^-2 s^-1 were examined after exposure to cold stress for 0, 24 h/1 day, 48 h/2 days, 72 h/3 days, 96 h/4 days, 120 h/5 days, and 144 h/6 days using chlorophyll fluorescence. Our mock treatment comprised a group of plants that remained at the normal growing temperature (28°C) throughout the experimental period. To observe the physiological features of these seedlings, we used samples collected before cold stress was induced, as well as from plants after 4 days of stress, and then after recovery under normal conditions for 5 days. Fresh weights (FWs) were recorded after recovery from cold stress, and dry weights (DWs) were measured after the samples were dried at 80°C for 2 days.

For monitoring the expression of cold-inducible marker genes, seedlings (selected at 10 DAG, or 10 DAG) were hydroponically cultured in Yoshida solution and exposed to 4°C for 0, 1, 3, 6, 12, or 24 h. Primers of OsZFP182/LOC_Os03g60560 and OsWYRKY71/LOC_Os02g08440 were used for RT and qRT-PCR analyses at a final concentration of 10 pmol, with 3 μL (equivalent to 30 ng of total RNA) of cDNA as template (Supplementary Table S1). The internal controls were primers of rice ubiquitin 5 (OsUbi5) and rice actin 1 (RAc1) (Supplementary Table S1). An RNeasy Mini Plant Kit (Qiagen, Germany) was used for total RNA isolation and an RT Complete Kit (Biofact, Korea) was used for cDNA synthesis according to the manufacturers’ instructions. Primers were designed with Gene Runner software³ and NCBI primer blast⁴. The amplified products were resolved on a 1% agarose gel.

An uptake assay was conducted to determine the relative concentration of H₂O₂, using Amplex^® Red reagent (10-acetyle-3, 7dihydroxyphenoxazine; Molecular Probes/Invitrogen, United States) (Mohanty et al., 1997). Leaf tissues (0.1 mg μL^-1) were homogenized in a standard MS medium (Murashige and Skoog, 1962) and then incubated under darkness for 30 min with horseradish peroxidase (0.2 U mL^-1) and Amplex^® Red reagent (1 μM). The H₂O₂ released from these tissues was detected by a SpectraMax 250 Microplate Reader (Molecular Devices Inc., United States) with absorbance measured at 560 nm (Kumar et al., 2014).

We downloaded raw data for five GSE data series (i.e., GEO accession number GSE6901, GSE33204, GSE37940, GSE38023, and GSE31077) that are related to cold-stress responses, as indicated from the NCBI GEO⁵ (Barrett et al., 2011). Details are presented in Supplementary Table S2. The data were normalized using an Affy Package encoded by R language, and the intensity values were transformed into the log₂ scale as we have previously described (Cao et al., 2012). This allowed us to generate log₂ fold-change values for cold-stressed samples. Similar fold-changes were revealed for other stress conditions. For each data series, we used those fold-change data to perform a KMC analysis to identify genes that were consistently up-regulated under all cold-stress conditions. The KMC analysis of meta-expression data for abiotic stresses – salt, drought, cold, heat, submergence, and anaerobic conditions – grouped all of the candidate genes into 12 clusters. From these, we selected 502 genes that were up-regulated by cold-stress treatment but not during the recovery period. Heatmap images were produced using Mev software (Chu et al., 2008).

To examine GUS expression patterns, we germinated seeds from two promoter trap lines in an MS medium for 7 days. These lines were obtained from a mixed pool of PFG T-DNA tagging lines from POSTECH in Korea (Lee et al., 2004; Jung et al., 2005, 2006, 2015; Hong et al., 2017; Wei et al., 2017). The resultant plantlets were then exposed to cold stress (4°C) for 0 or 24 h. Afterward, whole seedlings from all treatment groups were soaked for 8 h in a GUS-staining solution before their roots were photographed with a camera (Canon EOS 550D; Cannon, Tokyo, Japan).

To identify any consensus CREs in the promoters of our cold-inducible genes, we extracted 2-kb upstream sequences of ATG for LOC_Os01g31370 and LOC_Os03g49830, which were validated in our current GUS assays. We also used the sequence for LOC_Os10g41200, which was previously reported to be a cold-inducible promoter based on the promoter-GUS system (Rerksiri et al., 2013; Jeong and Jung, 2015) from PLANTPAN⁶ (Chang et al., 2008). Several MEME searches were performed with those sequences in the FASTA format via the Web server hosted by the National Biomedical Computation Resource⁷. We looked for up to five CREs with an option of 12 maximum motif widths. Using the MAST, we then searched DNA sequences for matches to the putative TOMTOM within a set of promoter sequences (Bailey et al., 2006).

To analyze the biological significance of selected candidate genes, we employed the GO enrichment tool installed in the Rice Oligonucleotide Array Database⁸ (Jung et al., 2008a; Cao et al., 2012). For this, we uploaded 502 genes showing upregulation in both japonica and indica cultivars under cold stress. A fold-enrichment value higher than the standard (‘1’) meant that the selected GO term was over-represented more than was expected. Terms with >2-fold enrichment values and p-values < 0.05 were also used as criteria for choosing the most significant GO terms in the ‘Biological Process’ category.

The rice MapMan classification system covers 36 BINs, each of which can be extended in a hierarchical manner into subBINs (Usadel et al., 2005; Urbanczyk-Wochniak et al., 2006). By applying diverse MapMan tools, a significant gene list selected from high-throughput data analysis can be integrated to diverse overviews. Here, we uploaded locus IDs from the RGAP for 502 DEGs with a value of ‘3,’ which indicated upregulation under cold stress. Finally, we used four overviews – Metabolism, Regulation, Transcription, and Proteasome – installed in the MapMan toolkit.

To evaluate the functional significance of our candidate genes, we compared our list with the one from OGRO, which summarizes rice genes with known functions (Table 1; Yamamoto et al., 2012).

Plants from an Ox line for OsWYRKY71 (OsWYRKY71-Ox) under the control of CaMV35S promoter (Kim et al., 2016) and from the WT (Japonica cv. Dongjin) were grown for 10 days in plastic boxes containing soil. To test their tolerance, we then exposed them to cold stress (4^oC) for 5 days and then returned them to normal growing conditions for 6 days of recovery. Survival rates were determined at the end of this experimental period. Cold stress analysis of OsWYRKY71-Ox lines was done with three replicates.

Cold stress adversely affects plant growth and yield, and rice isparticularly susceptible at the seedling stage (Zhang et al., 2014). Ouranalysis involved 10-day-old ‘DJ’ (japonica) and ‘IR64’(indica) plants exposed to 4°C for 4 days.Afterward, they recovered for 5 days at 28°C. Their phenotypes are shown in Figure 1A. At the end of this experimental period, the survival rate was 30.5% for ‘DJ’ versus 0.0% for ‘IR64,’indicating that the former was more old-tolerant (Figure 1B). The F_W value was 162 mg higher for ‘DJ’ while its DW was 29 mg higher than for ‘IR64’ (Figure 1C). Prolonged cold stress also negatively affected photosynthetic efficiency, with both cultivars showing significant reductions after 24 h (Figure 1D). The decline in efficiency after 48 h was more severe for ‘IR64’ than for ‘DJ.’

The accumulation of ROS, including H₂O₂, is a major indicator of the plant response to various abiotic stresses. We found that ‘IR64’ had higher H₂O₂ concentrations than did ‘DJ’ after 3 and 24 h of cold treatment (Figure 1E).

We also evaluated the expression patterns of two well-known cold stress-responsible genes, OsZFP182 and OsWRKY71 (Huang et al., 2007; Kim et al., 2016) and found that, as expected, their expression was significantly induced, and to nearly the same extent, in both cultivars (Figure 1F). This demonstrated that the tool of global transcriptome data can be broadly applied for determining and, ultimately, improving cold tolerance in rice.

As a quantitative trait, tolerance to cold stress is governed by different sets of genes, and through diverse mechanisms. We used meta-expression analysis with transcriptome data and downloaded information about global candidate genes from the NCBI GEO for series GSE37940 and GSE38023 (Zhang F. et al., 2012; Zhang H. et al., 2012). After normalizing these data, we generated 63 comparisons for cold-stress treatment, as well as 49 comparisons for drought stress, 6 for high temperatures, and 4 for submergence (Supplementary Table S2). Our KMC analysis with the resultant fold-change data revealed 502 genes that were significant up-regulated upon cold stress but not under recovery conditions (Figure 2). From this, we prepared 27 comparisons with two japonica cultivars – ‘C418’ (a japonica restorer line for hybrid rice production and cold sensitive) and ‘Li-Jiang-Xin-Tuan-Hei-Gu’ (‘LTH,’ cold tolerant genotype) – and 36 comparisons with five indica cultivars – ‘IR24’ (photoperiod-insensitive, high yielding and cold sensitive variety), ‘IR64’ [variety with moderate tolerance toward toxicity in response to various molecules including salt, alkali, iron, and boron as well as deficiencies in phosphorus and zinc, but sensitivity to cold stress], ‘K354’ (a BC2F6 introgression line as a progeny of C418 and cold tolerant variety), ‘Huahui 1’ (’HH1,’ insect-resistant variety as a progeny of Minghui 63), and ‘Minghui 63’ (‘MH,’ heat tolerant variety and a parental line of HH1). Their upregulation was conserved between japonica and indica cultivars. All of these genes provide potential for a broader range of applications to enhance cold tolerance in rice. These 502 DEGs were used for further analysis of the cold-stress response (Supplementary Table S3).

Promoter traps employing the GUS reporter gene system have been used to identify promoters involved in regulating tissue-specific and stress-responsive expression patterns (Jung et al., 2005, 2006). Our meta-expression analysis identified the top 50 genes showing >3.5-fold upregulation by cold stress when compared with the control (Supplementary Table S3). We then searched and identified 52 potential promoter trap lines of 43 candidate genes and examined GUS expression patterns in 7-day-old seedlings. The lines for two genes (PFG 3A-50649 for LOC_Os01g31370 and PFG 1C-08613 for LOC_Os03g49830) displayed GUS expression in the roots after plants were exposed to stress for 24 h (Figure 3 and Supplementary Figure S1). This cold-related expression was also verified by qRT-PCR (Supplementary Figure S1). Previous studies using a promoter-GUS vector or promoter trap system have confirmed the upregulation of LOC_Os10g41200 in response to cold stress (Su et al., 2010; Jeong and Jung, 2015). Our findings demonstrated that the promoter trap system, when combined with qualified genome-wide transcriptome data, is a very effective way for quickly identifying the activity of an endogenous promoter. This enables researchers to develop novel promoters for future applications.

To identify the cis-regulatory regulatory elements (CREs) associated with the response to cold, we used promoter regions in 2-kb sequences upstream of ATG of the two cold-inducible genes (LOC_Os01g31370 and LOC_Os03g49830) that had been validated through GUS trap assays and also included the promoter region of LOC_Os10g41200, which have previously been reported as a cold-inducible gene using GUS reporter systems (Su et al., 2010; Rerksiri et al., 2013; Jeong and Jung, 2015). Through in silico analysis of CREs, we revealed the presence of common 51 CREs in the promoter regions from the PLANTPAN 2.0 database⁹ (Chow et al., 2016) and MEME tool¹⁰ (Bailey et al., 2006). Selected promoter regions and CREs are summarized in Supplementary Table S4. Of these, we have more interest in five unique CREs: DRECRTCOREAT (RCCGAC), ABREMOTIFIOSRAB16B (AGTACGTGGC), ABADESI2 (GGACGCGTGGC), GARE2OSREP1 (TAACGTA), and ANAERO3CONSENSUS (TCATCAC) (Figure 4 and Supplementary Table S4). DRECRTCOREAT is a core motif of dehydration-responsive element/C-repeat (DRE/CRT) found in the promoters of genes from various species. Previous studies reported that OsDREB1A, AtDREB1A and ZmDREB1A bound to (G/ACCGAC) with the different efficiency by competitive DNA binding assays (Sakuma et al., 2002; Dubouzet et al., 2003; Qin et al., 2004) and OsDREB gene encodes transcription activators that function in drought, salt and cold-responsive gene expression (Dubouzet et al., 2003). However, although the Aloe DREB1 can bind to the DRE, it may also bind to other CREs effectively, which can function in a new cold-induced signal transduction pathway (Wang and He, 2007). It has been known that phytohormones including ABA, auxin, gibberellic acid (GA), salicylic acid (SA) and ethylene are related to the cold responses positively or negatively (Miura and Furumoto, 2013; Verma et al., 2016). Among the ABA-responsive CREs, we found that ABREMOTIFIOSRAB16B and ABADESI2 earlier identified from rice Osrab16B promoter and wheat histone H3 promoter were related to ABA-regulated transcription (Terada et al., 1993; Ono et al., 1996; Busk and Pagès, 1998). In addition, GARE1OSREP1 is involved in Gibberellin-responsive element (GARE) found in rice OsREP-1 promoter (Ogawa et al., 2003; Sutoh and Yamauchi, 2003). ANAERO3CONSENSUS found in promoters of anaerobic genes is involved in the fermentative pathway and related to anaerobic response (Mohanty et al., 2005). In summary, DRECRTCOREAT might be related to cold-preferred expression, and ABREMOTIFIOSRAB16B, ABADESI2 and GARE1OSREP1 might be associated with crosstalk between phytohormones and cold stress-preferred expression. The other CREs not mentioned here might have novel roles in driving cold stress-preferred expression and further experiments will be required to clarify our estimation.

To determine the functions of 502 DEGs up-regulated by cold stress in rice roots, we studied their GO terms within the ‘biological process’ category. In all, 15 terms were highly over-represented in our gene list, with p-values < 0.05 and fold-enrichment values of >2-fold (Figure 5 and Supplementary Table S5). We have also previously reported this (Jung et al., 2008b). The terms included ‘L-phenylalanine catabolic process’ (19.9-fold enrichment), ‘response to water’ (16.2), ‘phenylpropanoid metabolic process’ (15.6), ‘oxylipin biosynthetic process’ (12.9), ‘activation of protein kinase C activity by GPCRP signaling pathway’ (9.7), ‘phospholipid metabolic process’ (8.1), ‘gibberellin metabolic process’ (7.3), ‘response to stress’ (7.1), ‘lipid catabolic process ’ (6.1), ‘protein amino acid dephosphorylation’ (5.3), ‘trehalose biosynthetic process’ (5.2), ‘cytochrome complex assembly’ (4.7), ‘lipid biosynthetic process’ (4.4), ‘regulation of transcription’ (3.2), and ‘protein ubiquitination’ (3.2).

Of these, ‘L-phenylalanine catabolic process’ was the most significantly enriched by cold stress while another critical component in that response was ‘phenylpropanoid metabolic process.’ Transcriptome profile analysis of maize (Zea mays) seedlings in response to cold stress has shown that 31 DEGs for phenylalanine metabolism are induced (Shan et al., 2013). Transcript and metabolic profiling of Arabidopsis thaliana (Charlton et al., 2008) has indicated that phenylpropanoids, along with Lys, Met, Trp, Tyr, Arg, Cys, and the polyamine biosynthetic pathway, are important metabolites that are highly accumulated in response to cold stress. Profiling of maize seedling transcripts by Shan et al. (2013) has also revealed the induction of 54 DEGs for phenylpropanoid metabolism. All of these results suggest that the phenylpropanoid metabolic pathway is activated when various plant species are exposed to cold stress.

Metabolic profiling of Camellia sinensis in response to cold (Wang X.C. et al., 2013) has shown that expression is increased for genes involved in the signal transduction mechanism. Three oxylipin biosynthetic-related genes and two trehalose biosynthetic genes are highly expressed in cold-tolerant Elymus nutans (Fu et al., 2016). Moreover, transcriptomics profiling of Lotus japonicus under cold stress has demonstrated that those conditions lead to the upregulation of the phospholipid metabolic process (Calzadilla et al., 2016).

Transcriptome profiling has presented the upregulation of GA metabolism in cold-stressed ‘Meyer’ zoysiagrass (Wei et al., 2015) and greater than threefold induction of gibberellin 2-beta-dioxygenase genes in cassava, which is also related to responses to abiotic and biotic stimuli (An et al., 2012). All of these reports indicate that the gibberellin metabolic pathway is activated during periods of cold stress.

Genes for ‘lipid catabolic process,’ ‘protein amino acid dephosphorylation,’ ‘cytochrome complex assembly,’ ‘regulation of transcription,’ and ‘protein ubiquitination’ also have important roles in the abiotic-stress response (see, e.g., data in Figure 5). For example, in A. thaliana, several lipid catabolism enzymes in rice (in particular, phospholipids A and D) are activated by low temperatures, as manifested by the heightened accumulation of fatty acids (Wang et al., 2006; Usadel et al., 2008). Serine phosphorylation or dephosphorylation is involved in cold activation signaling of Arabidopsis ICE1, and its Ox in Isatis tinctoria confers cold tolerance (Chinnusamy et al., 2003; Xiang et al., 2013). Campos et al. (2003) have reported that a cold-tolerant genotype of Coffea sp. copes with chilling through an enhanced lipid biosynthetic process. Regulation of transcription is also important for cold tolerance. For example, in Arabidopsis, ICE1 and an R2R3-type MYB control the transcriptional regulation of DREB TFs within the mechanism for cold tolerance (Agarwal et al., 2006; Miura et al., 2007). We also identified ‘Protein ubiquitination’ as another important GO term that is also linked with cold tolerance. For example, Arabidopsis HOS1 mediates the ubiquitination and degradation of ICE1 and negatively regulates the response to cold stress (Dong et al., 2006). In summary, the biological processes that we identified here as being closely associated with the cold-stress response provide novel and informative resources for improving our knowledge about regulatory factors involved in the molecular mechanism(s) that enable plants to cope in a low-temperature environment.

The MapMan program is very effective for visualizing diverse overviews associated with high-throughput transcriptome data (Jung and An, 2012). We uploaded Locus IDs for 502 DEGs for the cold-stress response (Supplementary Table S3) to various overviews installed in that program. Among them, 79 elements were assigned to the ‘RNA’ category, 58 to ‘protein,’ 36 to ‘signaling,’ 25 to ‘miscellaneous function’ (‘misc’), 22 to ‘hormone metabolism,’ 17 to ‘stress,’ 14 to ‘development,’ 13 to ‘transport,’ 10 each to ‘lipid metabolism’ and ‘cell wall,’ 7 to ‘secondary metabolism,’ and a smaller number to other functional groups (Supplementary Table S6). Another 154 genes did not have assigned MapMan terms. In particular, the identification of 17 cold stress-regulated elements supports our proposal that they have potential significance for enhancing tolerance when our candidate genes are expressed.

To investigate the significant metabolic pathways involved in the response to cold stress, we analyzed the Metabolism overview associated with 502 DEGs (Figure 6). Among the 44 elements found there, secondary metabolism included six for phenylpropanoids; nine for lipid metabolism, e.g., phospholipid biosynthesis and lipid degradation; 10 for cell wall metabolism, including cellulose synthase and modification; three for mitochondrial electron transport; seven for major carbohydrate (CHO) metabolism; four for minor CHO metabolism; as well as several others related to this stress, such as amino acid, nitrogen, and nucleotide metabolisms (Figure 6A and Supplementary Table S6). These results implied that a rice plant triggers those metabolic pathways as part of its stress response. Similar to our findings from the GO enrichment analysis, ‘L-phenylalanine catabolic process,’ ‘L-phenylalanine metabolic process,’ and category ‘secondary metabolism’ (including ‘phenylpropanoid metabolism’) were over-represented.

Our Regulation overview of 502 DEGs demonstrated that 73 TFs, 30 genes related to protein modification, 21 associated with protein degradation, and 22 related to hormone metabolism were up-regulated in rice during periods of cold stress (Figure 6B). Of these, the TFs were the most abundant, meaning that they are largely involved in regulating the response and tolerance of rice to such conditions. Therefore, those genes should be considered candidates for further study to regulate the cold-stress response in rice. Accordingly, we found 13 WRKY TFs, 10 MYB and four MYB-related TFs, 10 Apetala2/Ethylene Responsive Element Binding Proteins (AP2/EREBPs), five Basic Helix-Loop-Helix (bHLH) genes, five Constans (CO)-like zinc finger family TFs, five C2H2 zinc finger family TFs, and other TFs for this response (Figure 6C and Supplementary Table S6).

In plants, the WRKY TFs have been more actively studied than others, and most of them have positive roles in the cold-stress response in various plant species, including Ipomoea batatas, where the function of a WRKY TF was first described (Ishiguro and Nakamura, 1994). This TF contains a WYRKY domain and a zinc-finger motif. Marè et al. (2004) have reported the role of Hv-WRKY38 in the cold-stress response by Hordeum vulgare, and Ox of WYRKY76 and WYRKY71 has been shown to increase cold tolerance in rice (Yokotani et al., 2013; Kim et al., 2016). Likewise, Ox of CsWYRKY46 in Cucumis sativus regulates tolerance to chilling and freezing (Zhang et al., 2016), and the cold-inducible BcWYRKY46 from Brassica campestris enhances cold tolerance in transgenic tobacco (Nicotiana tabacum) (Wang et al., 2012). In contrast, OsWYRKY45 and OsWRKY13 negatively regulate cold tolerance in rice (Qiu et al., 2008; Tao et al., 2011), while WYRKY34 mediates the cold sensitivity of mature pollen in A. thaliana (Zou et al., 2010) CsWRKY2, a novel WRKY gene from Camellia sinensis, is involved in cold stress responses (Wang Y. et al., 2016).

Like WRKY TFs, MYB TFs have important roles in cold tolerance. They include OsMYB4 OsMYB2 and MYBS3 in rice (Vannini et al., 2004; Su et al., 2010; Yang et al., 2012), MYB15 and HOS10 in Arabidopsis (Zhu et al., 2005; Agarwal et al., 2006), and GmMYBj1 in soybean (Su et al., 2014); and TaMYB3R1 in Triticum aestivum (Cai et al., 2015). Whereas all of those TFs have positive effects, MYBC1 in Arabidopsis negatively regulates cold tolerance (Zhai et al., 2010).

The AP2/EREBP TFs also enhance cold tolerance. They include JcDREB, JcCBF2, BnaERF-B3-hy15, DEAR1, ZmDREB1A, OsDREB1D, and ZmDBP4 analyzed in Arabidopsis (Qin et al., 2004; Tsutsui et al., 2009; Zhang et al., 2009; Wang et al., 2010, 2014; Tang et al., 2011; Xiong et al., 2013); and JERF1, OsDREB1, and AtDREB1A in tobacco (Kasuga et al., 2004; Li et al., 2005; Wu et al., 2007).

A major TF family of other TFs involved in cold tolerance is bHLH. ICE1, ICE2, VabHLH1, and OrbHLH001 analyzed in Arabidopsis (Chinnusamy et al., 2003; Fursova et al., 2009; Li et al., 2010; Xu et al., 2014) and OsbHLH1 in rice (Wang et al., 2003) are involved in cold tolerance. Next, HOS1, a member of the CO-like zinc finger family, regulates cold tolerance in Arabidopsis via CONSTANS degradation (Jung et al., 2012), while OsZFP245, a member of the C2H2 zinc finger family, confers such tolerance in rice (Huang et al., 2009).

Related to protein degradation, signal transduction, and hormone metabolism, a few studies have been conducted. Therefore, future analyses of uncharacterized TFs and the regulatory elements associated with protein degradation, signal transduction, and hormone metabolism identified in this study might shed the light on the effective methods for improving cold tolerance in rice.

To evaluate the significance of our candidate genes, we searched the literature to determine if functions for them have been reported previously. This was accomplished with the online OGRO database, which provides a thorough summary of rice genes that have been characterized through molecular and genetic techniques (Yamamoto et al., 2012). That summary presents the roles of 49 genes according to three agronomic trait categories: morphological, physiological, and resistance/tolerance. The functional identification of genes related to resistance/tolerance traits is the most abundant, with 27 genes being part of that category, including 12 genes involved in cold tolerance; 16, drought tolerance; 11, salinity tolerance; six, blast resistance; five, bacterial blight resistance; two, soil stress tolerance; one each for sheath blight resistance and insect resistance; and four for other stress resistances (Figure 7). Of these, 17 genes are partially responsible for at least two traits in that resistance/tolerance category. OsMAPK5 and OsWRKY45 are involved in tolerance to both biotic stress (bacterial blight and blast) and abiotic stress (drought, salinity, and cold). Others include OsMYB2, ZFP182, OsDREB1A, OsDREB1B, and OsDREB1C, for responses to drought, salinity, and cold; OsbZIP52/RISBZ5 and OsCAF1B, cold and drought; OsTPP1, cold and salinity; and OsCPK4, OsCDPK7, and OsNAC045, drought and salinity. The results from our transcriptome analysis had also suggested that these last three are active in the cold-stress response. We found it interesting that genes induced by low temperatures also function in other abiotic-stress responses. This implies that regulation of those responses is very complex and that intensive crosstalk might occur among them.

Regarding morphological traits, 13 genes are related to dwarfism, five to rooting, four to culms/leaves, three to seeds, three to shoots/seedlings, two to panicles/flowers, and three to other plant components (Figure 7). These results indicate that the cold stress-responsive genes studied here might also affect various traits, e.g., dwarfism, that can inhibit or delay normal growth. Regarding physiological traits, we found that two genes each are related to flowering, germination dormancy, and source activity, while one is related to sterility, and one to other traits (Figure 7). Because our findings demonstrate an interaction between cold stress and diverse morphological/physiological traits, we suggest that future studies should screen mutants and focus on their morphological and physiological phenotypes while also screening phenotypes under cold-stress conditions.

Among the cold-inducible genes identified in our study,OsWRKY71 is induced by cold stress(Figure 1F). As we have reported previously(Kim et al., 2016), its Ox leads to cold tolerance(Figure 8). The survival rate is 19% higherfor OsWRKY71-Ox lines than for the WT, and the transgenics also have 30% higher FWs and 60% higher DWs. EstimatingF_v/F_m values is a good way to depict photosynthetic efficiency under cold stress. Our data indicated that, after 96 h of chilling treatment, this efficiency in OsWRKY71-Ox lines decreased from 0.8 to 0.5 while that value in the WT declined from 0.8 to 0.3. Therefore, the Ox lines are 25% more efficient and OsWRKY71 confers cold tolerance.

The response to low temperatures can be divided into four steps: perception of cold stress, signaling cascades for the response, regulation of gene expression, and protection from freezing damage. Our proposed model (Figure 9) is based on published physiological and biochemical aspects as well as reports of functions for genes involved in the relevant signaling and transcriptional pathways.

We theorize that the first reaction by a plant to chilling is to increase membrane rigidity. This is followed by the generation of ROS, then regulation of phosphate homeostasis and activation of calcium receptors and histidine kinases. The ensuing signal transduction cascades are coupled with signal perception. Examples include MAP kinase cascades and a two-component signaling system by histidine kinase. The former is more likely because the cascades of MAP kinase (OsMAPK5/LOC_Os03g17700.1), MAP kinase kinase (MAP2K; OsMKK4/LOC_Os02g54600.1), and MAP kinase kinase kinase (MAP3K, seven members in Figure 9), are stimulated in response to cold stress, making them the most probable candidates for this pathway. Of them, it has been known that OsMAPK5 positively regulates tolerance to cold temperatures and other sources of stress (Xiong and Yang, 2003).

For the latter possibility, the processes might be more complex. In response to cold, plants use Ca²⁺ as a signal. Although we did not yet identify the histidine kinase genes in rice showing significant induction under cold stress, the signal received by a Ca²⁺ channel might bind to a Ca²⁺ sensor, such as calmodulin (CaM), and CaM-like protein might stimulate Ca²⁺/CaM-dependent protein kinases as suggested in Figure 9. Thereafter, gene expression is regulated by TFs through a process that incorporates CBF/DREB-dependent or -independent pathways.

In the case of the CBF/DREB-dependent pathway, a signal from the map kinase cascades is recognized by ICE1, which encodes a bHLH TF that activates the expression of DREB genes in the downstream pathway by directly binding the promoter regions. This results in stimulation of cold stress-responsive genes that are required for altering cellular metabolism. OsbHLH148 or RERJ1 are probable candidate genes, having the same roles as ICE1 in Arabidopsis, i.e., OsbHLH148 is involved in drought tolerance and RERJ1 functions in normal plant growth and development (Seo et al., 2011). OsDREB1A, OsDREB1B, and OsDREB1C have roles in tolerance to cold, drought, and salinity by triggering the expression of target genes (Ito et al., 2006). Regarding the CRT/DREB-independent pathway, TFs such as OsWRKY71 (Kim et al., 2016), OsWRKY76 (Yokotani et al., 2013), OsbZIP52 (Liu et al., 2012), ZFP182 (Huang et al., 2012), and OsMYB2 (Yang et al., 2012) are components of the trait for cold stress response. For example, a rice line that over-expressor of OsbZIP52 displays a cold-sensitive phenotype (Liu et al., 2012) and the application of such stress induces the expression of OsbZIP52, which then negatively affects the extent of that tolerance.

Although the functions of most genes for cold tolerance have not yet been defined, other types of TFs identified in our meta-expression and MapMan analyses might also be important for regulating tolerance, as indicated by the TF overview presented by MapMan (Figure 6 and Supplementary Table S6).

Among other processes, HOS1, encoding the ring type E3 ligase, participates in the degradation process of ICE1 that is stimulated at low temperature, resulting in inactivation of the CRT/DREB-dependent transcription regulation pathway (Chinnusamy et al., 2003; Dong et al., 2006). Likewise, OST1, encoding the well-known Ser /Thr protein kinase, is activated in response to cold and phosphorylates ICE1, leading to its stability and transcriptional activity (Ding et al., 2015). However, OST1 also hinders the interaction between HOS1 and ICE1, subsequently leading to the degradation of ICE1 under cold stress when HOS1 is suppressed. OsCAF1B, with RNase D activity, functions in post-transcriptional regulation and may affect various pathways for cold tolerance (Chou et al., 2014). OsTPP1 has a role in resistance to abiotic stress. At low temperatures, it also positively regulates the expression of tolerance genes by participating in the glucose deprivation signaling pathway (Ge et al., 2008). Despite these numerous reports, however, all of these hypotheses must still be verified through further experiments.

Cold stress is one of the main environmental factors that adversely affect plant growth and yield. Thus, it is important that we understand this stress signaling and its regulatory network if we are to develop cultivars with greater tolerance. To this end, we have produced a hypothetical model that considers our current findings as well as data derived from earlier research.