Comparing time-series transcriptomes between chilling-resistant and -susceptible rice reveals potential transcription factors responding to chilling stress

Low temperature is one of the most important environmental factors that inhibits rice growth and grain yield. Transcription factors (TFs) play crucial roles in chilling acclimation by regulating gene expression. However, transcriptional dynamics and key regulators responding to low temperature remain largely unclear in rice. In this study, a transcriptome-based comparative analysis was performed to explore genome-wide gene expression profiles between a chilling-resistant cultivar DC90 and a chilling-susceptible cultivar 9311 at a series of time points under low temperature treatment and recovery condition. A total of 3,590 differentially expressed genes (DEGs) between two cultivars were determined and divided into 12 co-expression modules. Meanwhile, several biological processes participating in the chilling response such as abscisic acid (ABA) responses, water deprivation, protein metabolic processes, and transcription regulator activities were revealed. Through weighted gene co-expression network analysis (WGCNA), 15 hub TFs involved in chilling conditions were identified. Further, we used the gene regulatory network (GRN) to evaluate the top 50 TFs, which might have potential roles responding to chilling stress. Finally, five TFs, including a C-repeat binding factor (OsCBF3), a zinc finger-homeodomain protein (OsZHD8), a tandem zinc finger protein (OsTZF1), carbon starved anther (CSA), and indeterminate gametophyte1 (OsIG1) were identified as crucial candidates responsible for chilling resistance in rice. This study deepens our understanding in the gene regulation networks of chilling stress in rice and offers potential gene resources for breeding climate-resilient crops.


Introduction
Higher plants encounter a wide variety of environmental stresses in their growth and development (Nawaz et al., 2023).The chilling stress (0°C-15°C) is one of main abiotic stresses bringing serious negative impacts on plant development and causing huge economic losses in crops (Ding et al., 2019).Encountering with low temperatures, crops, especially in the tropical and subtropical climatic zones, show irreversible physiological and biochemical perturbations, including leaf falling, incomplete ripening, seedling wilting, and even death (Lukatkin, 2005;Xu et al., 2008;Lukatkin et al., 2012;Sun et al., 2017).
To gain deeper insights into gene networks responding to the chilling stress, plenty of differentially expressed genes (DEGs), especially for transcription factors (TFs), have been identified when treated with chilling stress (Seki et al., 2001).These TFs play key roles in responding to chilling in plants.For example, overexpression of the C-repeat binding factor1 (CBF1) triggered the expression of cold-responsive (COR) genes and, subsequently, the identification of ICE (inducer of CBF expression) function led to the establishment of the ICE-CBF-COR pathway to enhance the chilling resistance in Arabidopsis (Thomashow, 1999;Chinnusamy et al., 2006;Nie et al., 2022).Overexpression each of osmotic stress/ abscisic acid (ABA)-activated protein kinase 6 (OsSAPK6), ideal plant architecture 1 (IPA1), and OsCBF3 forms a chilling-induced OsSAPK6-IPA1-OsCBF3 signal cascade, and enhances rice chilling tolerance (Jia et al., 2022).Other TFs, such as MaNAC1, OsWRKY63, and ZmMYB31, are also involved in chilling tolerance by regulating the expression of downstream genes as well (Li et al., 2019;Zhang et al., 2022;Yin et al., 2024).In addition, plant hormones usually show a close connection with chilling stresses.ABA regulates several biochemical functions including fruit ripening, leaf senescence, and stomatal closure (Kumar et al., 2022).Studies also show that chilling stress impacts the development and stress resistance of plants by influencing the accumulation of endogenous ABA or inducing key molecular hub genes in ABA synthesis and transduction pathways (Guan et al., 2023).In tomato, overexpression of abscisic acid-insensitive 5 (MaABI5-like) increases fatty acid desaturation and enhances cold tolerance (Song et al., 2022).In melon seedlings, abscisic acidresponsive element-binding factor 1 (CmABF1) and CmCBF4 cooperatively regulate putrescine synthesis to help melon seedlings to response to chilling stress (Li et al., 2022).
Rice, which is cultivated across an expanse of 164 million hectares in over 100 countries, serves as a major staple food crop and provides dietary energy for half of the population in the world (Shafi et al., 2023).Since it is a tropical and subtropical crop, the growth of rice is sensitive to the chilling stress, which affects the seedling growth, heading, grain filling, and seed setting of rice.Therefore, understanding gene networks responsive to the chilling stress will provide crucial clues for breeding climate resilient rice breeding.DC90 is a chromosome segment substitution line that is introgressed with the chilling-tolerance locus CTS-12 from the wild rice, while 9311 is a chilling-sensitive cultivar.After 5 days chilling treatment, there was no significant difference between DC90 and 9311, and phenotypic differences were observed to develop gradually between DC90 and 9311 during the recovery period.After 7 days recovery, the seedlings of 9311 were completely wilted, and the seedlings of DC90 were survived after the chilling stress (Cen et al., 2018, Cen et al., 2020).In previous study, DEGs in DC90 and 9311 under 72h chilling treatment (time-course RNA-seq at 0h, 12h, 24h, and 72h) and 24h recovery (at 84h and 96h) have been characterized by Cen et al (Cen et al., 2020).Nevertheless, transcriptional dynamics and crucial TFs responding to low temperature between DC90 and 9311 at various time points remain unclear.Thus, we aimed to conduct weighted gene coexpression network analysis (WGCNA) and gene regulatory network (GRN) analyses and identified potential TFs responsible for chilling resistance in rice.
In this study, we re-analyzed time-course transcriptomes from DC90 and 9311.The expression data were compared at 0h, 12h, 24h, 72h, 84h, and 96h between two rice cultivars through coexpression and regulatory networks.We revealed key genes that had positive and negative regulatory effects on the low temperature and constructed a gene co-expression network for chilling stress in rice.Further, we discovered that crucial genes involved in the ABA synthesis and signal-transduction pathway that were induced by low temperature, suggesting ABA might plays a significant role in chilling signaling cascades.Importantly, several TFs potentially responsible for chilling resistance and their relationship with each other in the regulatory network were identified.Collectively, our findings identified a hub-TFs regulatory network of chilling stress in rice.These data might provide clues for breeding climateresilient crops.

Materials and methods
2.1 Plant materials and data sources DC90 and 9311 are two rice cultivars that are resistant and sensitive to the chilling stress, respectively (Cen et al., 2018, Cen et al., 2020).All rice seedlings (three-leaf stage) were grown in Guangxi University, China, and then exposed to low temperature at 10°C/8°C (day/night) (artificial light was supplemented with day-13h/night-11h circadian rhythm).The seedling leaves were collected at 12h, 24h, and 72h under chilling stress treatments, and at 84h and 96h for recovery stages; 0h was set as control.The collected samples were immediately ground with liquid nitrogen with three biological replicates.The paired-end reads were generated from Illumina 2500 platform and the RNA-seq datasets were deposited in NCBI under the accession number PRJNA601963.

Co-expression network recognition based on the expression level of DEGs
Co-expression network analysis was performed using the R package WGCNA (Langfelder and Horvath, 2008).Briefly, the FPKM matrix of all DEGs was used as input data and these datasets of all time points with three replicates were filtered to remove some genes at low-expression level (FPKM <1) of any sampling replicates.For these resulting genes, we firstly estimated the soft threshold power by the pickSoftThreshold function, which provided the best-weighted coefficient b = 8 for the scalefree network construction.Second, a correlation matrix was constructed and a topological overlap matrix was calculated based on this correlation matrix.Finally, these genes were classified into different modules according to the topological overlap dissimilarity by average hierarchical clustering algorithm.These modules with different colors were then calculated based on a dynamic tree with minimum group size of 30 genes and merging threshold function of 0.25.Hub genes can represent the maximum correlation points with each module and they are usually represented by the kME value.The threshold used to screen hub genes in each module was |kME| > 0.85.The intricate networks were achieved based on the kME value by employing Cytoscape v3.8.2 (Gnanaraj, 2008).

Construction of GRN
According to the PlantTFDB v5.0 (https://planttfdb.gaolab.org/download.php)annotations, 261 of 3,590 genes encoded TFs were identified.These 261 genes were defined as potential regulators based on the machine learning.We used the R packages GENIE3 ( Van and Geurts, 2018) to measure the confidence levels ("weights" ≥ 0) and the top 10,000 inferred regulatory relationships were found from GRN as potential GRN.The top 50 TFs with the heaviest summed weights in the GRN were retained as the candidates.The networks were drawn by Cytoscape v3.8.2 software (Gnanaraj, 2008).

Results
3.1 Identification of DEGs under chilling stress between two rice cultivars DC90 and 9311 are two rice cultivars that are resistant and sensitive to the chilling stress at seedling stage, respectively.To understand their distinguished responses to chilling stress at transcriptional level, samples were collected for time-course transcriptomes (Cen et al., 2020).For the chilling treatment (see Methods), samples were collected at 0h, 12h, 24h, and 72h, while samples were subsequently collected at 84h and 96h for recovery (Cen et al., 2020).To explore the transcriptome dynamics responding to chilling stress in DC90 and 9311, we analyzed the complete gene expression in seedlings at the 6 time points with three biological replicates.Based on principal component analysis, we found that all time points in two rice cultivars could be divided into two main components (Figure 1A).The first principal component (PC1) accounted for 61% of the variance and was explained by samples with chilling stress and samples without treatment or after recovery.The second principal component (PC2) accounts for 16% of the variation and was interpreted by the genetic differences between DC90 and 9311.This result indicates that transcriptomes of the chilling stress do not only reflect the chilling treatment versus the recovery condition but also present a unique transcriptome composition in different genetic backgrounds.
To identify DEGs responding to chilling stresses, pairwise comparisons among the same 6 points of DC90 relative to 9311 were performed using the DESeq2 program (Love et al., 2014).A total of 3,590 genes were obtained showing a clear distinction between different treatment points and the detailed information for those DEGs were listed in Supplementary Table S1.Among these 6 time points, the greatest number of DEGs was found after 72h chilling treatment (1,204 upregulated and 952 downregulated) (Figure 1B).Most of these DEGs were unique at 72h, including 897 unique upregulated genes and 764 unique downregulated genes (Figures 1C, D), indicating that a longer duration of exposure to the chilling stress led to profound changes of rice transcriptome profile.According to the multiple-point comparison, we detected 54 upregulated DEGs and 51 downregulated DEGs that commonly changed at 6 time points between DC90 and 9311 (Supplementary Tables S2, S3).Globally, the number of upregulated DEGs was more than the number of downregulated genes under chilling treatment (0h, 12h, 24h, and 72h), while the number of downregulated DEGs was more than the number of upregulated genes at each time point under the chilling recovery stages (84h and 96h).These results presented a huge difference of DEGs in response to chilling stress between these two rice cultivars.

Co-expression analysis reveals co-regulated modules of chilling responsive genes
In order to investigate key regulators responsive to chilling stress in rice, lowly expressed gene (the FPKM value <1) were filtered, and remaining DEGs were further employed for WGCNA analysis (Langfelder and Horvath, 2008).We constructed the analysis of network topology for soft-thresholding powers from 1 to 40 and identified a relatively balanced scale independence and mean connectivity (Figure 2A).In the current research, b = 8 was selected with scale free (R 2 = 0.9030) to ensure the network topology.Then a gene tree based on hierarchical cluster were established, and a total of 12 co-expression modules with were obtained except the gray module which was contained the genes (n = 524) that were not effectively clustered (Figure 2B).All genes in the same module shared the similar expression patterns and the number of DEGs in these modules ranged from 35 DEGs (the dark orange module) to 1,291 (the magenta module) (Supplementary Table S4).We identified eigengenes (ME, defined as the first principal component of genes in module) as representative patterns of module and evaluated the similarity of each module by correlating their eigengenes (Figure 2C).
According to these eigengenes expression patterns, we found that at 3 time points after chilling stress, expression levels of eigengenes including salmon, green yellow, cyan, and light green modules were showed significantly higher in DC90 when compared to 9311.At the same time, expression levels of DEGs in these module eigengenes were higher under chilling stage (12h, 24h, and 72h) compared to unchilling stage (0h, 84h, and 96h), indicating that they were consistently stimulated in DC90 and might play crucial roles in response to chilling stress (Figures 2D-G).In addition, some modules also correlated with specific genotypes according to expression profiles of eigengenes (Figures 2H, I).The eigengenes of DC90 showed higher expression levels in tan modules but showed lower expression levels in light yellow modules when compared to 9311, suggesting that differences in gene expression between these two modules are mainly due to variety differences, rather than being significantly correlated with cold stress.

Gene ontology enrichment analysis of DEGs in different modules
To assess if these DEGs have any biological function in response to chilling stresses, we used gene ontology functional enrichment analyses for DEGs corresponding to each module and the top 15 GO terms in each module were given (Figure 3A; Supplementary Figure S1).For the six modules in Figures 2D-I, we have discovered many terms potentially involved in chilling stress response.In the salmon module, the most enriched significant processes contained transcription regulator activity, temperature stimulus, cold, hormones, hormone-mediated signaling pathways and cellular responses to hormone stimulus (Figure 3A).These results suggested that plant hormones might play key roles in response to chilling stress in rice.In the green yellow module, the most enriched significant processes were responsive to translation, peptide biosynthetic process, peptide metabolic process, and protein metabolic process, suggesting that genes in this module were likely involved in protein synthesis and metabolism (Figure 3B).In the cyan module, transcription regulator activity was the most significant function-enriched GO term, indicating that the TFs might play crucial roles in response to chilling stress (Figure 3C).In the light green module, the most significant processes were enriched for sequence-specific DNA binding, response to stress, and response to stimulus, while in the tan module, oxidoreductase activity, metal ion binding, and cation binding were enriched (Figures 3D, E).Although the genes in light yellow modules were not significantly enriched in any GO functions, we found that the GO annotations of many genes in this module were related to calmodulinbinding protein, oxidoreductase, and heavy metal-associated (Supplementary Table S4).Taken together, the analysis reveals that chilling stress influences the expression of genes that regulate plant hormone levels, protein contents, and TF activities, and these DEGs may be severed as a versatile resource for breeding chilling-resistant plants, thus we keep salmon, green yellow, cyan, and light green modules as the "rice chilling resistance response" modules for subsequent analyses.

ABA biosynthesis and signaling genes actively respond to the chilling stress
In the GO analysis, we found that ABA was a highly enriched term, which is consistent with previous study (Cen et al., 2020).To understand roles of ABA in rice chilling stresses, the intricate synthesis and signal-transduction pathway of ABA responding to low temperature were further analyzed (Figure 4).In the ABA biosynthesis pathway, 9-cis-epoxycarotenoid dioxygenase (OsNCED3/5) displayed a high expression level in DC90 than that in 9311 after 12h chilling treatment, whereas aldehyde oxidase (OsAAO3;1) showed a low expression level at all chilling hours in two rice cultivars (Figure 4A).Of particular significance, genes situated at the inception of the ABA signal-transduction cascade exhibited obvious induction in the wake of chilling signals.
(Figure 4B).Among them, OsPP2C11, OsPP2C27, and OsPP2C41 showed high expression profiles in 9311 at 72h than that in DC90 that might because of PP2Cs as a negative regulator of ABA signaling pathway.We also investigated expression levels of genes involved in other hormone pathways, such as brassinosteroids (BR), cytokinin (CTK), ethylene (ETH), gibberellin (GA), auxin (IAA), and jasmonic acid (JA), but no significant differences were found between these two varieties under chilling stress (Supplementary Figures S2-S7).These results indicate that under low temperature stress, ABA synthesis and signal transduction pathways play key roles in responding to chilling stress in DC90 versus 9311.

Identification of hub TFs responding to the chilling stress
TFs served as key regulators to regulate plant development and responses toward the chilling stress (Kaufmann et al., 2010).To investigate hub TFs and genes involved in chilling resistance response, we calculated the module membership values (kME) among the four "rice chilling resistance response modules" based on WGCNA packages.We chose candidate TFs with the |kME| >0.85 in these modules, and this analysis retained a lot of hub genes in each module (Supplementary Table S5).Based on WGCNA analysis, totally, we identified 15 hub TFs from four modules that belong to MYB, ethylene-responsive element binding factor (ERF), zinc finger and INDETERMINATE DOMAIN protein (IDD), lateral organ boundaries domain (LBD), and other TF families (Figure 5A).These TFs formed a complex regulatory network, which may play crucial roles in chilling resistance in rice.We also investigated the expression profiles of these 15 TFs in DC90 and 9311 among six different points (Figure 5B).We found 13 of these 15 TFs were highly induced at 12h or 24h after chilling stress in DC90 compared with that in 9311, while the other two TFs were highly accumulated 72h chilling stress in DC90.All these 15 TFs expressed higher under chilling stress (12h, 24h, and 72h) in DC90 than that of 9311, suggesting that these TF might serve as key TFs to regulate the chilling resistance.

Inferring the TFs regulatory network by machine learning
WGCNA is an effective method to explore the gene coexpression networks and hub genes; however, the regulatory relationships between different genes are not easy to evaluate through this method.An analytical pipeline (see Methods) suggested 15 TFs served as reliable candidates for the chilling response network (Figure 5), but the regulation direction between these TFs was not captured.Therefore, we employed the GENIE3 package (Van and Geurts, 2018), which was used to predict the regulatory relationships of different genes and TFs based on a machine learning algorithm.
According to the random forest method, 50 TFs, which rapidly respond to the chilling stress, with the potential regulatory relationships were explored by GENIE3 (Supplementary Table S6).After associating these 50 TFs with 15 hub TFs presented in former modules (salmon, green yellow, cyan, and light green) related to chilling response, we finally retained five TFs, OsCBF3 (LOC_Os09g35030, an ERF family member), OsZHD8 (LOC_Os04g35500, a ZF-HD family member), OsTZF1 ( L O C _ O s 0 5 g 1 0 6 7 0 , a C3 H f a m i l y m e m b e r ) , OsC SA (LOC_Os01g16810, an MYB family member), and OsIG1 (LOC_Os01g66590, an LBD family member), as key regulators controlling gene expression networks responding to chilling stress in rice (Figure 6A; Supplementary Table S7).Among these TFs, OsCBF3 is a well-studied TF in rice chilling stress.In rice, OsSAPK6 interacts with IPA1 and phosphorylates IPA1 to stabilize it, allowing IPA1 to accumulate and activate the expression of OsCBF3 downstream to enhance rice chilling stress resistance (Jia et al., 2022).Thus, this analysis indicates our analysis is efficient and reliable.
Subsequently, we extracted the sub-regulatory network based on these five TFs of the entire gene regulator network (GRN).This subset of data contained 141 genes for salmon, greenyellow, cyan, and lightgreen modules accounting for approximately 59.75% of the total 236 genes of the entire sub-network (Figure 6A).A lot of genes were simultaneously regulated by multiple TFs (Figure 6B).For example, OsCBF3 in cyan module was calculated to be regulated by OsCSA, and the other four TFs (OsTZF1, OsZHD8, OsIG1, and OsCSA) in the salmon module were also predicted to regulate each other, indicating a complex feedback regulation network among the crucial regulators (Figure 6C).
In our analysis, we uncovered that five TFs might play as key regulators that fitness the resistance of chilling stress in DC90 compared to 9311.To validate the accuracy of our results, we used transcriptomes derived from the chilling resistant (Volano) and chilling susceptible (Thaibonnet) varieties treated with chilling stresses (4°C for 0h, 2h, and 10h) conditions (Figure 6D).We found that OsTZF1 and OsCBF3 were highly expressed at 2h, while OsCSA, OsIG1, and OsZHD8 were showed higher expression at 12h in the chilling resistant cultivar Volano than Thaibonnet.Similar result was also detected in the tolerance glutinous rice 89-1 and the susceptible indica rice R07 (Figure 6E).Taken together, we suppose that these five TFs serve as major and common TFs responding to the chilling stress and play crucial roles in rice chilling resistance.

Discussion
Chilling stress often leads to a severe decline for rice growth distribution.Therefore, exploring new rice chilling tolerance regulatory genes/TFs will provide key clues for breeding climateresilient rice.To this end, understanding the transcriptome dynamics under chilling stress in rice is crucial.In our study, we investigated a set of time-course transcriptomes in rice of DC90 (a chilling-resistant cultivar) and 9311 (a chilling-susceptible cultivar) treated with a low-temperature stress.We revealed 12 modules based on the DEGs and a set of hub genes/TFs involved in chilling stress.
Chilling stress can induce alterations of the gene expression profiles in plants.We have found that a lot of DEGs showed rapidly changes at transcriptional levels in response to chilling stimuli in DC90 and 9311.However, it is difficult to explore crucial genes, which might involve in the chilling conditions because of numerous members of DEGs.WGCNA was confirmed as an effective method to help researchers find key hub genes related to specific traits among a large concentration of DEGs (Kuang et al., 2021;Chen et al., 2024;Hu et al., 2024;Xie et al., 2024).In this research, we have divided the DEGs involved in chilling stress into 12 co-expression modules using the WGCNA algorithm to find key genes.Performing functional enrichment annotation on these modules, we found some significant GO terms of "response to temperature stimulus," "response to hormone," "protein metabolic process," and "transcription regulator activity."These results revealed complicated networks in rice seedlings responsible for or adaptation to the chilling stress.
Based on WGCNA, our analysis revealed some "rice chilling resistance response" modules (salmon, green yellow, cyan, and light green modules), and 15 hub TFs from these modules involved in chilling conditions in rice.Among these 15 TFs, some of them have been validated and confirmed as key TFs related to the chilling stress and other abiotic stress, such as OsDREB1C, OsIDD14, OsIDD12, and OsMYB91.OsDREB1C serves as a DREB homolog encoding an SLK protein which has been proved to response to chilling acclimation (Mao and Chen, 2012).IDD14 protein interacts with IDD12 and IDD13 to form a transcriptional complex, activating MAPK expression (Cui et al., 2022).OsMAPK3 and OsMAPK6 positively regulate rice chilling tolerance and OsMAPK3 can effectively respond to ABA stimulation, so we speculate that IDD12/14 might response to low temperature through MAPK cascade signaling pathway (Liu et al., 2023).In addition, overexpression of OsMYB91 enhances tolerance to abiotic stress and alters the sensitivity of seed germination to ABA (Zhu et al., 2015).Thus, our data support that these 15 TFs are potential regulators in chilling stress response.However, it is noted that roles of these TFs in low temperature remain to be further validated.
hormones, such as BR, ETH, JA, and ABA, have been reported to play key roles in the chilling stress (He et al., 1990;Cao et al., 2012;Xiang et al., 2017;Shu et al., 2022).We found ABA was a crucial factor in response to chilling stress in DC90 versus 9311.Previous study showed that ABA response genes were regulated by OsCBF3 to resist low temperature (Jia et al., 2022;Li et al., 2022).The chilling stress dramatically increases endogenous ABA level in rice, which results in a significant upregulation of genes involved in ABA biosynthesis and signaling pathways (Xiang et al., 2017).The silencing of ABI5, which is composed of PYR/PYL/RCAR receptors in the core ABA signaling pathway, significantly reduces cold resistance in both rice and banana (Song et al., 2022).OsNAC5 stabilizes the OsABI5 protein by inhibiting proteasome-mediated degradation to improve the cold resistance of rice (Li et al., 2024).In our data, the ABA pathway genes OsNCED1, OsPP2C2/5/7, and OsbZIP77 showed higher expression levels at 12h under low temperature in DC90 than in 9311.Although OsPYL/RCR10, OsAAO3;1, and OsbZIP46/69 all showed lower expression profiles in these two species under chilling stimulate, they were upregulated at 84h (recovery stage), suggesting that they might be some negative regulatory factors involved in chilling stress.
Although WGCNA is an effective method for analyzing gene co-expression networks, the regulatory relationships between coexpressing genes need to be explored.Here, we used a machine learning algorithm GENIE3, which was based on the inference of a collection of regression trees to predict regulatory networks.In previous studies, many nitrogen deficiency-responsive genes and TFs were predicted to function as key regulators within the network in response to nitrogen deficiency based on GENIE3 in rice (Ueda et al., 2020).In this study, we further used GENIE3 to evaluate the top 50 TFs as the candidates.Combining WGCNA and GENIE3 tow algorithms, five crucial TFs (OsCBF3, OsZHD8, OsTZF1, OsCSA, and OsIG1) associated with the low-temperature resistance in rice were predicted.OsCBF3/DREB1A has been proven to bind with multiple genes and TFs to regulate the cold resistance of downstream genes (Jia et al., 2022).OsZHD8 can bind to the promoter of OsDREB1B, which is specifically induced by low temperature, and shows an increasing expression after ABA treatment in rice (Figueiredo et al., 2012), suggesting that it might have function on the chilling stress.Overexpression of OsTZF1 causes ABA-induced growth inhibition of rice seedlings compared with WT, suggesting OsTZF1 positively regulates ABA response (Zhang et al., 2012).Meanwhile, the high expression level of OsTZF1 enhances resistance to hydrogen peroxide, salt, and drought stresses, suggesting that OsTZF1 might be a multiple abiotic stress resistant regulators and probably play an important role in resisting chilling stress (Jan et al., 2013).The R2R3 MYB TF OsCSA can increase adaptations to adverse environments by balancing glucose and ABA homeostasis (Sun et al., 2021).OsIG1 plays an important role in regulating rice floral organs development (Zhang et al., 2015), and its role in cold stress and ABA signaling pathways is currently unclear.Altogether, our results suggest the combination WGCNA and GENIE3 would be an efficient way to predict hub/key regulators in complicated networks.
In the original article, combining analysis of the differentially accumulated metabolites (DAMs) and DEGs, it has been predicted that galactinol, b-alanine, glutamate, naringenin, serotonin, ABA, and OsNCED3 might be involved in the chilling tolerance variation of 9311 and DC90 (Cen et al., 2020).However, core TFs involved in regulating the biological processes of chilling stress remain unclear.In our study, crucial genes/TFs of low temperature tolerance were identified through WGCNA in combination with GENIE3 network based on DEGs.In addition, we identified a lot of TFs and constructed a potential TF regulatory network.Our results not only reveal common candidates (such as the ABA pathway, and OsNCED3), but also uncover more crucial TFs, which may reflect the application of different prediction methods.These two sets of results could be complemented to each other and will provide more candidate genes for rice breeding and deepen our understanding mechanisms underlying how the rice resists low temperature.

Conclusions
In this study, we performed a transcriptome analysis of rice seedlings in two cultivars under chilling stress at 6 time points.We identified 3,590 DEGs including chilling-responsive genes that were dynamically regulated under low temperature.Specifically, abscisic acid biosynthesis and signaling, coupled with transcription regulator activity, emerges as closely intertwined with the lowtemperature resistance of rice.The integration analysis of hub gene networks has shed light on the key TFs triggered in response to chilling conditions.Through the correlation investigation of WGCNA and GRN, five TF candidates were predicted involved in the processes of chilling response.Our results provide invaluable gene resources for breeding climate resilient rice.
FIGURE 1 Transcriptomic analyses of DC90 and 9311 under 72h chilling treatment and 24h recovery.(A) Principal component analysis (PCA) of DC90 and 9311 samples at 6 time points (0h, 12h, 24h, 72h, 84h, and 96h).(B) Statistics of differentially expressed genes (DEGs) numbers in DC90 relative to 9311 at 6 time points.Each dot presents a DEG.Dots in blue show downregulated genes, and dots in red represent upregulated genes.The left coordinate (y-axis) represents the log 2 fold change values of DEGs and the right coordinate (y-axis) represents the numbers of DEGs.(C) Upset plot shows unique and shared upregulated genes at each time point in DC90 relative to 9311.(D) Upset plot shows unique and shared downregulated genes in DC90 relative to 9311.
FIGURE 2 Weighted gene co-expression analysis (WGCNA) of DEGs between DC90 and 9311 under 72h chilling treatment and 24h recovery.(A) Scale independence and mean connectivity of the network for various soft-thresholding powers.(B) A hierarchical cluster tree of the DEGs distributed in different modules.The assigned modules are depicted by color bands and branches.Genes are represented by the tips of the branches.(C) The correlation between co-expression networks based on eigengenes.The intensity of color on the heatmap signifies the extent of overlap and the darker yellow denoting greater correlation.(D-I) Expression patterns of co-expressed genes in six modules, including (D) salmon, (E) green yellow, (F) cyan, (G) light green, (H) light yellow, and (I) tan modules.The expression quantity was processed with z-score normalization, and the relative expression levels are shown in colorful scale.
FIGURE 3GO enrichments for co-expression genes in five selected modules, including (A) salmon, (B) green yellow, (C) cyan, (D) light green, and (E) tan modules.
FIGURE 5 Identification of hub transcription factors (TFs) in three selected modules involved in chilling stress and their co-expression network.(A) Visualization of the 15 key TFs co-expression networks.Hub TFs of salmon, green yellow, and cyan modules are indicated with salmon, green yellow, and cyan color, respectively.Gray lines indicate the correlations between hub TFs.The size of the points represents the level of |KME|, and the width of the gray line represents the strength of co-expression intensity.(B) The heatmap of these 15 hub TFs in DC90 and 9311.The values in the heatmap represent the z-scores of FPKM at 6 time points in two cultivars.The red and blue colors indicate a high and low expression profiles, respectively.