Transcriptomic Profiling of Cold Stress-Induced Differentially Expressed Genes in Seedling Stage of Indica Rice

Cold stress significantly constrains the growth, development, productivity, and distribution of rice, particularly the indica cultivar, known for its susceptibility to cold, limiting its cultivation to specific regions. This study investigated the genes associated with cold responsiveness in the roots of two indica cultivars, SQSL (cold-tolerant) and XZX45 (cold-susceptible), through transcriptome dynamics analysis during the seedling stage. The analysis identified 8144 and 6427 differentially expressed genes (DEGs) in XZX45 and SQSL, respectively. Among these DEGs, 4672 (G2) were shared by both cultivars, while 3472 DEGs (G1) were specific to XZX45, and 1755 DEGs (G3) were specific to SQSL. Additionally, 572 differentially expressed transcription factors (TFs) from 48 TF families, including WRKY, NAC, bHLH, ERF, bZIP, MYB, C2H2, and GRAS, were identified. Gene Ontology (GO) enrichment analysis revealed significant enrichment of DEGs in the G3 group, particularly in the “response to cold” category, highlighting the crucial role of these specific genes in response to cold stress in SQSL. Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated pronounced enrichment of DEGs in the G3 group in metabolic pathways such as “Pyruvate metabolism”, “Glycolysis/Gluconeogenesis”, and “Starch and sucrose metabolism”, contributing to cold tolerance mechanisms in SQSL. Overall, this study provides comprehensive insights into the molecular mechanisms underlying cold responses in the indica cultivar, informing future genetic improvement strategies to enhance cold tolerance in susceptible indica rice cultivars.


Introduction
Rice (Oryza sativa L.) plays a critical role as a staple crop, providing more than half of the world's population with their primary source of carbohydrates [1]. Its cultivation encompasses diverse ecological conditions, ranging from upland to deep-water environments and spanning various altitudes, climates, and soil types [2]. However, rice production faces numerous challenges due to abiotic stresses, including cold, high temperatures, drought, salinity, and flooding. Cold stress, in particular, poses a significant obstacle to rice farming, leading to chlorosis, reduced seedling growth, stunting, withering, decreased tillering, and ultimately reduced productivity in susceptible cultivars [3]. The global impact of low-temperature stress on rice is substantial, affecting approximately 15 million hectares in 24 countries. In South and Southeast Asia alone, about 7 million hectares of land are unsuitable for rice cultivation due to the detrimental effects of cold stress [4]. Unanticipated cold weather events, resulting from extreme climatic conditions, further compound this challenge [5]. Cold stress has a significant impact on rice yield in countries such as Australia, Korea, China, Japan, and India, where dry-season boro rice is particularly vulnerable during the seedling stage, affecting approximately 4 million hectares of rice-growing areas and causing growth delays. Moreover, these delays coincide with high-temperature stress during the flowering stage, leading to diminished yields [6]. Consequently, the development of high-yielding rice varieties with cold stress tolerance assumes paramount importance in enhancing rice production in these vulnerable regions. Additionally, a comprehensive understanding of the molecular mechanisms underlying chilling stress tolerance is of critical significance.
Plants have evolved diverse mechanisms to counteract stresses, adapting to different stress types, intensities, and ecological contexts [7]. Various stressors such as wounds, pathogen invasion, UV radiation, flooding, drought, salinity, heat, and cold stress, either individually or in combination, trigger an increase in reactive oxygen species (ROS) levels [8].
Maintaining ROS homeostasis assumes a crucial role in regulating the expression of coldresponsive genes in rice [9]. Previous investigations have revealed the involvement of cold stress binding factors (CBFs) in facilitating ROS detoxification [10]. Furthermore, a ROSmediated regulatory module has been identified as an early component of the chilling stress response pathway in rice [11]. Reactive oxygen species are considered a converging point for multiple gene networks in response to stress [8]. Transcription factors (TFs) serve as key regulators, orchestrating the expression of downstream target genes within a network. Prior research has underscored the significance of WRKY, MYB, bHLH, bZIP, AREB/ABF, NAC, and DREB1/CBF in governing the expression of genes responsive to salinity, temperature, and drought stress in rice [12]. Heat shock factor (HSFs) genes have also emerged as crucial nodes for cross-talk in the rice stress response [13]. Recent discoveries have demonstrated that a common set of TFs, coupled with stress signaling networks and upstream regulatory gene regulons, govern plant responses to multiple stresses [10,[14][15][16]. Target TFs operate in a network mode to regulate molecular, biochemical, genetic, and physiological processes in response to various stressors [17]. The binding of TFs in conjunction with cis-elements provides a structural basis for generating distinct patterns of gene expression [18]. In a recent study involving the cloning of the COLD1 quantitative trait locus, a mechanism underlying chilling tolerance was elucidated, whereby the QTL interacts with the alpha subunit of G-protein to activate the Ca 2+ channel. This activation enhances the GTPase activity of G-protein, ultimately leading to the expression of chilling tolerance in japonica rice [19].
The transcriptome encompasses the entirety of transcripts present within a cell, encompassing their types and quantities, during a specific developmental stage or physiological condition [20]. With the advent of next-generation sequencing (NGS), RNA sequencing (RNA-seq) has emerged as an innovative technology for transcriptome analysis, enabling precise identification of differentially expressed genes (DEGs) and potential molecular mechanisms [20][21][22]. RNA-seq primarily relies on NGS platforms that encompass comprehensive systems for library construction, sequencing, and analysis [22]. Over the past decade, RNA-seq has been extensively employed to compare the transcriptomes of coldtolerant and cold-sensitive plants under low-temperature stress conditions [23][24][25][26][27].
In this study, we conducted an investigation into the transcriptomic dynamics of the root tissues in two indica rice cultivars, SQSL and XZX45, under conditions of optimal temperature and cold stress using high-throughput RNA-seq technology. The aim was to gain insights into the intricate molecular responses occurring at the transcriptome level in these cultivars when subjected to cold stress. By employing a comprehensive RNA-seq approach, we sought to elucidate the complex regulatory mechanisms and identify key genes and pathways involved in the indica rice root's response to cold stress.

Phenotyping of Contrasting Genotypes for Cold Stress Tolerance and Recovery
Two cultivars, SQSL and XZX45, were subjected to a controlled low-temperature treatment at 8 • C for a duration of 5 days, followed by a subsequent 5 days recovery treatment. Under the imposed low-temperature stress conditions, the cold-susceptible genotype XZX45 exhibited notable withering in both its belowground and aboveground organs. Conversely, in comparison to XZX45, the cold-tolerant genotype SQSL displayed sustained growth and maintained normal leaf coloration, with its belowground and aboveground components demonstrating substantial recovery, approaching near-optimal levels ( Figure 1A). Following the 5 days low-temperature stress period and subsequent 5 days recovery treatment, SQSL exhibited an impressive average survival rate of 99%, while the survival rate of XZX45 was a mere 8% ( Figure 1B).

Phenotyping of Contrasting Genotypes for Cold Stress Tolerance and Recovery
Two cultivars, SQSL and XZX45, were subjected to a controlled low-temperature treatment at 8 °C for a duration of 5 days, followed by a subsequent 5 days recovery treatment. Under the imposed low-temperature stress conditions, the cold-susceptible genotype XZX45 exhibited notable withering in both its belowground and aboveground organs. Conversely, in comparison to XZX45, the cold-tolerant genotype SQSL displayed sustained growth and maintained normal leaf coloration, with its belowground and aboveground components demonstrating substantial recovery, approaching near-optimal levels ( Figure 1A). Following the 5 days low-temperature stress period and subsequent 5 days recovery treatment, SQSL exhibited an impressive average survival rate of 99%, while the survival rate of XZX45 was a mere 8% ( Figure 1B). (B) Survival rates of SQSL and XZX45 plants after a 5-day recovery period subsequent to a 5-day cold treatment (5 °C). Statistical significance is denoted as follows: ns (not significant) for p ≥ 0.05, and *** for p < 0.001.

Transcriptome Sequencing and Differentially Expressed Genes in Response to Cold Stress
We employed RNA-seq technology to investigate the transcriptome dynamics in the roots of SQSL and XZX45 under optimal temperature and cold stress conditions. A total of 36 samples, including three independent biological replicates at three different time points, namely 1, 3, and 5 days after cold treatment, were analyzed, resulting in 241.2 Gb of clean data (Table S1). The clean reads showed high quality, with average Q20 and Q30 values exceeding 97.55% and 90.59%, respectively (Table S1). After mapping the reads to the rice reference genome of cv. Nipponbare, an average mapped read ratio of 85.36% (ranging from 70.20% to 94.60%) was achieved (Table S1). Read counts mapped to each gene were determined using featureCounts [28], and transcripts per million (TPM) values were calculated. Genes with a TPM ≥ 1 were considered expressed ( Figure S1).
To assess the global differences in transcriptomes between SQSL and XZX45 under cold stress, Pearson correlation coefficients (PCC) were calculated based on TPM values of expressed genes for all samples from the three time points. The PCC values of the three biological replicates of SQSL and XZX45 exceeded 0.9 (Figure 2A,B). Principal component . Statistical significance is denoted as follows: ns (not significant) for p ≥ 0.05, and *** for p < 0.001.

Transcriptome Sequencing and Differentially Expressed Genes in Response to Cold Stress
We employed RNA-seq technology to investigate the transcriptome dynamics in the roots of SQSL and XZX45 under optimal temperature and cold stress conditions. A total of 36 samples, including three independent biological replicates at three different time points, namely 1, 3, and 5 days after cold treatment, were analyzed, resulting in 241.2 Gb of clean data (Table S1). The clean reads showed high quality, with average Q20 and Q30 values exceeding 97.55% and 90.59%, respectively (Table S1). After mapping the reads to the rice reference genome of cv. Nipponbare, an average mapped read ratio of 85.36% (ranging from 70.20% to 94.60%) was achieved (Table S1). Read counts mapped to each gene were determined using featureCounts [28], and transcripts per million (TPM) values were calculated. Genes with a TPM ≥ 1 were considered expressed ( Figure S1).
To assess the global differences in transcriptomes between SQSL and XZX45 under cold stress, Pearson correlation coefficients (PCC) were calculated based on TPM values of expressed genes for all samples from the three time points. The PCC values of the three biological replicates of SQSL and XZX45 exceeded 0.9 (Figure 2A,B). Principal component analysis (PCA) was conducted, indicating significant differences in transcriptomes among different genotypes and treatments ( Figure 2C). Hierarchical cluster analysis further confirmed the distinct clustering of the three biological replicates of SQSL and XZX45, resulting in the division of the 36 samples into 12 groups, highlighting the significant differences in transcriptomes among the experimental groups ( Figure 2D). analysis (PCA) was conducted, indicating significant differences in transcriptomes among different genotypes and treatments ( Figure 2C). Hierarchical cluster analysis further confirmed the distinct clustering of the three biological replicates of SQSL and XZX45, resulting in the division of the 36 samples into 12 groups, highlighting the significant differences in transcriptomes among the experimental groups ( Figure 2D). Differentially expressed genes (DEGs) under cold stress conditions were identified using a threshold of FDR < 0.05 and |log2foldchange| > 1. In SQSL, a total of 16,283 DEGs were identified, with 10,298 (5587 up-and 5341 down-regulated), 11,072 (4989 up-and 6083 down-regulated), and 11,554 (5134 up-and 6420 down-regulated) DEGs found in SQSL1d (SQSL after 1 day of cold treatment), SQSL3d, and SQSL5d, respectively. Meanwhile, 19,040 DEGs were identified in XZX45, with 11,178 (6000 up-and 5178 down-regulated), 14,626 (6875 up-and 7751 down-regulated), and 15,192 (7021 up-and 8171 downregulated) differentially expressed genes in XZX451d, XZX453d, and XZX455d, respectively ( Figure 3).
Comparing the number of DEGs between SQSL and XZX45, it was evident that XZX45 had notably more DEGs than SQSL. While the number of DEGs in SQSL and XZX45 after 1 day of cold treatment did not show a significant difference, substantial differences were observed after 3 and 5 days of cold treatment ( Figure 3B,C). These results indicate that cold stress had a more pronounced and sustained impact on gene expression levels in XZX45 compared to SQSL, suggesting that SQSL exhibited stronger tolerance and adaptability to cold stress than XZX45. A total of 9899 genes were detected in both genotypes across the three time points, with 8144 genes in SQSL and 6427 genes in XZX45 Differentially expressed genes (DEGs) under cold stress conditions were identified using a threshold of FDR < 0.05 and |log 2 foldchange| > 1. In SQSL, a total of 16,283 DEGs were identified, with 10,298 (5587 up-and 5341 down-regulated), 11,072 (4989 up-and 6083 down-regulated), and 11,554 (5134 up-and 6420 down-regulated) DEGs found in SQSL1d (SQSL after 1 day of cold treatment), SQSL3d, and SQSL5d, respectively. Meanwhile, 19,040 DEGs were identified in XZX45, with 11,178 (6000 up-and 5178 downregulated), 14,626 (6875 up-and 7751 down-regulated), and 15,192 (7021 up-and 8171 downregulated) differentially expressed genes in XZX451d, XZX453d, and XZX455d, respectively ( Figure 3).
Comparing the number of DEGs between SQSL and XZX45, it was evident that XZX45 had notably more DEGs than SQSL. While the number of DEGs in SQSL and XZX45 after 1 day of cold treatment did not show a significant difference, substantial differences were observed after 3 and 5 days of cold treatment ( Figure 3B,C). These results indicate that cold stress had a more pronounced and sustained impact on gene expression levels in XZX45 compared to SQSL, suggesting that SQSL exhibited stronger tolerance and adaptability to cold stress than XZX45. A total of 9899 genes were detected in both genotypes across the three time points, with 8144 genes in SQSL and 6427 genes in XZX45 ( Figure 3B). Among these DEGs, we identified 4672 DEGs common to both genotypes (G2), with 1755 genes specifically differentially expressed in SQSL (G3) and 3472 genes specifically differentially expressed in XZX45 (G1) (Figures 3B and S2). Overall, these findings provide insights into the transcriptional dynamics and differential gene expression patterns between SQSL and XZX45 under cold stress conditions, highlighting the contrasting responses and potential molecular mechanisms underlying their cold stress tolerance. (G2), with 1755 genes specifically differentially expressed in SQSL (G3) and 3472 genes specifically differentially expressed in XZX45 (G1) (Figures 3B and S2). Overall, these findings provide insights into the transcriptional dynamics and differential gene expression pa erns between SQSL and XZX45 under cold stress conditions, highlighting the contrasting responses and potential molecular mechanisms underlying their cold stress tolerance.

Role of Transcription Factors under Cold Stress
To identify cold-responsive transcription factors (TFs) in the regulatory network of cold stress, we utilized the PlantTFDB database [29], which contains information on known or annotated TF genes in the rice genome. Out of the 1,759 TF genes belonging to 56 different TF families, a total of 572 differentially expressed TFs from 48 TF families, including WRKY, NAC, bHLH, ERF, bZIP, MYB, C2H2, and GRAS, were identified ( Figure  4, Table S2). Among these TFs, 192 were present in G1, 294 in G2, and 86 in G3 ( Figure 4, Table S2).

Role of Transcription Factors under Cold Stress
To identify cold-responsive transcription factors (TFs) in the regulatory network of cold stress, we utilized the PlantTFDB database [29], which contains information on known or annotated TF genes in the rice genome. Out of the 1,759 TF genes belonging to 56 different TF families, a total of 572 differentially expressed TFs from 48 TF families, including WRKY, NAC, bHLH, ERF, bZIP, MYB, C2H2, and GRAS, were identified ( Figure 4, Table S2). Among these TFs, 192 were present in G1, 294 in G2, and 86 in G3 ( Figure 4, Table S2).
ing exclusively induced and Os01g0935200 being repressed after 1, 3, and 5 days of cold treatment (Table S2).
These findings provide insights into the differential expression pa erns and specific TF families involved in the cold stress response, highlighting the potential roles of WRKY, NAC, bHLH, ERF, bZIP, MYB, C2H2, and GRAS TFs in the regulatory network underlying cold stress tolerance in rice.

Gene Ontology and Pathway Enrichment Analysis of DEGs
Gene Ontology (GO) enrichment analysis was conducted to evaluate the enrichment of differentially expressed genes (DEGs) in groups G1, G2, and G3, where the statistical significance threshold was set at p.adjust < 0.05 (Tables S3-S5). Our analysis identified a total of 354 GO terms encompassing diverse biological processes (BP), cellular components (CC), and molecular functions (MF). Specifically, within these GO terms, 231 were associated with BP, 40 with CC, and 82 with MF. Notably, we present the top 18 GO terms for each group ( Figure 5A-C).

Gene Ontology and Pathway Enrichment Analysis of DEGs
Gene Ontology (GO) enrichment analysis was conducted to evaluate the enrichment of differentially expressed genes (DEGs) in groups G1, G2, and G3, where the statistical significance threshold was set at p.adjust < 0.05 (Tables S3-S5). Our analysis identified a total of 354 GO terms encompassing diverse biological processes (BP), cellular components (CC), and molecular functions (MF). Specifically, within these GO terms, 231 were associated with BP, 40 with CC, and 82 with MF. Notably, we present the top 18 GO terms for each group ( Figure 5A-C).
To investigate the biological processes occurring during cold stress treatment, we utilized the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database [30]. Specifically, we performed KEGG pathway analysis and compared gene clusters' functional profiles for the differentially expressed genes (DEGs) in groups G1, G2, and G3, where the statistical significance threshold was set at p.adjust < 0.05 (Table S6). The analysis revealed the involvement of 413 DEGs across 35 pathways in the G1, G2, and G3 groups (Table S6).
Interestingly, we observed that "alpha-Linolenic acid metabolism" (dosa00592) was a common pathway between the G1 and G2 groups, "Galactose metabolism" (dosa00052) was a common pathway between the G1 and G3 groups, and "Motor proteins" (dosa04814) was a common pathway between the G2 and G3 groups ( Figure 5D, Table S6).

Analysis of the DEGs in GO Category: "Response to Cold" (GO:0009409) in G3 Group
We have identified a total of 1755 differentially expressed genes (DEGs) specific to the cold-tolerant genotype SQSL, referred to as the G3 group ( Figure 3B). Our Gene Ontology (GO) enrichment analysis revealed a significant enrichment of DEGs in the G3 group associated with the biological process "response to cold" (GO:0009409) ( Figure 4C), which was not observed in the G1 and G2 groups ( Figure 4A,B). Notably, we identified 27 DEGs enriched in the "response to cold" term (GO:0009409), among which 14 DEGs, including Os12g0137500, Os05g0459000, Os06g0164400, Os03g0719100, Os02g0686100, and Os09g0522000, were exclusively induced, while 13 DEGs, including Os01g0857200, Os01g0249200, Os08g0499300, Os03g0452300, Os02g0781400, Os08g0151800, Os03g0785900, and Os08g0441500, were repressed following cold stress treatment ( Figure 6A, Table 1). These findings suggest that these specific DEGs associated with the "response to cold" (GO:0009409) are actively involved in the response to cold stress in the cold-tolerant genotype SQSL. Additionally, we discovered that the OsCOLD1 (Os04g0600800) gene, which encodes a regulator of G-protein signaling localized on the plasma membrane and endoplasmic reticulum (ER), contributes to cold tolerance in rice (Table 1) [19]. Our findings revealed that OsCOLD1 exhibited higher expression levels in SQSL compared to XZX45 after 3 and 5 days of cold treatment ( Figure 6B, Table 1). In response to cold stress, Drought Responsive Element Binding (DREB) genes are rapidly induced. DREB factors, which belong to the APETALA2/ethylene-responsive factor (AP2/ERF) transcription factor subfamily, can recognize dehydration-responsive element/C-repeat (DRE/CRT) cis-elements. This recognition triggers the expression of cold stress-related genes, thereby enhancing cold tolerance [10]. In our study, we observed up-regulation of OsDREB1B (Os09g0522000), encoding the dehydration-responsive element-binding protein 1B, which confers chilling tolerance in rice. Specifically, OsDREB1B was up-regulated during the later response phases at 3 and 5 days of cold treatment in SQSL (Table 1, Figure 6A). Furthermore, the expression level of OsDREB1B in SQSL was higher than that in XZX45 at these time points ( Figure 6B), suggesting its potential role in conferring cold tolerance during the late response phase in rice.
We also identified the gene OsSIZ2 (Os03g0719100), encoding SUMO E3-ligases, involved in stress response and stress adaptation [32]. In our study, OsSIZ2 was up-regulated at 1, 3, and 5 days of cold treatment in SQSL. Notably, the expression level of OsSIZ2 in SQSL was higher than that in XZX45 after 5 days of cold treatment (Table 1, Figure 6B). DEGs between SQSL and XZX45 at different stages of cold stress. The data are presented as means ± SEM (standard error of the mean) from three independent replicates. Statistical significance is denoted as follows: ns (not significant) for p ≥ 0.05, * for p < 0.05, ** for p < 0.01, *** for p < 0.001, and **** for p < 0.0001. Moreover, several genes, including OsbHLH108 (Os06g0164400), OsABA3 (Os06g0670000), OsWRKY30 (Os08g0499300), OsRAB16B (Os11g0454200), and two TAP42like family genes (Os11g0141000, Os12g0137500), exhibited exclusive differential expression pa erns under 1, 3, and 5 days of cold stress treatment. These genes also showed differential expression levels between SQSL and XZX45 ( Figure 6A,B, Table 1). These results indicate that the DEGs enriched in the "response to cold" (GO:0009409) gene DEGs between SQSL and XZX45 at different stages of cold stress. The data are presented as means ± SEM (standard error of the mean) from three independent replicates. Statistical significance is denoted as follows: ns (not significant) for p ≥ 0.05, * for p < 0.05, ** for p < 0.01, *** for p < 0.001, and **** for p < 0.0001. To further explore the expression patterns of these DEGs, we compared their expression levels between the cold-susceptible genotype XZX45 and the cold-tolerant genotype SQSL following 1, 3, and 5 days of cold treatment. Among these DEGs, we identified Os-TRXh1 (Os07g0186000) ( Figure 6A,B, Table 1), a member of the thioredoxin family encoding the h-type Trx protein. Previous studies have demonstrated the involvement of h-type Trx protein-encoding genes in the cold stress response of rice seedlings [24,31]. Our results revealed that the expression of OsTRXh1 was up-regulated consistently at 1, 3, and 5 days of cold stress treatment (Table 1). Interestingly, we observed no significant difference in expression levels between SQSL and XZX45 after 1 day of cold stress. However, at 3 and 5 days of cold treatment, the expression levels of OsTRXh1 were significantly higher in SQSL compared to XZX45 ( Figure 6B).
Additionally, we discovered that the OsCOLD1 (Os04g0600800) gene, which encodes a regulator of G-protein signaling localized on the plasma membrane and endoplasmic reticulum (ER), contributes to cold tolerance in rice (Table 1) [19]. Our findings revealed that OsCOLD1 exhibited higher expression levels in SQSL compared to XZX45 after 3 and 5 days of cold treatment ( Figure 6B, Table 1). In response to cold stress, Drought Responsive Element Binding (DREB) genes are rapidly induced. DREB factors, which belong to the APETALA2/ethylene-responsive factor (AP2/ERF) transcription factor subfamily, can recognize dehydration-responsive element/C-repeat (DRE/CRT) cis-elements. This recognition triggers the expression of cold stress-related genes, thereby enhancing cold tolerance [10]. In our study, we observed up-regulation of OsDREB1B (Os09g0522000), encoding the dehydration-responsive element-binding protein 1B, which confers chilling tolerance in rice. Specifically, OsDREB1B was up-regulated during the later response phases at 3 and 5 days of cold treatment in SQSL (Table 1, Figure 6A). Furthermore, the expression level of OsDREB1B in SQSL was higher than that in XZX45 at these time points (Figure 6B), suggesting its potential role in conferring cold tolerance during the late response phase in rice.
We also identified the gene OsSIZ2 (Os03g0719100), encoding SUMO E3-ligases, involved in stress response and stress adaptation [32]. In our study, OsSIZ2 was up-regulated at 1, 3, and 5 days of cold treatment in SQSL. Notably, the expression level of OsSIZ2 in SQSL was higher than that in XZX45 after 5 days of cold treatment (Table 1, Figure 6B).

Validation by qRT-PCR
We conducted validation of the RNA-seq data's accuracy and reliability by assessing the expression levels of 7 randomly selected genes. The expression profiles of these genes, determined through quantitative real-time polymerase chain reaction (qRT-PCR), exhibited a high degree of similarity with the results obtained from RNA-seq analysis ( Figure S3, Table S7).

Discussion
Given the constraints imposed by low environmental temperatures, plants have evolved diverse response mechanisms to cope with cold stress. Rice, originating from tropical or subtropical regions, is particularly susceptible to cold stress, which severely limits its productivity, an important consideration for ensuring sustainable rice production [33][34][35]. Asian cultivated rice (Oryza sativa) has undergone domestication from its wild relatives, namely Oryza nivara and O. rufipogon. It comprises two major subspecies, indica (O. sativa ssp. indica) and japonica (O. sativa ssp. japonica) [34,35]. The temperate japonica cultivars, a specific type of japonica, are commonly cultivated in regions with lower average temperatures and exhibit higher tolerance to chilling stress compared to indica cultivars [34,35]. While there have been several studies investigating the molecular genetics underlying cold tolerance during the seedling stage between japonica and indica cultivars [23][24][25]27,36,37], limited research has focused solely on indica cultivars. In this study, two indica cultivars, SQSL, which was classified as very highly tolerant, and XZX45, which was highly susceptible to chilling stress (Figure 1), were utilized to investigate the mechanisms involved in the cold stress response of rice using RNA-seq analysis.
A total of 9899 DEGs in both the contrasting genotypes were identified, of which 8144 genes and 6427 genes were observed in XZX45 and SQSL, respectively (Figure 3), which suggests that cold stress had a more pronounced and sustained impact on gene expression levels in XZX45 than in SQSL, indicating that the tolerance and adaptability of SQSL to cold stress were stronger than those of XZX45. Previous studies on cold and salinity tolerance also revealed a relatively different number of DEGs in contrasting genotypes [23][24][25][26]. Of 9899 DEGs, we found 4672 DEGs (G2) were common to SQSL and XZX45, 1755 (G3) and 3472 genes (G1) were only differentially expressed in SQSL and XZX45, respectively, which suggested that these 1755 specific genes may be continuously leading the tolerance to cold stress in SQSL ( Figure 3B).
Transcription factors (TFs) significantly regulate various aspects of plant development, including propagation, maturation, and responses to unfavorable environmental conditions such as drought, cold, salinity, and high temperature [38,39]. Numerous studies have demonstrated the crucial roles played by different types of transcription factors (TFs) in various functional pathways within plants. For instance, WRKY TFs play critical roles in plant responses to biotic and abiotic stresses [40]; NAC TFs are involved in controlling numerous aspects of plant development, such as embryo development, leaf senescence, lateral root formation, and response to abiotic stresses [41], bZIP TFs regulate processes including pathogen defense, light and stress signaling, seed maturation, and flower development [42], bHLH TFs have been found to be involved in phytochrome signaling activity [43], ERF TFs have been implicated in diverse responses to environmental stimuli [44], and MYB TFs have been identified as key players in the control of anthocyanin biosynthesis [45]. We have identified a total of 572 differentially expressed TFs related to 48 TF families, including WRKY, NAC, bHLH, ERF, bZIP, MYB, C2H2, and GRAS ( Figure 4, Table S2). Notably, the bZIP family had the highest number of TF DEGs in G1 (19), while the WRKY and bHLH families had the highest numbers in G2 (34 WRKYs) and G3 (9 bHLHs), respectively ( Figure 4, Table S2). These results suggest that these TFs might be involved in response to cold stress in rice.

Plant Materials and Cold Treatment
The cold-tolerant genotype SQSL and cold-susceptible genotype XZX45 were used in this study. Mature non-dormant seeds of the two indica cultivars were distilled with 2% H 2 O 2 for 30 min and then rinsed three times with distilled water. Rice seeds were transferred to the germination box in a growth chamber for 3 days at 22 • C in the dark. The germinated seeds were transferred to 96-well plates and then grown hydroponically. The 96-well plates were placed in a plant growth chamber (14 h-light/10 h-dark conditions) with temperatures of 28 • C and 25 • C for the light and dark conditions, respectively. The three-leaf seedlings were transferred to a growth chamber and exposed at a temperature of 8 • C (treatment) or 30 • C (control) with a 12 h light:12 h dark photoperiod as the temperature treatment. Rice roots at 1, 3, and 5 days later under control and cold treatment were harvested with three technological and biological replicates and then stored at −80 • C for RNA extraction and sequencing.
In this study, we utilized the cold-tolerant genotype SQSL and the cold-susceptible genotype XZX45. To initiate the experiment, mature non-dormant seeds of the two genotypes were subjected to a 30-min treatment with 2% H 2 O 2 , followed by rinsing with distilled water three times. Subsequently, the rice seeds were transferred to a germination box within a growth chamber and kept in the dark at a temperature of 22 • C for a period of 3 days.
Upon germination, the sprouted seeds were carefully transferred to 96-well plates for hydroponic cultivation. The 96-well plates were placed in a plant growth chamber (14 h-light/10 h-dark conditions) with temperatures of 28 • C and 25 • C for the light and dark conditions, respectively. The environmental conditions were maintained consistently throughout the experiment.
After the seedlings had developed three leaves, they were transferred to a growth chamber for the subsequent temperature treatment. The temperature conditions were set at 8 • C for the cold treatment group and 30 • C for the control group, both under a 12-h light and 12-h dark photoperiod. These conditions represented the respective cold and optimal temperature regimes.
To analyze the effects of the temperature treatment on gene expression, rice roots were harvested at 1, 3, and 5 days after the initiation of the control and cold treatments. Each sampling time point was conducted with three replicates, encompassing both technological and biological replicates. The harvested root samples were immediately stored at −80 • C for subsequent RNA extraction and sequencing procedures.

RNA Extraction and Sequencing
To prepare the root samples for RNA extraction, fresh root samples were ground using mortars pre-cooled with liquid nitrogen. Total RNA extraction was carried out using Trizol reagents (CAT. No. 15596018, Invitrogen, State of California, USA) according to the manufacturer's protocol. The extracted RNA samples were further purified using the Rneasy Mini Kit (CAT. No. 74104, Qiagen, Dusseldorf, Germany) to remove contaminants and impurities. The concentration and quality of the purified RNA samples were assessed using a Colibri instrument (LB 915, Titertek Berthold, Pforzheim, Germany), which measures parameters such as RNA concentration, purity, and integrity. For sequencing, a total of 36 libraries were prepared, with each library representing a separate sample or condition. The libraries were sequenced using the BGISEQ-500 sequencer, which is a high-throughput sequencing platform capable of generating large amounts of sequencing data.

RNA-Seq Analysis
The raw sequencing reads obtained from the BGISEQ-500 sequencer underwent data preprocessing steps. Firstly, adapter trimming and removal of low-quality reads were performed using the fastp (v0.23.2) [47]. The resulting clean reads were then aligned to the genome sequence of Oryza sativa ssp. japonica cv. Nipponbare (IRGSP-1.0 2022-09-01, https://rapdb.dna.affrc.go.jp) using the HISAT2 (v2.2.1) [48]. To determine the read counts mapped to each gene, the featureCounts, tool was utilized [28]. Additionally, the transcripts per million (TPM) values for each gene were calculated using in-house R scripts. Genes with a TPM value of at least 1 were considered as expressed. Differential expression analysis was performed using the DESeq2 R package (v1.40.1) [49]. The analysis aimed to identify genes showing significant differential expression under cold stress conditions. The criteria for differential expression were set as a false discovery rate (FDR) of less than 0.05 and an absolute log2foldchange greater than 1. To gain insights into the functional annotations and pathways associated with the differentially expressed genes (DEGs), Gene Ontology (GO) enrichment analysis and KEGG pathway enrichment analysis were conducted using the ClusterProfile R package (v 4.7.1) [50]. The enrichment analysis was implemented with the false discovery rate (FDR) correction. To identify transcription factors (TFs) among the DEGs, the PlantTFDB database (http://planttfdb.gao-lab.org) [29] was utilized.

qRT-PCR Validation
The RNA-seq data were validated via RT-qPCR analyses using an Applied Biosystems StepOnePlus™ Real-Time System with the SYBR1 Select Master Mix (2×) (ABI, Los Angeles, CA, USA). The reaction protocol was as follows: 95 • C for 10 min; 45 cycles at 95 • C for 10 s, 60 • C for 30 s, and 72 • C for 20 s; and then 72 • C for 5 min. All the specific primers used are presented in Table S7. Three technical replicates were included per sample. The rice Ubiquitin (Os03g0234200) gene was used as an internal standard. The relative expression values of the different genes were calculated using the 2 −∆∆Ct method.

Conclusions
In this study, we conducted an investigation into the transcriptome dynamics of the roots of two indica rice cultivars, SQSL and XZX45, under optimal temperature and cold stress conditions using RNA-seq technology. Our analysis focused on identifying differentially expressed genes (DEGs), the transcription factors associated with cold response, conducting GO analysis, KEGG pathway analysis, and analyzing the expression profiles of cold-response-related genes.
Through our comprehensive analysis, we discovered crucial DEGs and transcription factors that are specifically related to cold response. These findings provide insights into the molecular mechanisms underlying cold tolerance during the seedling stage in SQSL and XZX45. Notably, a significant proportion of the DEGs specific to SQSL were enriched in the "response to cold" gene ontology category (GO:0009409). This enrichment suggests that these DEGs likely play a critical role in the cold response mechanisms in rice.
Overall, our study contributes to a systematic understanding of the genetic and molecular basis of cold responses in rice. The knowledge gained from this study will serve as a foundation for future endeavors aimed at improving the cold tolerance of cold-sensitive rice varieties.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12142675/s1. Figure S1: The TPM distribution across all samples. Figure S2: Heatmap depicting the expression levels of the DEGs in G1, G2, and G3. The heatmap is constructed based on TPM (Transcripts Per Million) values across all samples. The color spectrum ranging from blue to red represents the gene expression levels from low to high. Figure S3: qRT-PCR verification of expression profiles obtained through RNA-seq. Table S1: Summary of RNA-Seq reads and their mapping on the rice genome. Table S2: Summary of transcription factors (TFs) analysis for DEGs in G1, G2, and G3. Table S3: GO enrichment for DEGs in G1. Table S4: GO enrichment for DEGs in G2. Table S5: GO enrichment for DEGs in G3. Table S6: KEGG pathway enrichment comparison among G1, G2, and G3. Table S7: Primers used for qRT-PCR validation of selected genes.
Author Contributions: T.Y.: investigation, writing-original draft. W.T. and G.Z.: conceptualization, supervision. T.Y., M.S., R.S. and X.W.: bioinformatics analysis, including data processing, statistical analysis, and interpretation of the results. X.L. (Xuedan Lu), Y.X., H.D., X.L. (Xiong Liu) and G.Z.: writing-review and editing. All authors made significant contributions to the article by providing intellectual input and reviewing and editing the manuscript. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: All raw sequencing reads generated in this study have been deposited in the public database of the National Genomics Data Center under the data accession number PRJCA017633. All the other data are available from the corresponding authors upon request.

Conflicts of Interest:
The authors declare that they have no competing interest or conflict of interest that could be perceived as having influenced the research conducted in this study.