Systematic Analysis of Cold Stress Response and Diurnal Rhythm Using Transcriptome Data in Rice Reveals the Molecular Networks Related to Various Biological Processes

Rice (Oryza sativa L.), a staple crop plant that is a major source of calories for approximately 50% of the human population, exhibits various physiological responses against temperature stress. These responses are known mechanisms of flexible adaptation through crosstalk with the intrinsic circadian clock. However, the molecular regulatory network underlining this crosstalk remains poorly understood. Therefore, we performed systematic transcriptome data analyses to identify the genes involved in both cold stress responses and diurnal rhythmic patterns. Here, we first identified cold-regulated genes and then identified diurnal rhythmic genes from those (119 cold-upregulated and 346 cold-downregulated genes). We defined cold-responsive diurnal rhythmic genes as CD genes. We further analyzed the functional features of these CD genes through Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses and performed a literature search to identify functionally characterized CD genes. Subsequently, we found that light-harvesting complex proteins involved in photosynthesis strongly associate with the crosstalk. Furthermore, we constructed a protein–protein interaction network encompassing four hub genes and analyzed the roles of the Stay-Green (SGR) gene in regulating crosstalk with sgr mutants. We predict that these findings will provide new insights in understanding the environmental stress response of crop plants against climate change.


Introduction
Rice (Oryza sativa L.) is a model crop plant and a staple crop across the globe, particularly Asia. Owing to dramatic global climate change and population growth, it has become necessary to increase the production of staple food sources [1,2]. Because rice is a cold-sensitive crop, cold stress encompassing chilling (0-20 • C) and freezing temperature (below 0 • C) is considered a major factor that leads to drastic reduction in rice productivity [3][4][5]. Therefore, various studies on enhancing cold stress tolerance have been reported, with the identification of important factors affecting cold stress response in rice [6]. However, there is evidence that the cold stress response is associated with various biological processes, including biological clocks, disease resistance, metabolism, hormonal balance, and aging [7][8][9]. This complexity is not only a major obstacle to increasing cold stress tolerance but also to the regional adaptation of resistant varieties.
The circadian clock caused the diurnal rhythmic expression of the genes to cope with the earth's rotation and daily repeated environmental fluctuations. Organisms have evolved adaptive mechanisms associated with the light/dark environmental alternations and temperature changes, which manifests as a profound change in metabolism, physiology, and behavior occurring between the day and night cycles in most organisms [10,11]. Circadian clock greatly affect plants, particularly regarding synchronizing biological processes and increasing the efficiency of photosynthesis owing to the non-motile feature of these organisms [12][13][14]. It is involved in various physiological regulatory roles, including flowering time, by controlling the expression and photostability of CONSTANS proteins, which induce FLOWERING LOCUS proteins for flowering [15].
With the recent introduction of high-throughput technology such as microarrays, several reports provided insights into the molecular components of the crosstalk between the circadian clock and cold stress response [16]. There are reports on Arabidopsis that the disruption of the circadian clock changes cold-response-associated transcriptomes, and the knockout of the three pseudo-response regulator genes (PRR), PRR5, PRR7, and PRR9, affects the expression of stress response genes [17,18]. Several cold-related and circadian rhythm genes have been reported in rice. A serine-threonine protein kinase without lysine kinase, OsWNK1, exhibited differential expression under various abiotic stressors, including cold, and also exhibited circadian rhythm [19]. Aquaporin in rice roots is responsible for cold stress-induced acclimation and a cold-induced MYB transcription factor, CMYB1, exhibited a circadian rhythm expression pattern in rice leaves [20,21]. However, the crosstalk between cold stress response and the circadian clock in rice remains unclear from a systems biology approach.
To understand the interplay between the circadian clock and cold stress response in rice, we performed a systemic analysis using two large transcriptome datasets related to the diurnal rhythm and abiotic stress such as drought, salt, cold, and submergence. We identified cold-regulated genes and then identified cold-responsive diurnal rhythmic genes (hereafter CD genes; 119 and 346 genes for up/down-regulation, respectively). Next, we analyzed the functional features of these CD genes through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, after a literature search for functionally characterized genes. Furthermore, we constructed a protein-protein interaction (PPI) network and validated it using the stay-green (sgr) mutant for the regulation of the crosstalk between the cold-response and circadian clock. Based on these results, we proposed that a hypothetical molecular network mediates interplay among various biological processes, including the circadian clock, cold stress response, hormone signaling, and senescence. Our systems biology approach regarding the cold stress response and circadian clock identified some details of the molecular network which could be used to improve the productivity and regional adaptation of rice.

Genome-wide Identification of Cold Stress Response Genes Exhibiting Diurnal Rhythm Expression Patterns Using Meta-Expression Datasets
For systematic analyses aimed at revealing the unknown molecular mechanisms of cold stress response and circadian rhythm, we identified candidate genes responsive to cold stress and diurnal rhythm expression patterns (hereafter cold-responsive diurnal rhythmic genes; CD genes). We identified rice genes unique to cold stress response using a meta-expression dataset consisting of four abiotic stress responses such as drought, salt, cold, and submergence. Consequently, we identified 885 and 572 cold-response genes that were up/downregulated at least 2-fold, respectively (Table S1). Next, we then analyzed the diurnal expression patterns of these genes using Agilent 44K array data of rice leaves under natural field conditions, encompassing the vegetative, reproductive, and ripening stages, at 2-h intervals for 2 days in each developmental stage [22]. Finally, we identified 119 and 346 up/downregulated CD genes, respectively, exhibiting diurnal rhythm expression patterns ( Figure 1, Table S2). These genes were used for further analyses.
In addition to GO enrichment, we performed KEGG enrichment analysis to identify the metabolic pathways involved in cold response and the circadian clock mechanism. Interestingly, only one pathway, photosynthesis-antenna proteins, was enriched, which was consistent with the enriched GO term in downregulated CD genes ( Figure 2B). Using the KEGG mapper software, we identified three light-harvesting proteins involved in this pathway ( Figure 2C) [51]. These results were consistent regarding photosynthesis and the light harvesting process mainly associated with the cold stress response and circadian clock, which was also confirmed with the MapMan software (3.6.0 RC1; Figure S1A, Table S4). Along with the photosynthesis process, GO and MapMan results also implied a variety of potential biological processes involved in this crosstalk ( Figure S1B). Among the upregulated CD genes, seven were present in the network nodes, two of which were functionally characterized. Brassinozole-resistant1 (OsBZR1) is a key transcription regulator in rice with activity regulated by 14-3-3 interactions. The other gene is a mitogen-activated protein kinase in rice (OsMAPK5) that was reported as an ABA-induced regulator that can modulate disease resistance and abiotic stress tolerance [26,30]. Among the downregulated CD genes, two known genes were identified in the network nodes: the waxy gene in rice (Wx), which is associated with the amylose content in rice endosperms, and the rice histone deacetylase gene (OsHDA702), which is associated with plant growth, particularly in the root [44,46]. Among the interactors, 42 known genes were identified, five of which were related to flowering time, including phytochromes A, B, and C (OsPhyA, OsPhyB, OsPhyC); four were related to biotic stress resistance, including calcium-dependent protein kinase 10 (CPK10), somatic embryogenesis receptor-like kinase 1 (SERK1), WRKY62, and disease resistance-responsive gene 8 (DR8); four were related to dwarfism, including decreased DNA methylation (DDM1a), dwarf 1 (D1), pyruvate dehydrogenase kinase 1 (PDK1), and HDA704; and three were related to chloroplast development, including virescent3 (V3), thioredoxin m (Trxm), and spo0B-associated GTP-binding protein (ObgC). Flowering time control gene in rice (rFCA) and pseudo-response regulator (OsPRR37) in rice were also identified. Other identified genes are listed in Table 2.

Construction of Protein-protein Interactions of CD Genes Reveals the Molecular Network Interplay with Various Biological Processes
To further explore the molecular mechanisms regulating crosstalk between the circadian clock and cold stress response, we constructed a hypothetical protein-protein interaction (PPI) network of up/downregulated CD genes using a rice interaction viewer (http://bar.utoronto.ca/interactions/ cgi-bin/rice_interactions_viewer.cgi) [52]. This viewer generated a network of CD genes comprising 764 nodes and 1340 edges ( Figure S2). However, owing to its complexity, we reconstructed the network to only show the interactions between genes for which some information was available. The reconstructed network comprised 208 nodes and 248 edges (Figure 3, Table S5), including both up-and downregulated CD genes.
Among the upregulated CD genes, seven were present in the network nodes, two of which were functionally characterized. Brassinozole-resistant1 (OsBZR1) is a key transcription regulator in rice with activity regulated by 14-3-3 interactions. The other gene is a mitogen-activated protein kinase in rice (OsMAPK5) that was reported as an ABA-induced regulator that can modulate disease resistance and abiotic stress tolerance [26,30]. Among the downregulated CD genes, two known genes were identified in the network nodes: the waxy gene in rice (Wx), which is associated with the amylose content in rice endosperms, and the rice histone deacetylase gene (OsHDA702), which is associated with plant growth, particularly in the root [44,46]. Among the interactors, 42 known genes were identified, five of which were related to flowering time, including phytochromes A, B, and C (OsPhyA, OsPhyB, OsPhyC); four were related to biotic stress resistance, including calcium-dependent protein kinase 10 (CPK10), somatic embryogenesis receptor-like kinase 1 (SERK1), WRKY62, and disease resistance-responsive gene 8 (DR8); four were related to dwarfism, including decreased DNA methylation (DDM1a), dwarf 1 (D1), pyruvate dehydrogenase kinase 1 (PDK1), and HDA704; and three were related to chloroplast development, including virescent3 (V3), thioredoxin m (Trxm), and spo0B-associated GTP-binding protein (ObgC). Flowering time control gene in rice (rFCA) and pseudo-response regulator (OsPRR37) in rice were also identified. Other identified genes are listed in Table 2.

Validation of the Cold Stress Response and Circadian Clock Network with a Case Study Using the Staygreen (sgr) Mutant
Among the functionally characterized genes exhibiting crosstalk between cold stress response and circadian clock, sgr mutant maintains greenness during leaf senescence [37]. Under cold stress for 4 days and recovery for 4 days, we found that the sgr mutant exhibited more cold stress tolerant phenotype than the control plant, Dongjin rice (DJ) ( Figure 4A). Moreover, we confirmed that this tolerant phenotype is not originated from the developmental effect using reproducible experiment with Hwacheong-wx rice in the same genetic background with sgr mutant ( Figure S3). Upregulation of the SGR gene under cold treatment and its diurnal expression patterns were confirmed by qRT-PCR ( Figure 4B). The expression patterns of the four hub genes in our PPI network, i.e., OsMAPK5 (LOC_Os03g17700), OsSnRK1a (LOC_Os05g45420), OsPhyB (LOC_Os03g19590), and OsHDA702 (LOC_Os06g38470) ( Figure 4C) were then examined. All four genes were upregulated in the sgr mutants compared with DJ after cold treatment for 4 days. Two genes were also significantly changed. Consistent with the results of our PPI network analysis, sgr mutants also exhibited altered expression levels of four hub genes in response to cold stress. These results suggest that several Figure 3. Hypothetical protein-protein interaction (PPI) network associated with CD genes. The PPI network constructed from the Rice Interactions Viewer by querying all CD genes. The network was simplified by removing all genes that were not functionally characterized. The full network is shown in Figure S2.

Validation of the Cold Stress Response and Circadian Clock Network with a Case Study Using the Stay-green (sgr) Mutant
Among the functionally characterized genes exhibiting crosstalk between cold stress response and circadian clock, sgr mutant maintains greenness during leaf senescence [37]. Under cold stress for 4 days and recovery for 4 days, we found that the sgr mutant exhibited more cold stress tolerant phenotype than the control plant, Dongjin rice (DJ) ( Figure 4A). Moreover, we confirmed that this tolerant phenotype is not originated from the developmental effect using reproducible experiment with Hwacheong-wx rice in the same genetic background with sgr mutant ( Figure S3). Upregulation of the SGR gene under cold treatment and its diurnal expression patterns were confirmed by qRT-PCR ( Figure 4B). The expression patterns of the four hub genes in our PPI network, i.e., OsMAPK5 (LOC_Os03g17700), OsSnRK1a (LOC_Os05g45420), OsPhyB (LOC_Os03g19590), and OsHDA702 (LOC_Os06g38470) ( Figure 4C) were then examined. All four genes were upregulated in the sgr mutants compared with DJ after cold treatment for 4 days. Two genes were also significantly changed. Consistent with the results of our PPI network analysis, sgr mutants also exhibited altered expression levels of four hub genes in response to cold stress. These results suggest that several biological processes, including cold response and senescence, are linked through the hub genes predicted by our PPI network model.

Discussion
Cold threatens the normal growth and yield requirements of various major crops [92]. The circadian clock rhythmically controls the behavior and physiological activities of plants for improved environmental adaptation [16,93]. The interplays between the circadian clock, light-quality, and coldresponse have been studied intensively in the model plant, Arabidopsis thaliana [13,94]. Circadian clock components such as CIRCADIAN CLOCK ASSOCIATED 1, TIMING OF CAB EXPRESSION1, LATE ELONGATED HYPOCOTYL and Pseudo-Response Regulators (PRR) genes mutually interact with temperature signals [95]. Moreover, the photoreceptors: ZEITLUPE, phytochromes, and cryptochromes which are regulated by GIGANTEA protein, integrate light-quality information with circadian clock [96][97][98][99]. These clock components have a role in regulating the expression of C-repeat

Discussion
Cold threatens the normal growth and yield requirements of various major crops [92]. The circadian clock rhythmically controls the behavior and physiological activities of plants for improved environmental adaptation [16,93]. The interplays between the circadian clock, light-quality, and cold-response have been studied intensively in the model plant, Arabidopsis thaliana [13,94]. Circadian clock components such as CIRCADIAN CLOCK ASSOCIATED 1, TIMING OF CAB EXPRESSION1, LATE ELONGATED HYPOCOTYL and Pseudo-Response Regulators (PRR) genes mutually interact with temperature signals [95]. Moreover, the photoreceptors: ZEITLUPE, phytochromes, and cryptochromes which are regulated by GIGANTEA protein, integrate light-quality information with circadian clock [96][97][98][99]. These clock components have a role in regulating the expression of C-repeat binding factors that bind to the promoter of the cold-regulated genes [100]. In rice, the circadian clock genes have been reported using ortholog search, but it has been mainly studied in depth related to the agronomic traits especially, flowering time [101]. The research of the trait along with the geographical distribution of the rice suggested the relationship between the circadian clock and cold stress [67], but, still, there have been limited reports to shape the molecular network as in the arabidopsis.
Therefore, it is of great developmental significance to explore the possible correlations between cold stress response and circadian clock through systematic analysis in rice.
We first identified 465 CD genes, including 119 upregulated and 346 downregulated genes. These were identified through global transcriptome analyses using public resources, which provided a more extensive scope of utilities that could possibly enhance cold tolerance in rice. These genes will serve as the basis for a series of studies to follow.
To understand the meaningful biological processes underlying the crosstalk between cold stress response and circadian clock in rice, we performed GO and KEGG enrichment analyses. Based on the visualized results, we found that photosynthesis and light harvesting are key biological processes in rice associated with cold stress response. Previous studies also suggested the importance of photosynthesis in cold stress response and circadian clock. For instance, several genes of Zea mays cold-induced genes associated with photosynthesis during long-term cold treatment in maize seedlings [102]. Transgenic rice plants overexpressing choline oxidase A exhibited expanded Fv/Fm when exposed to cold stress [103]. Besides, some researchers indicated that photosynthesis is coupled with diurnal rhythms through light harvesting and CO2 fixation rates by ribulose 1,5-bisphosphate carboxylase/oxygenase (Rubisco) [104]. However, research on photosynthesis and the crosstalk effects of cold stress and circadian clock remains limited. We also found that the GO term 'polysaccharide catabolic process' was highly enriched in the upregulated CD genes. This could be explained by the energy supply required for plants under cold stress.
Our network analysis suggests a hypothetical molecular mechanism underlying the diverse phenomenon related to crosstalk between cold stress response and circadian clock regulation. In this network, OsPhyA, OsPhyB, OsPhyC, and OsPRR37 are key components associated with the regulation of circadian clock. These components interact with both OsMAPK5, a component of Ca 2+ signaling and modulator of disease resistance that is upregulated in response to cold and exhibits diurnal rhythmic expression, and OsHDA702, a component of DNA remodeling that is downregulated in response to cold and exhibits rhythmic expression [44,105]. OsSnRK1a has roles in regulating various metabolism-related processes in plants and is also related to ABA signaling [106,107]. Moreover, D1/RGA1 is related to these four hub genes. Studies of daikoku 1 (d1) mutant using microarray analysis revealed that D1/RGA1 has a role in regulating multiple abiotic stresses, including drought, cold, heat, and salinity [108]. Angel et al. (2016) also revealed decreased drought sensitivity in the erect leaves of d1 mutant owing to the increased photo-avoidance and more effective light harvesting [109,110]. More interestingly, sgr mutant with cold tolerant phenotype is exhibited as erect leaves which is similar to d1 mutant, and the expression level of four hub genes are all upregulated in sgr mutant relative to wild type. Thus, the identification of these components implies that cold-induced SGR and four hub genes (OsMAPK5, OsSnRK1a, OsHDA702, and OsPhyB) might coregulate the physiological response to cold in a circadian clock-dependent manner.
Overall, we hypothesize that the interconnection between different physiological responses, such as cold, and circadian clock is achieved by manipulating the expression of the hub genes (OsMAPK5, OsSnRK1a, OsHDA702, and OsPhyB). Cumulatively, we constructed a framework diagram of multiple elements that are interconnected through the hub genes ( Figure 5). More experimental evidence is required to validate this hypothesis; however, we anticipate that our data will serve as a bridge for integrating initial studies on cold stress and circadian clock.

Collection of Microarray Data and Meta-expression Analysis
To understand the crosstalk between cold stress and circadian clock regulation, we used two different meta-expression datasets: Affymetrix array gene expression dataset for abiotic stresses (E-MEXP-2401, GSE16108, GSE18930, GSE21651, GSE24048, GSE25176, GSE26280, GSE33204, GSE37940, GSE38023, and GSE6901) and a diurnal dataset comprising Agilent 44K array gene expression data from 202 leaf samples collected from nine developmental stages (GSE36040). As previously described [111,112], public transcriptome data were downloaded from the GEO database [113]. After collecting the expression datasets, data were normalized with the Affy and Limma packages using the R programming language [114,115]. After normalization, intensity values were converted to log2 scale, and the log2 fold-change values in response to cold stress (compared with their values in an untreated control) were calculated. Cold responsive genes were defined as those with a  1.5 or −1.5 log2 fold change, and a p-value of  0.05. For the selected genes, log2 fold-change data were visualized with those in other abiotic stress treatments using the MeV software [116]. KMC clustering in MeV was used to identify cold-induced or cold-repressed diurnal rhythmic genes.

Literature Search for Functionally Characterized Genes
The OGRO database was searched to determine the functionally characterized cold-induced diurnal rhythmic genes [23]. Detailed information is presented in Table 1.

Collection of Microarray Data and Meta-expression Analysis
To understand the crosstalk between cold stress and circadian clock regulation, we used two different meta-expression datasets: Affymetrix array gene expression dataset for abiotic stresses (E-MEXP-2401, GSE16108, GSE18930, GSE21651, GSE24048, GSE25176, GSE26280, GSE33204, GSE37940, GSE38023, and GSE6901) and a diurnal dataset comprising Agilent 44K array gene expression data from 202 leaf samples collected from nine developmental stages (GSE36040). As previously described [111,112], public transcriptome data were downloaded from the GEO database [113]. After collecting the expression datasets, data were normalized with the Affy and Limma packages using the R programming language [114,115]. After normalization, intensity values were converted to log2 scale, and the log2 fold-change values in response to cold stress (compared with their values in an untreated control) were calculated. Cold responsive genes were defined as those with a ≥ 1.5 or ≤−1.5 log2 fold change, and a p-value of ≤ 0.05. For the selected genes, log2 fold-change data were visualized with those in other abiotic stress treatments using the MeV software [116]. KMC clustering in MeV was used to identify cold-induced or cold-repressed diurnal rhythmic genes.

Literature Search for Functionally Characterized Genes
The OGRO database was searched to determine the functionally characterized cold-induced diurnal rhythmic genes [23]. Detailed information is presented in Table 1.

GO and KEGG Enrichment Analyses
The rice oligonucleotide array database was used to retrieve GO information for our candidate genes [117]. Fold-enrichment values were obtained by dividing the query number by the query expected value. Significant GO terms were defined as those with a fold-enrichment value of >2 and a hypergeometric p-value of <0.05. The ClusterProfiler software package was then used to perform KEGG enrichment analysis [118]. To use the enrichKEGG function in this package, we used input data containing clustering information and the ID of the rice annotation project database. Additionally, we selected an adjusted cutoff p-value of <0.05, thereby selecting the organism code data and filtering the results. To visualize these results, we used the dotplot function of R studio (version: 1.2.5042) and modified the graph using the ggplot2 software package (version: 3.3.0). Illustrator software was used to refine the presentation.

MapMan Analysis
We used the MapMan toolkit (3.6.0 RC1) to functionally classify our candidate genes [119]. We assigned different colors for diurnal rhythmic genes that were up/downregulated (up, red; down, green) in response to cold stress and visualized all candidate genes simultaneously. In our analysis, we used metabolism and regulation overviews. Detailed information regarding the MapMan analysis is presented in Table S4.

Protein-Protein Interaction Network Construction
A hypothetical protein-protein interaction network was generated to investigate crosstalk between the cold stress response and circadian rhythm. A network file of the queried genes (up/down-CD genes) was obtained from the rice interaction viewer (http://bar.utoronto.ca/interactions/cgi-bin/rice_ interactions_viewer.cgi). Next, a base file containing interaction information was constructed. A data file containing additional information on network nodes was also generated. These data were visualized using the Cytoscape software (3.3.0) [120,121]. We used a p-value threshold of <0.05 and enabled the "fusion option", which merges nodes with similar information to simplify our network.

Plant Material and Stress Treatment
Japonica rice (Oryza sativa L.) (cv. Dongjin and Hwacheong-wx) seeds and sgr defective mutants (cv. Hwacheong-wx) were independently germinated on Murashige and Skoog (MS) medium for one week in an incubator at 28 • C/22 • C (day/night). To simulate cold stress, we exposed 7-day-old seedlings to 4 ± 1 • C for 4 days in a fridge and allowed them to recover for 4 days. For RNA extraction, we collected the leaves of Dongjin and sgr at two time points: 7-day-old seedlings before cold treatment and after 4 days of cold treatment. For confirming diurnal rhythm, we germinated and incubated seeds of japonica rice, cv. Dongjin, on MS medium for a week and transplant the seedlings to a growth chamber (28 • C/25 • C day/night, 14/10 hrs light/dark, and 80% humidity). We then sampled the leaves every two hours in the growth chamber for a day.

RNA Extraction and Quantitative RT-PCR (qRT-PCR) Analysis
Leaf samples from Dongjin and sgr before and after cold treatment (four days) were immediately frozen in liquid nitrogen. After total RNA was isolated using the RNAiso (Takara Bio, Shiga, Japan), first-strand cDNA was synthesized using the SuPrimeScript RT Premix (with oligo (dT), 2×) (GeNet Bio, Daegu, Korea). The synthesized cDNAs were amplified using SYBR Green I (GeNet Bio, Korea), and qRT-PCR was performed on the Rotor Gene Q instrument system (Qiagen, Hidden, Germany). To normalize the amplified transcripts, we used a primer pair for rice ubiquitin 5 (OsUbi5/Os01g22490) [122,123]. All the primers for this analysis are summarized in Table S6.