Comparative transcriptomic analysis of germinating rice seedlings to individual and combined anaerobic and cold stress
BMC Genomics volume 24, Article number: 185 (2023)
Rice is one of the most important cereals consumed worldwide. Two major abiotic factors affecting rice plants in different growth stages are flooding stress and cold stress. These abiotic stresses can take place independently or simultaneously and significantly affect rice plants during germination and seedling growth. Fortunately, a wide array of phenotypic responses conferring flooding stress and chilling stress tolerance exist within the rice germplasm, indicating the presence of different molecular mechanisms underlying tolerance to these stresses. Understanding these differences may assist in developing improved rice cultivars having higher tolerance to both stresses. In this study, we conducted a comparative global gene expression analysis of two rice genotypes with contrasting phenotypes under cold stress, anaerobic stress, and combined cold and anaerobic stress during germination.
The differential gene expression analysis revealed that 5571 differentially expressed genes (DEGs), 7206 DEGs, and 13279 DEGs were identified under anaerobic stress, cold stress, and combined stress, respectively. Genes involved in the carbohydrate metabolic process, glucosyltransferase activity, regulation of nitrogen compound metabolic process, protein metabolic process, lipid metabolic process, cellular nitrogen compound biosynthetic process, lipid biosynthetic process, and a microtubule-based process were enriched across all stresses. Notably, the common Gene Ontology (GO) analysis identified three hub genes, namely Os08g0176800 (similar to mRNA-associated protein mrnp 41), Os11g0454200 (dehydrin), and OS10g0505900 (expressed protein).
A large number of differentially expressed genes were identified under anaerobic, cold conditions during germination and the combination of the two stress conditions in rice. These results will assist in the identification of promising candidate genes for possible manipulation toward rice crops that are more tolerant under flooding and cold during germination, both independently and concurrently.
Rice is a staple food crop for more than half of the world’s population. Many rice-growing regions, especially in South and Southeast Asia, are in flood-prone lowland areas . Flash flooding due to unpredicted rainfall is a common phenomenon in such areas, severely affecting rice crops in various stages of their development, requiring different strategies to address the problems [2,3,4,5,6,7]. Rice is adapted to aquatic ecology; therefore, it has a unique mechanism of germinating and developing a coleoptile underwater . The coleoptile elongation is achieved by rapid elongation of basal cells (up to 200 um in 12 h) immediately after emerging from an embryo . However, not all rice accessions possess vigorous growth underwater, and the germination potential varies greatly among rice genotypes under flooded conditions.
Another abiotic stress that is detrimental to rice germination and growth is cold stress. Cold stress alone or its combination with unexpected heavy rains at the beginning of rice-growing seasons are common in rice-growing areas, especially in temperate regions, including the U.S. The recommended date of rice planting in Texas is March 15 to April 15 (https://beaumont.tamu.edu/elibrary/Newsletter/2010_April_Newsletter.pdf), but the yield can be significantly improved if early planting is successful. The main challenges for early planting in Texas are cold weather and unexpected rain. These stresses cause poor emergence and poor establishment of the germinating rice seeds. Robust seedling establishment is necessary to improve yield stability under stress-prone conditions, and tolerance of anaerobic germination (AG) and cold temperature are the most desirable traits. Understanding the molecular mechanisms underlying tolerance to these stresses and subsequently using this knowledge to assist in developing rice varieties tolerant to these conditions can potentially improve rice seedling establishment, which translates to more robust and productive rice growth.
Many studies on QTL mapping and genome wide association studies (GWAS) on AG tolerance in rice have been reported [10,11,12,13,14,15,16,17,18,19,20,21,22]. Notably, a gene underlying one of the QTLs, qAG-9-2, has been cloned as trehalose-6-phosphate-phosphatase (OsTPP7), which is involved in starch mobilization during germination . Previous biochemical and enzymatic experiments have demonstrated the involvement of starch breakdown mechanisms and the induction of amylase during low oxygen germination in rice [24,25,26]. A positive correlation between coleoptile length and total amylolytic activities has been reported . Under anoxic conditions, the level of ethanol is elevated in the fast-growing coleoptiles, which suggests that energy is generated from the fermentative metabolism pathway [28,29,30]. Recent studies have also demonstrated that maintaining a high rate of energy production as well as fluxes between the glycolytic and fermentative pathways is crucial for survival during anaerobic conditions [31, 32].
A few studies on genome-wide transcriptome analysis to investigate the gene expression analysis of rice coleoptiles under hypoxic or anoxic conditions have been reported. The results of the studies showed that some common molecular mechanisms are involved in coleoptile growth, including carbohydrate metabolism, fermentation, hormone induction, and cell division and expression [33,34,35]. A study on the tip and basal segments of oxygen-deprived coleoptiles of rice showed that the gene expression is region-specific. Furthermore, the metabolic activities are different along the length of the coleoptile under normoxic (air), hypoxic (3%), and anoxic conditions . Tolerant plants have developed the ability to generate ATP in the absence of oxygen through fermentative metabolism and to develop specific morphologies (such as air channels and enhanced shoot elongation) that improve the entrance of oxygen [29, 36,37,38,39,40,41,42,43]. Moreover, gene expression studies on plants exposed to low oxygen conditions revealed the up-regulation of genes coding for signal transduction components , transcription factors [45, 46], non-symbiotic hemoglobin , ethylene biosynthesis , nitrogen metabolism , and cell wall loosening .
Similarly, several studies on the identification of cold tolerance QTLs and genes during germination and seedling stages using biparental QTL mapping and GWAS have been reported [51,52,53]. A QTL for cold tolerance, COLD1, was identified on the long arm of chromosome 4. The causal gene encodes a regulator of G-protein signaling, which together with the G-protein α subunit plays a key role in activating the Ca2+ channel for sensing low temperature and accelerating G-protein GTPase activity, contributing to the adaptation of japonica rice to chilling conditions . Similarly, it was reported that the transcription factor bZIP73 may have contributed to japonica adaptation to cold climates . A gene underlying a QTL for cold tolerance at the booting stage, CTBa, encoding a conserved leucine-rich repeat, has also been cloned .
Several transcriptomic studies to investigate the underlying mechanisms of cold tolerance in rice using microarray technologies have been reported [57,58,59,60,61]. These studies revealed that a large number of cold-responsive genes could be categorized into three groups. The first group is the signaling components that regulate gene expressions in cold stress responses, which include protein kinases and transcription factors. Secondly, the functional components that directly protect plant cells against cold stress. These include enzymes in metabolic pathways, aquaporins and proteins involved in photosynthesis. Lastly, the small noncoding RNAs such as the microRNAs (miRNAs).
A comprehensive understanding of the molecular mechanisms of anaerobic stress and cold adaptation is an important area of research. The tolerance or sensitivity to these abiotic factors is a complex phenomenon involving responses at the molecular, biochemical, and physiological levels. However, limited research has been performed in comparative genome-wide transcriptomics analysis across these stresses, especially when they occur together. The goal of this study was to evaluate the differentially expressed genes in two contrasting rice cultivars in response to anaerobic, cold, and a combination of anaerobic and cold stress during germination to better understand the underlying molecular mechanisms and metabolic pathways that contribute to tolerance.
Results and discussion
Overview of cDNA library sequencing
Tissues from four replicates per sample were harvested to construct cDNA libraries of germinating rice seeds from two contrasting genotypes, Darij (tolerant to anaerobic germination (AG) and cold stress conditions) and line 4610 (susceptible), grown under control, AG stress, cold stress, and combined cold and AG stress conditions. The cDNA libraries were then sequenced with the HiSeq 4000 using 2X75 . An average of 17.62, 23.00, and 20.53 million clean reads were generated from AG stress, cold stress, and combined AG and cold stress and mapped to the Nipponbare reference genome with Hisat2 , with 91.62% to 96.78% clean reads from all samples successfully mapped to the reference genome (Supplementary Tables 1, 2, 3, 4 and 5). The total gene counts were normalized using the DESeq2 package, and the expression signals of each gene were calculated. The PCA analysis based on regularized log transformation of total gene counts with DESeq2 showed 58%, 77%, and 59% variance explained by PC1 under AG, cold and combined AG and cold stress, respectively, and 34%, 16%, and 35% variance explained by PC2 under AG, cold and combined AG and cold stress, respectively (Supplementary Figure 1).
Based on the phenotypic evaluation, the two selected lines showed contrasting differences under AG stress, cold stress and combined AG and cold stress (Table 1). Under control conditions, the mean coleoptile length of Darij was 1.54 cm and 4610 was 1.2 cm at 4 days after sowing (DAS) and 2.48 cm and 3 cm, respectively, at 7 DAS. There was no statistically significant difference (P value < 0.05) in the coleoptile length between these two lines under control conditions. In contrast, a statistically significant difference (P value < 0.05) was observed among these lines under stress treatments. Under AG stress condition, the mean coleoptile length of Darij was 2.2 cm and 4610 was 0.25 cm at 4 DAS. Under cold conditions, the coleoptile length of Darij was 0.41 and 4610 was 0.025 cm at 21 DAS and 0.86 cm and 0.14 cm, respectively, at 28 DAS. The mean plumule length of Darij and 4610 after 4 days of recovery at 30°C was 2.63 cm and 0.63 cm, respectively. Under combined stress of AG and cold, the mean coleoptile length of Darij was 0.43 cm and 0.68 cm at 28 DAS and 35 DAS, respectively, whereas the susceptible line, 4610, did not germinate at all, even at 35 DAS. All variables (coleoptile length under AG stress, shoot length under cold exposure, plumule growth recovery, and coleoptile length under the combination of AG and cold stress) showed superior performance of the tolerant line (Darij) over the susceptible line (4610). The tolerant line was superior for shoot length under AG stress, cold stress, and the combined stress, and the seedling growth was increased in the tolerant line even after 21 days of cold exposure. Likewise, the plumule growth recovery of the tolerant line was also much better than the susceptible one. On the other hand, the seedling growth of the susceptible line was negligible under cold stress and showed no germination under combined AG and cold stress. Similarly, the tolerant line showed significantly longer coleoptile growth under the submerged condition in 4 DAS. Alpi and Beevers, 1983  and Kawai and Uchimiya, 2000  reported that when flooding occurs after direct seeding, tolerant lines germinate faster. At the same time, their coleoptiles grow faster, enabling them to reach the more aerated zone and oxygen referred to as the snorkel effect. The significant reduction of coleoptile growth under cold and submerged conditions during the initial stage of germination in rice is a limiting factor for good crop establishment [66, 67]. Thus, selecting genotypes with longer coleoptile lengths is critical for cold and flooding tolerance during germination [68, 69].
Differential expression analysis
To further examine whether the gene expression was affected by treatment only or both genotype and treatment, multi-factorial linear model testing was performed on our whole transcriptome dataset using DESeq2. Applying an FDR corrected p-value of < 0.05, the separate analysis of Darij and 4610 under AG stress versus control showed that the total number of differentially expressed genes (DEGs) was 16,383 in Darij and 15,435 genes in line 4610, with 9,798 DEGs in common (Fig. 1a). In this analysis, the total number of upregulated genes was 8,082 in Darij and 7,579 in 4610, while the total number of downregulated genes was 8,744 in Darij and 7,856 in 4610 (Supplementary Table 6). Under cold stress, 12,101 and 16,240 DEGs were identified in Darij and 4610, respectively, compared to the control, with 7,962 DEGs in common. The total number of upregulated genes in Darij and 4610 were 5,612 and 8,010, respectively, and the total number of downregulated genes in Darij and 4610 were 6,489 and 8,230, respectively. For the combination of AG and cold stress, 19,040 and 21,906 DEGs were identified in Darij and 4610, respectively, with 14,382 DEGs in common. The total number of upregulated genes in Darij and 4610 were 8,435 and 10,590, respectively, and the total number of downregulated genes in Darij and 4610 were 10,605 and 11,316, respectively. The total number of DEGs were higher in both lines under AG stress and under the combination of AG and cold stress. The number of DEGs was higher in the susceptible line, 4610, under cold stress and the combination of AG and cold stress. Similar findings were previously reported [70, 71]. Da Maia et al. 2016 revealed 19 times more DEGs in cold-sensitive genotype than in cold-tolerant genotype in response to cold, and Pan et al. 2020 identified more downregulated DEGs in cold-sensitive rice. This could be due to severe damage imposed on the susceptible lines at the cellular and molecular level than on the tolerant lines.
We found statistically significant interaction effects indicating that the gene expression under different stresses is regulated differently depending on the genetic background, i.e., variation of coleoptile elongation in diverse rice accessions could be potentially determined by genotype-specific gene expression. Considering both genotype and treatment effects, 5,571 DEGs were identified under AG stress where 2,294 genes were up-regulated, and 3,277 genes were down-regulated. Under cold stress, 7,206 DEGs were identified, where 3,276 genes were up-regulated, and 3,930 genes were down-regulated. Under combined stresses, 13,279 DEGs were identified with 6,282 up- and 6,997 down-regulated genes (Supplementary Table 7).
From the pairwise comparison, the highest number of common DEGs (5,010) were found between cold and combined AG and cold stress, followed by DEGs (2,399) between AG and combined AG and cold stress and DEGs (1,403) between AG stress and cold stress (Fig. 1). The results of shared differentially expressed genes indicated a greater relationship between cold stress and combined AG and cold stress, followed by AG stress and combined AG and cold stress. Among all DEGs, 992 genes were shared among all three stresses, of which 430 genes were commonly up-regulated, 195 genes were commonly down-regulated, and 367 of them were up-regulated under one condition and down-regulated under the other two conditions or inversely. Among 5,010 common DEGs between cold and combined AG and cold stress, 2,112 genes were up-regulated, 2,726 genes were down-regulated, and 172 genes were up-regulated under one condition and down-regulated under the other condition. Among 2,399 common DEGs between the AG and combined AG and cold stress, 515 genes were up-regulated, 901 genes were down-regulated, and 983 genes were up-regulated under one condition and down-regulated under the other condition. Among 1,403 common DEGs between AG and cold stress, 289 genes were up-regulated, 629 genes were down-regulated, and 485 genes were up-regulated under one condition and down-regulated under the other condition. The total number of differentially expressed genes was higher in combined cold stress and AG stress than under individual stress indicating that the combined stress could have caused severe damage to the cellular and metabolic processes compared to individual stress, leading to greater complexity of gene regulation.
Common and unique DEGs in the tolerant and susceptible line
In broad GO annotation, we observed some similarities and differences among the tolerant and susceptible lines. The GO enrichment analysis identified a total of 309, 278 and 289 enriched GO terms in the tolerant line, Darij, and a total of 318, 291 and 297 enriched GO terms in the susceptible line, 4610, under AG stress, cold stress and the combination of AG and cold stress, respectively (Supplementary Table 8). There were 247, 220 and 242 common GO terms between the tolerant line and susceptible line under AG stress, cold stress and combination of AG and cold stress, respectively. In the tolerant line, a total of 62, 58 and 47 unique GO were detected under AG stress, cold stress, and a combination of AG and cold stress, respectively. In the susceptible line, a total of 71, 71 and 55 unique GO terms were identified under AG stress, cold stress, and a combination of AG and cold stress, respectively. Under AG stress, the common GO terms included the carbohydrate metabolic process, the nitrogen compound metabolic process, the carboxylic acid metabolic process, the glucose metabolic process, glycolysis, thylakoids, and photosynthesis. On the other hand, the unique GO terms in the tolerant line included the cellular aromatic compound metabolic process, the aromatic compound biosynthetic process, ATPase activity, transferase activity and response to abiotic stimulus; while the unique GO terms in the susceptible line included the microtubule-based process, GTPase activity, the regulation of signaling process, cytoskeletal and peroxidase activity. The molecular functions related to 6-phosphofructokinase activity, ATPase activity, acetyltransferase activity, and transferase activity were uniquely enriched in the tolerant line, while response to oxidative stress, the regulation of signaling process, and GTPase regulator activity were uniquely enriched in susceptible line.
During flooded conditions, rice plants are threatened by a reduction in oxygen availability either due to slower diffusion of oxygen through water or continuing competition for oxygen with respiring microorganisms [36, 72]. Major complications associated with low oxygen stress are related to energy crises caused by the inhibition of mitochondrial oxidative phosphorylation and subsequent reduction in ATP production. The shortage of ATP causes impaired functioning of the plasma membrane H+ pumping ATPase, which leads to hypoxia-induced acidification of cytoplasm . The adaptability of plants to low oxygen stress is generally considered in terms of anaerobic carbon metabolism . The enrichment of the cellular nitrogen metabolic process and regulation of the nitrogen metabolic process in this study show the role of anaerobic carbon metabolism in AG stress conditions. Several researchers have also reported the contribution of nitrogen metabolism in cellular acclimation to low oxygen stress in plants [75, 76]. Low oxygen stress significantly affects several nitrogen metabolism processes, including nitrogen absorption, nitrate and nitrite reduction, nitric oxide (NO) production, ammonium assimilation, and amino acid metabolism . Several studies have shown that nitrate fertilization improves hypoxia tolerance and limits hypoxia stress in the roots [77,78,79]. The higher number of genes involved in the cellular nitrogen compound metabolic process, NAD or NADH binding, oxidoreductase activity, and glucosyltransferase activity might contribute to tolerance. O-glucosyltransferase 2 was previously shown to be involved in Arabidopsis seed germination . Fructokinase (OsFK2) was previously reported to be induced by anoxia . Fructokinase plays a key role in driving the hexose (fructose and UDP-glucose) released by sucrose synthase into glycolysis which helps in increasing tolerance under anaerobic conditions; this might be the reason for increased 6-phosphofructokinase activity in the tolerant line, Darij.
The biological processes related to cellular lipid metabolic process, lipid biosynthetic process, and molecular function related to lipid binding were also uniquely enriched in the tolerant line. The number of genes enriched in lipid localization and lipid transport were higher in the tolerant line than the susceptible line. The hydrogenation of unsaturated fatty acids and the biosynthesis of lipids were hypothesized to be adaptive mechanisms under anaerobic stress because they function as terminal acceptors of respiratory electrons and protons [82, 83] and lipid biosynthesis and accumulation were found to bypass the major problem of anaerobiosis in Echinochloa phyllopogon. However, the rice coleoptile requires anaerobic synthesis of lipids. This synthesis is an adaptive mechanism to achieve the turnover of saturated fatty acids, phospho-, glyco- and neutral lipids under extreme conditions of oxygen deficiency, thereby stabilizing cell membranes under strict anoxia .
The enzyme scavengers, and peroxidase activity enriched exclusively in susceptible line 4610, suggests that reactive oxygen species (ROS) metabolism was higher in this line. Additionally, the metabolism carried out by mitochondria is directly connected to oxygen availability. Anoxia-germinated seedlings of rice develop coleoptiles with intact structure and functional activity of mitochondria capable of oxidative phosphorylation . The GO process involving the cellular component for the mitochondrial membrane is also exclusively enriched in the susceptible line. This may suggest that there were higher activities of the genes related to mitochondria in susceptible genotypes. Oxygen deprivation might have affected mitochondrial ultrastructure leading to susceptibility of genotypes under submerged conditions.
The genes involved in photosynthesis were found to be downregulated under AG stress in both lines; however, the number of downregulated genes in the susceptible line was higher. This agrees with the hypothesis that photosynthetic ability is affected under stresses, including flooding stress . Interestingly, several researchers have reported the upregulation of photosynthesis-related genes under flooding. Nanjo et al. (2011) reported the upregulation of photosynthesis-related genes after imposing flooding at 2 d old seedling stage of soybean . The upregulation of photosynthetic-related genes in Arabidopsis after imposing low oxygen stress for 3h and in poplar after imposing flooding stress for 5 h was also reported . However, Lee et al. (2014) reported the downregulation of photosynthesis-related genes in the leaves of rape seedlings after 3 d of water-logged conditions . This suggests that flood stress stimulates the expression of photosynthesis-related genes during initial stress response and represses the expression of these genes at later stages. Reduced photosynthesis during flooding stress could be due to the production of ROS . High levels of ROS can damage cells through peroxidation of lipids, oxidation of proteins, and other pathways, finally leading to cell death [90, 91]. Moreover, the number of enriched genes in the biological function of heme binding is higher in the tolerant line than in the susceptible line. Non-symbiotic hemoglobins induced under hypoxia were reported to play a role in modulating low oxygen responses as low oxygen concentration induces class-1 hemoglobin gene GLB1 in Arabidopsis [89, 92] which was identified to increase survival under severe hypoxia . Also, non-symbiotic hemoglobins may act as NO-detoxifying agents under hypoxia, where conspicuous amounts of NO are generated .
Under cold stress, some of the common GO categories were carbohydrate metabolic process, nitrogen compound metabolic process, lipid metabolic process, carboxylic acid biosynthetic process, cytoskeletal, and microtubule-associated complex. GO terms such as glucan metabolic process, lipid localization, phospholipid metabolic process, mitochondrial part, ribosome were found to be unique in the tolerant line, while GO terms such as glycolysis, chromosome organization, ATPase activity, heat shock protein binding, and cellular response to stress were found to be unique in the susceptible line. All the common GO terms have a larger number of enriched genes in the susceptible line. High expression of cellulose synthase genes was reported to be important for tolerance to biotic  and abiotic stresses in rice plants . Microtubule organization for cellulose synthase positioning and cellulose deposition has been demonstrated to be essential for biomass production during salt stress in Arabidopsis . In this study, we identified microtubule-based process, microtubule cytoskeleton, microtubule-associated complex, and microtubule motor activity in both tolerant and susceptible lines. However, the number of enriched genes in the susceptible line was higher than in the tolerant line. Additionally, we identified enrichment of genes related to lipid localization, lipid transport, and phospholipid metabolic process exclusively in the tolerant line, which might partly be the reason for higher tolerance. Cell membranes are the major target for temperature susceptibility in plants , and a strong association between temperature and fatty acid content of plant membranes has been demonstrated . An increasing level of saturated fatty acids has been demonstrated to be correlated with higher sensitivity to cold temperatures in Arabidopsis plants , and unsaturated fatty acids, due to their lower melting point, were reported to maintain membrane fluidity .
We also identified a uniquely enriched GO related to cellular component, clathrin coat, in the tolerant line, which was previously reported to play an important role in Ca2+ signaling pathways across the plasma membranes . Q7XKE9 (Clathrin light chain 1) protein was found to be a differentially expressed protein in cold tolerant lines . Chilling stress manifests in a decrease of photosynthesis, an increase of reactive oxygen species (ROS), over-reduction of chloroplast electron transport chain, and damage to membrane integrity which causes cellular dehydration, osmotic imbalance, and other deleterious effects [103, 104]. Our study shows that the tolerant line has more genes expressed for photosynthesis and maintenance of structural membrane integrity. Previous studies have shown that circadian rhythms influence photosynthetic processes and play a vital role in coordinating photosynthetic activity with the diurnal changes in light. The studies have also demonstrated that chilling temperatures disrupt photosynthetic processes in warm-climate plants [105,106,107,108]. Our results confirm previous findings that tolerant lines have developed responses to coping with the changes in membrane compositions to sustain photosynthesis under cold stress conditions.
Under the combination of AG and cold stress conditions, the common GO terms were nitrogen compound metabolic process, carbohydrate metabolic process, signaling process, microtubule-based process, and cytoskeletal process. The unique GO terms in the tolerant line included isoprenoid metabolic process, glucosyltransferase activity, and glutamine family amino acid metabolic process, while the unique GO terms in the susceptible line were GTP binding, GTPase activity and response to stress. Under the combination of AG and cold stress, the common GO terms with higher number of enriched genes in the tolerant line were oxidoreductase activity, peroxidase activity, microtubule motor activity, cellular ketone metabolic process, organic acid metabolic process, oxoacid metabolic process, carboxylic acid biosynthetic process, and lipid biosynthetic process. Other GO terms had more enriched genes in the susceptible line, including cytoskeleton, mitochondrial envelope, microtubule cytoskeleton, hydrolase activity, oxidoreductase activity, protein serine/threonine kinase activity, antioxidant activity, nitrogen compound metabolic process, protein metabolic process, regulation of the cellular metabolic process, regulation of the biosynthetic process, regulation of nitrogen compound metabolic process, regulation of nucleobase, nucleoside, nucleotide, and the nucleic acid metabolic process. The common GOs identified between AG, cold and combined stress included carbohydrate metabolic process, glucosyltransferase activity, regulation of nitrogen compound metabolic process, protein metabolic process, lipid metabolic process, cellular nitrogen compound biosynthetic process, lipid biosynthetic process, microtubule-based process, whereas isoprenoid biosynthetic process, isoprenoid metabolic process, cell redox homeostasis, and the glutamine family amino acid metabolic process were uniquely enriched only in combined stress of AG and cold conditions. Isoprenoids, naturally synthesized volatile compounds, were found to play significant roles in defense mechanisms against different abiotic stresses, including drought, heat, and cold [109, 110]. Isoprenoids alter the ROS response by direct reactions of isoprenoids with ROS, indirect modification of ROS, and stabilization of lipid membranes developing the capability of plants to cope with oxidative stress inside the cells. The induction of 35 genes related to metabolism, transport, signal transduction, and stress responses, including transcription factor genes such as DREB1A, IRO2, and NAC5, with external application of glutamine, was reported . Glutamate was reported to amplify its signal and interact with other signaling pathways to regulate metabolism, growth, and defense responses in rice.
The comparison of enriched GO terms in the tolerant line across different stresses identified the enrichment of ATPase activity, NADP or NADPH binding, photosynthesis and thylakoid part exclusively under AG stress, enrichment of phospholipid metabolic process, starch synthase activity, mitochondrial lumen, regulation of signaling process and GTPase regulator activity exclusively under cold stress, and enrichment of isoprenoid metabolic process, protein-DNA complex, DNA repair, and cell redox homeostasis exclusively under the combination of AG and cold stress. Similarly, the comparison of enriched GO terms in the susceptible line depicted enrichment of photosynthesis, carboxylic acid transport, glucan metabolic process, glucosyltransferase activity, sulfur compound biosynthetic process exclusively under AG stress, enrichment of tetrapyrrole metabolic process, porphyrin metabolic process, hydrolase activity, and ATPase activity exclusively under cold stress, and enrichment of lipid binding, transferase activity, and cellular protein complex assembly under the combination of AG and cold stress.
Genes and transcription factors
Under AG stress, the unique downregulated genes and transcription factors with the highest expression values were OS03G0309400 (pectin esterase) with -11.70 log twofold change (log2FC) value, followed by OS03G0427300 (glutelin), OS05G0305300 (pentatricopeptide), OS04G0121300 (OsSub35 - Putative Subtilisin homologue), OS10G0443800 (MYB family transcription factor), OS07G0147550 (photosystem II 10 kDa polypeptide, chloroplast precursor) ranging from -10.1 to -11.1 log2FC. Whereas the upregulated genes and transcription factors with the highest expression value (ranging from 11.6 to 30 log2FC) were OS10G0475000 (alcohol oxidase) with 30 log2FC value, followed by OS08G0463500 (C2H2 zinc finger protein), OS04G0606450 (ATP-dependent family protein), OS04G0103900 (chalcone synthase), OS06G0297400 (cadmium tolerance factor), OS10G0558600 (pentatricopeptide), OS08G0152400 (cytochrome P450), and OS10G0389300 (red chlorophyll catabolite reductase) (Supplementary Table 7).
Under cold stress, the unique downregulated genes and transcription factors with the highest expression (-7.1 to 12.3 log2FC) were OS10G0178000 (LTPL132, protease inhibitor/seed storage/LTP family protein precursor) with -12.3 log2FC followed by OS09G0457400 (alpha-amylase precursor), and OS04G0295400 (jasmonate-induced protein). The upregulated genes and transcription factors with the highest expression (6 to 11.1 log2FC) were OS07G0639400 (peroxidase precursor), OS07G0520300 (cytochrome P450 domain-containing protein), OS04G0604800 (glycosyl hydrolases family), and OS03G0155500 (expansin precursor).
In combination with AG and cold stress, the unique downregulated genes and transcription factors with the highest expression (-9.5 to -15.3 log2FC) were OS11G0635300 (cytochrome P450), OS09G0551400 (serine/threonine-protein kinase receptor precursor), and OS07G0525900 (chalcone and stilbene synthases). On the other hand, the unique upregulated genes 7.5 to 12.03 log2FC) were OS01G0512200 (AP-3 complex subunit delta), OS07G0431160 (F-box domain-containing protein), OS11G0484700 (Protease inhibitor/seed storage/LTP family protein precursor), OS11G0575600 (lipoxygenase), and OS06G0688200 (beta-expansin precursor).
Under AG stress, transcription factors (TFs) related to aldehyde dehydrogenase, cytochrome P450, dehydration responsive element binding protein (DREB), ethylene response factor (ERF), heat stress transcription factor, heat shock protein, alpha expansin, beta expansin, trehalose-6-phosphate synthase, trehalose-6-phosphate phosphatase, ATP binding cassette (ABC) transporters, AP2 domain-containing proteins, bHLH, bZIP, MYB, C2H2 zinc finger proteins, and WRKY were differentially expressed in both tolerant and susceptible lines (Supplementary Table 6). Under cold stress, the same transcription factors were expressed except for trehalose-6-phosphate phosphatase and beta expansin. The number of TFs from the DREB family and heat stress transcription family was higher. In contrast, TFs from cytochrome P450, aldehyde dehydrogenase, and trehalose-6-phosphate synthase were lower under cold stress than AG stress. Similarly, the same TFs were expressed under combined AG and cold stress while TFs from cytochrome P450, ethylene response factor (ERF), aldehyde dehydrogenase, alpha expansin, and beta expansin were expressed higher than under AG and cold stress. In addition, triose phosphate translocator 2 (OsTPT2) was expressed under AG stress and combined AG and cold stress. The TF, BBR-BPC, was uniquely expressed in both tolerant and susceptible lines under cold stress. The TF, YABBY, was expressed only in the tolerant line under AG stress and cold stress but not under combined AG and cold stress. We found aldehyde dehydrogenase, cytochrome P450, dehydration-responsive element binding protein (DREB), ethylene response factor (ERF), heat stress transcription factor, heat shock protein, alpha expansin, beta expansin, trehalose-6-phosphate synthase, trehalose-6-phosphate phosphatase, ABC transporters, AP2 domain-containing proteins, bHLH, bZIP, MYB, C2H2 zinc finger proteins, and WRKY expressed in both tolerant and susceptible line under AG stress, cold stress and combined AG and cold stress.
The ATP binding cassette (ABC) protein family consists of membrane proteins that are active in the ATP-powered transport of structurally unrelated substrates across the membrane [112,113,114]. The role of ABC transporters in plants in response to abiotic stress has not been widely explored. Previously, Ospdr9, a gene encoding an ABC protein, was reported to be induced by hypoxia stress in rice . In recent years, ABC proteins have been identified as important elements in the transport of phytohormones, heavy metals, lipids, chlorophyll catabolites, and secondary metabolites . ABC proteins in plastids play a key role in the biosynthesis of membrane lipids like galactolipids in the thylakoid membrane . It has been previously reported that AP2/ERF genes, a large multigene superfamily of transcription factors, play an important role in waterlogging and flooding stress [117,118,119]. Hatorri et al. (2009) reported the deep-water response of two genes encoding ethylene response factors, SNORKEL1 and SNORKEL2, in rice . Under deep water conditions, ethylene accumulates in water and induces the expression of both genes, which ultimately lead to internode elongation via gibberellin. Hinz et al. (2010) reported that overexpression of RAP2.2 gene belonging to AP2/ERF transcription factor improves plant survivability under hypoxia stress conditions . MYB-related genes were found to be significantly upregulated in soybean seeding due to flooding stress .
In this study, we identified 13 genes of the MYB family upregulated in the tolerant line and downregulated in the susceptible line under AG stress (Supplementary Figure 2g). The MYB transcription factors contain the MYB domain that is highly conserved in all eukaryotes and involved in regulating a wide range of functions of MYB proteins . Several studies have shown the involvement of MYB transcription factors in regulatory networks associated with plant growth and development [121,122,123]. Many other studies have suggested the significant role of MYB proteins in regulating plant responses to abiotic stress [124,125,126,127]. A study conducted in O. coarcata transcriptomic changes under salt and submergence stresses (alone or combined, compared to control conditions) found several transcription factors, including NAC, WRKY, and MYB gene family members upregulated in leaves, indicating extensive transcriptional regulation in stress responses . MYB-related genes were found to be significantly upregulated in soybean seeding due to flooding stress . A few genes belonging to the MYB family have been previously reported to be induced under cold stress [129, 130]. OsMYBS3 was identified to be repressing the DREB1/CBF-dependent cold signaling pathway in rice . These genes might be involved in a new additional pathway that regulates cold adaptation in rice . Our study identified several differentially regulated genes of the MYB family induced under cold stress in both genotypes. This suggests that MYB family members may play an important role in the regulatory pathway of cold tolerance in rice. The P2/EREB family contains a large number of plant-specific TFs and participates in abiotic stress responses in plants . We also identified many genes in the AP2/EREB family differentially regulated under cold stress.
Many WRKY genes have been found to be transcriptionally regulated in response to abiotic stress treatments in rice [133, 134]. Three DEGs of WRKY genes (Os03g0758900, Os11g0490900 and Os01g0826400) regulated by cold or H2O2 stress in a rice variety have been reported . In our study, we found several genes of the WRKY family upregulated and downregulated in both genotypes under cold stress. Additionally, heat shock proteins have been reported to be involved in heat shock responsiveness in plants . Our study identified a few genes related to heat shock factor and heat shock protein differentially regulated. But the role of heat shock-related proteins in the cold stress response of plants remains obscure.
In this study, we observed 11 genes of Cyst2/His2 (C2H2) upregulated in the tolerant line and downregulated in the susceptible line under AG stress (Supplementary Figure 2i). These findings suggest that some members of C2H2-type zinc finger proteins (ZFPs) are important regulators of ROS signaling in rice under abiotic stress. The C2H2-type ZFPs have 176 members in Arabidopsis and 1,889 members in rice, and it constitutes one of the largest families of transcriptional regulators in plants [136, 137]. C2H2-type ZFPs are important components in regulating plant growth, development, hormone responses, and tolerance to biotic and abiotic stresses [136, 138, 139]. In rice, several members of C2H2-type ZFPs, such as ZFP182, ZFP245, ZFP252, and ZFP179, are involved in the responses of rice to drought, salinity, and oxidative stresses [140,141,142,143,144]. Many genes regulated by different transcription factors were previously reported to be induced under cold stress. Kanneganti et al. (2008) reported the expression of OsiSAP8 (Os06g0612800) under cold stress . Kong et al. 2006 have identified the expression of OsDOS (Os01g0192000) under drought stress . In our study, we found Os01g0192000 belonging to zinc finger protein upregulated in both genotypes under cold stress with higher log2FC in the susceptible line (3.35) than in the tolerant line (1.89). This shows that this gene may be involved not only in drought stress but also in cold stress.
NAC TF genes have been previously reported to encode a set of plant-specific TFs involved in responses to biotic and abiotic stress [147,148,149]. The expression of NAC genes was found to be upregulated during cold stress , but in our study, we found both upregulated and downregulated NAC genes in both genotypes. This shows NAC genes have diverse roles in cold stress responses in plants. Many genes related to bHLH protein and bHLH transcription factors have also been discovered to be differentially regulated in our study. ICE1, a bHLH protein in Arabidiopsis, has been reported to specifically bind to the CBF3 promoter region and promote the expression of CBF3, which leads to improved tolerance to cold stress . Differential regulation of Os030741100, a bHLH gene induced by cold stress  has also been reported. In this study, we identified more DEGs related to abiotic stress tolerance mechanisms in the stress-tolerant line than in the susceptible line. Functional analysis of the most promising candidate genes identified for stress tolerance in our study can be performed by both gain and loss of function analysis in future studies.
GO analysis of DEGs under different stresses
In order to identify the GO terms enriched by DEGs of each treatment, GO analysis was carried out. Considering the FDR of 0.05 and five minimum number of mapping entries per GO, a total of 232, 306, and 300 GOs were enriched for all DEGs identified for AG stress, cold stress and combined AG and cold stress (Supplementary Table 8). A total of 106, 355 and 331 GOs were enriched for upregulated DEGs of AG stress, cold stress, and combined AG and cold stress. A total of 280, 218 and 242 GOs were enriched for downregulated DEGs of AG stress, cold stress, and combined AG and cold stress. Pairwise comparison of enriched GOs across different stress was performed to find common GOs. A total of 65 and 92 GOs were discovered to be shared between AG stress and cold stress from upregulated and downregulated genes, respectively. A total of 41 and 80 GOs were identified to be shared between AG stress and combined AG and cold stress from upregulated and downregulated DEGs, respectively. A total of 219 and 147 GOs were shared between cold stress and combined AG and cold stress from upregulated and downregulated DEGs, respectively. On the other hand, a total of 108, 37, and 108 GOs were uniquely enriched for upregulated genes in combined stress, AG stress, and cold stress, and a total of 78, 171 and 42 GOs were uniquely enriched for down-regulated genes, respectively. Some of these uniquely enriched GOs from upregulated DEGs of one of these stresses were found to be enriched by down-regulated DEGs of other stresses or vice versa.
The top 10 GO terms of each biological process, molecular function, and cellular component main categories with the highest enrichment scores for each treatment are shown in Fig. 2. Among the top 10 GOs, the common GOs between AG stress, cold stress and combined AG and cold stress, were cellular process, metabolic process, cellular metabolic process, primary metabolic process, catalytic activity, cell part, and cell membrane-bounded organelle and intracellular membrane-bounded organelle. Between AG stress and cold stress, three GO terms were found to be common among top ten GO. i.e., binding, membrane-bounded organelle, and intracellular membrane-bounded organelle. Among cold stress and combined AG and cold stress, 17 GO terms including cellular macromolecule metabolic process, biosynthetic process, cellular biosynthetic process, macromolecule metabolic process, purine nucleotide binding, purine ribonucleotide binding, ribonucleotide binding, transferase activity, purine nucleoside binding, nucleoside binding, adenyl nucleotide binding, intracellular, intracellular part, intracellular organelle, organelle, cytoplasm, and macromolecular complex were found to be common among top 30 GO belong to biological process, cellular component, and molecular function.
The pairwise comparison of enriched GOs of upregulated genes of cold stress and combined AG and cold stress showed enrichment of nitrogen compound metabolic process, carbohydrate metabolic process, response to stress, cellular protein metabolic process, protein metabolic process, cellular ketone metabolic process, cellular nitrogen compound metabolic process, cellular amino acid metabolic process, and lipid metabolic process. Among 147 common GO between cold stress and combined AG and cold stress from downregulated DEGs, post-translational protein modification, phosphate metabolic process, protein amino acid phosphorylation, protein serine/threonine kinase activity, phospholipid binding, and lipid binding were significantly enriched. The GO terms such as transcription, response to stress, nitrogen compound metabolic process, cellular biosynthetic process, RNA metabolic process, carbohydrate metabolic process, photosynthesis, ATP binding, carboxylic acid binding, and cell wall were commonly enriched among upregulated DEGs AG stress and cold stress. Similarly, the GO terms such as carbohydrate metabolic process, alcohol metabolic process, protein modification process, phosphorylation, cellular nitrogen compound metabolic process, lipid binding, lipid localization, lipid transport, cellular ketone metabolic process, glycogen metabolic process, tryptophan metabolic process were commonly enriched among downregulated DEGs AG stress and cold stress. Among AG stress and combined AG and cold stress, the commonly enriched GOs from upregulated genes included nitrogen compound metabolic process, carbohydrate metabolic process, purine nucleotide binding, purine ribonucleotide binding, and response to stimulus, and the commonly enriched GOs from downregulated genes included carbohydrate metabolic process, alcohol metabolic process, protein modification process, phosphate metabolic process, cellular nitrogen compound metabolic process, and lipid binding.
The common GO enriched by upregulated DEGs in AG condition and downregulated DEGs in cold condition included trehalose metabolic process, trehalose biosynthetic process, glycoside metabolic process, disaccharide biosynthetic process, while the GO enriched by downregulated DEGs in AG condition and upregulated DEGs in cold condition included carbohydrate catabolic process, alcohol catabolic process, glucose catabolic process, fructose-bisphosphate aldolase activity. The GO terms enriched by upregulated DEGs in cold conditions and downregulated DEGs in combined AG and cold stress were chromatin assembly or disassembly, nucleosome organization, chromosome organization, and DNA packaging. The GO terms enriched by upregulated DEGs in cold conditions and downregulated DEGs in AG conditions were glycolysis, carbohydrate catabolic process, alcohol catabolic process, phosphotransferase activity, and phosphate group as acceptor. The GO terms enriched by upregulated DEGs in the combined condition of AG and cold stress and downregulated DEGs in AG condition were polysaccharide biosynthetic process, fatty acid metabolic process, lipid biosynthetic process, and cellular lipid metabolic process. The GO terms enriched by upregulated DEGs in AG and downregulated DEGs in the combined condition of AG and cold stress were exocyst and cell cortex. The GO terms enriched by downregulated DEGs in AG and upregulated DEGs in the combined condition of AG and cold stress were cellular lipid metabolic process, fatty acid metabolic process, polysaccharide biosynthetic process, dicarboxylic acid metabolic process, and oxidoreductase activity.
Glycosyl hydrolases family 16, Os06g0697000, was upregulated in both AG stress and cold stress but downregulated under combined AG and cold stress. In contrast, the opposite trend was shown by another glycosyl hydrolases family 16, Os06g0696566. Glycoside hydrolases were identified to be involved in the hydrolysis of glycosidic bonds in complex sugars such as cellulose and starch [152, 153]. Glycosyl hydrolase 14 family gene, Os10g0565200, encoding beta-amylase, was reported to be significantly upregulated under cold stress . Another glycosyl hydrolase gene, Os06g0675700, was significantly upregulated in the cold-tolerant line than in the cold-sensitive line . Transcriptomics analysis of anther of contrasting genotypes revealed Os3BGlu7  and GIF1 , encoding glycosyl hydrolase family 1 and glycosyl hydrolase family 32, respectively, were significantly upregulated in the cold susceptible genotype, LJ11 compared with cold-tolerant genotype, LJ25 under cold treatment . The glycosidase hydrolase genes, Os11g0701400, Os11g0701600, Os11g0701900, and Os11g0702200 were significantly upregulated in the tolerant line under submergence . Another member of the glycosyl hydrolase 14 family gene, Os10g0565200, was significantly upregulated under cold stress conditions , suggesting the potential involvement of the enzymes of the glycosyl hydrolase family on cold and AG tolerance.
Among the three common DEGs related to cellulose synthesis and representing cellulose synthase-like family A and family E, CSLA1(Os02g0192500) was upregulated under cold stress and combined AG and cold stress but downregulated under AG stress, whereas CSLE1 (Os09g0478100) and CSLE6 (Os09g0478300) were downregulated under all three stresses. High expression of cellulose synthase genes is reported to increase tolerance to biotic  and abiotic  stresses in plants. Endler et al. (2015) demonstrated that maintaining microtubule organization for cellulose synthase positioning and continuous cellulose deposition is critical for biomass production during salt stress in Arabidosis . Sperotto et al. (2018) used cellulose staining and showed higher cellulose deposition in cell walls during cold treatment compared to normal temperatures . This group also reported higher expression of cell-wall-related genes involved in cell-wall biosynthesis in the cold-tolerant line than in the cold-sensitive line, indicating the importance of cellulose synthase and cellulose-synthase like proteins on maintaining cellular integrity during cold stress. Kong et al. (2010) reported downregulation of proteins related to cell wall metabolism and structural modification, such as B-glucanase, B-glucosidase, B-galactosidase in wheat, suggesting the correlation of downregulated proteins with inhibition of cell wall elongation under flooding stress suppressing the growth of wheat . Similar trends were observed in woody plants such as gray poplar, where the expression of genes encoding cellulose synthases was strongly down-regulated (up to 44-fold) in roots under flooded condition .
Two common DEGs across all stresses were upregulated in cold and combined AG and cold stress and were downregulated under AG stress: one related to phospholipid synthesis which is cyclopropane-fatty-acyl-phospholipid (Os06g0574100) and the other, glyceraldehyde-3-phosphate dehydrogenase (Os02g0171100), was involved in glycolysis. Previous studies have demonstrated the upregulation of phosphometabolic processes under cold stress in Lotus japonicus, and Arabidopsis thaliana [161,162,163]. Changes in membrane phospholipids were reported to help adapt to low-temperature environments .
KEGG and MapMan analysis of DEGs
The KEGG pathway enrichment analysis was performed to determine the function of the differentially expressed genes (DEGs) [165, 166]. Among AG stress-responsive DEGs, 5571 DEGs were mapped to 102 biological pathways, 26 cellular components, and 71 molecular functions. Likewise, 7206 cold stress-responsive DEGs were mapped to 356 biological pathways, 89 cellular components, and 81 molecular functions, while 13279 combined stress-responsive DEGs were mapped to 500 biological pathways, 131 cellular components, and 112 molecular functions. The top 20 pathways enriched by AG, cold, and combined stress-responsive DEGs are shown in Table 2.
As a complementary approach to the GO and KEGG analysis, MapMan analysis was also employed to categorize DEGs in bins and molecular pathways. Overall, AG stress-responsive DEGs were mapped to 84 bins (Fig. 3a). Secondary metabolism, minor aldehyde metabolism (trehalose), lipid metabolism, fermentation, and amino acid metabolism were highly enriched by DEGs under AG. Cold stress-responsive DEGs were mapped to 98 bins (Fig. 3b). Cell wall modification, light reaction, minor aldehyde metabolism (raffinose family), minor aldehyde metabolism (starch degradation family), and nucleotide metabolism were highly enriched by DEGs under cold stress. DEGs from combined AG and cold stress were mapped to 102 bins (Fig. 3c). Cell wall modification, mitochondrial electron transport (ATP synthesis), minor aldehyde metabolism (raffinose family), and lipid metabolism were highly enriched by DEGs under combined stress conditions. Cell wall degradation, cell wall modification, and minor aldehyde metabolism (raffinose family) were commonly enriched from DEGs from cold stress and combined AG and cold stress. Minor CHO metabolism (trehalose) and Calvin cycle were commonly enriched from DEGs under AG and cold stress. Amino acid metabolism and secondary metabolism were commonly enriched from DEGs from AG stress and combined AG and cold stress. The molecular pathways related to amino acid metabolism (synthesis of tryptophan, methionine), cell wall (cellulose synthesis), cell wall modification, fermentation, glycolysis, lipid metabolism (steroids, squalene, glycolipid synthesis), minor CHO metabolism (raffinose family, trehalose), N-metabolism (ammonia metabolism, nitrate metabolism), light reaction, and photorespiration were commonly enriched from DEGs from all three stresses. The number of DEGs differed among the common pathways identified with MapMan across different stresses. Additionally, the DEGs upregulated in one stress were downregulated in another stress for the same pathway. Among the genes related to cell wall modification, six DEGs related to expansin precursor and glycosyl hydrolase family were common across all stresses. Expansin precursors, Os03g0102500, Os04g0228400, Os10g0555600, and Os10g0556100, were upregulated under cold stress and combined cold and AG stress. In contrast, they were downregulated under AG stress. Expansins facilitate cell wall loosening by disrupting the hydrogen bonds between cellulose and xyloglucan polymers . In a previous study, expansins were up-regulated at one or more time points of cold treatment in the tolerant line, NC567 and downregulated in susceptible line, Taiyan8 . The upregulation of expansins in the roots of the drought-tolerant line, Nootripathu compared to drought susceptible line, IR20, was also reported . Significant downregulation of the B-expansin gene (Os03g0106500), which regulates cell wall dynamics and cell size in rice, was also observed . Separate MapMan analysis of each genotype under all three different types of stresses showed large number of enriched genes related to phenylpropanoids, phenomics, terpenes and flavonoids in both Darij and 4610 under combined AG and cold stress followed by AG stress and cold stress (Supplementary Figures 3, 4). A large number of upregulated genes were enriched in both Darij and 4610 under combined AG and cold stress whereas a large number of downregulated genes were enriched in both Darij and 4610 under AG and cold stress. For the secondary metabolism category, the common pathways enriched in all three stresses were phenylpropanoids, simple phenols, mevalonate (MVA), lignin biosynthesis, non-mevalonate pathway, terpenoids, and carotenoids (Fig. 4). The glucosinolates category was uniquely enriched for DEGs under cold stress (Fig. 4b) and combined AG and cold stress (Fig. 4c). Anthocyanins and tocopherol biosynthesis were enriched for AG stress (Fig. 4a) and combined AG and cold stress. Tyrosine was uniquely enriched for combined AG and cold stress. Separate MapMan analysis of secondary metabolism of Darij and 4610 under three different stress conditions showed large number of enriched genes related to lignins and lignans, and MVA pathway in combined stress followed by AG stress and cold stress (Supplementary Figs. 5, 6).
Protein-protein interactions network of DEGs
A protein-protein interaction (PPI) network demonstrates various complex biological systems into a set of nodes and edges connecting the nodes. A total of 992 common DEGs identified among AG stress, cold stress and combined AG and cold stress were imported for protein-protein interactions. A high confidence score of 0.7 was used to construct the PPI. A hub gene was determined if the gene was connected to more than five genes. Three hub genes were identified in the data set (Fig. 5). Hub gene Os08g0176800 (similar to mRNA-associated protein mrnp 41), was found to be connected to 10 genes, including Os01g0159300 (zinc finger C3HC4 type domain containing protein), Os07g0274100 (calcineurin-related phosphoesterase-like protein), Os08g0505900 (putative DNA-damage-repair/toleration protein), and Os01g0159300 (zinc finger C3HC4 type domain-containing protein) (Fig. 6a). Another hub gene, Os11g0454200 (dehydrin) was found to be connected to seven other genes, including Os01g0705200 (late embryogenesis abundant protein 14), Os04g0589800 (late embryogenesis abundant, LEA), Os05g0542500 (LEA), and Os11g0454300 (LEA) (Fig. 6b). The last hub gene, Os10g0505900 (expressed protein) was found to be connected to seven genes, including Os01g0705200 (LEA protein 14), Os02g0158900 (CBS domain containing membrane protein, putative, expressed), Os05g0542500 (LEA), and Os11g0454200 (dehydrin) (Fig. 6c). The LEA protein family are unstructured polypeptides without a well-defined three-dimensional structure and have been reported to be significantly induced by abiotic stresses such as cold, drought, and salinity [171, 172]. The LEA protein family has been reported to be enhanced in drought tolerance by stabilizing cell membranes and protecting proteins . The LEA protein family represents stress-induced hyper-hydrophilic proteins that accumulate in response to cellular dehydration. Previous research has postulated a positive correlation between the expression of LEAs and abiotic stress tolerance in plants [174, 175]. In rice, LEA genes of LEA_2, LEA_3, and dehydrin exhibited a strong response to osmotic stress (PEG), salt, and ABA . OLE18, from the oleosin family, may have a structural role in stabilizing the lipid body during desiccation of the seed by preventing coalescence of the oil or providing recognition signals for specific lipase anchorage in lipolysis during seedling growth . RAB16B and RAB21, from the plant dehydrin family, were known to be water stress-inducible proteins [178, 179].
This study presents the findings of whole genome transcriptome profiles of two germinating rice genotypes that exhibit a contrasting response to flooding stress, cold stress, and combined flooding and cold stress. Our study highlights that the genes involved in carbohydrate metabolic process, glucosyltransferase activity, regulation of nitrogen compound metabolic process, protein metabolic process, lipid metabolic process, cellular nitrogen compound biosynthetic process, lipid biosynthetic process, and microtubule-based process were commonly enriched across all the stresses. Notably, we found that a significant number of DEGs have contributed to phenotypic differences in different stress conditions between the tolerant and susceptible genotypes. This study has provided the basis for identifying potential candidate genes involved in flooding, cold, and combined stress tolerance in germinating rice seed and possible use of the information for crop improvement.
Materials and methods
Plant materials and stress treatment
Two rice genotypes selected out of 257 accessions used for GWAS analysis and showing contrasting responses to anaerobic germination (AG) stress and cold stress  were selected. Darij was found to be tolerant, and 4610 was found to be susceptible, to both AG stress and cold stress during germination. All the seeds used in this study were freshly harvested and air-dried at 37°C for 5 days and then stored at 4°C prior to the experiment. In this experiment, we used seeds collected from single plants of each accession to maintain uniformity. Considering seed age effects on germination and vigor, we only used seeds that were about nine months old. The dormancy of seeds was broken by oven drying the seeds at 50°C for 5 days. The seeds were dehulled and surface sterilized with 70% ethanol for 5 mins, followed by rinsing with autoclaved distilled water three times. A total of ten replicates with 30 seeds/rep were used for each control and treatment in a completely randomized design for the phenotyping. The sterilized seeds were then germinated for 4 days on moist filter paper placed inside Petri dishes for phenotyping under control conditions. For submergence treatments, sterilized seeds were germinated in 8 cm deep water in glass beakers for 4 days.
For phenotyping under cold stress treatment, seeds were germinated in moist filter paper placed inside Petri dishes for 28 days in a refrigerator maintained at 13°C. A set of phenotypic traits, including shoot length under normal growth conditions (30°C) after 7 days after seeding (DAS), shoot length after 28 days of cold exposure (13°C), plumule recovery rate after 4 days of bringing the stressed seedlings back to normal temperature (30°C) was measured. Only seeds having coleoptile length half of the seed size were considered germinated . For combined AG and cold stress treatment, sterilized seeds were germinated in 8 cm deep water in glass beakers at 13°C temperature for 28 and 35 days.
Extraction of RNA and library construction
A total of four replicates with 30 seeds/rep were used for each control and treatment for the RNAseq samples. For preparing the samples for RNA extraction under cold stress, the surface-sterilized seeds were germinated in a growth chamber maintained at 13°C for 7 days in dark conditions. For the control samples, both the contrasting lines were germinated for 2 days in a growth chamber maintained at 30°C in dark conditions. By doing so, it is expected that the treated samples and the controls will have a similar physiological state since the growth of coleoptiles under cold stress are generally very slow. The germinating seeds with coleoptile length half the size of the seed were considered germinated. Germinating embryos from four biological replicates of each line grown under cold stress (7 DAS) and control conditions (2 DAS) were collected and immediately frozen in liquid nitrogen. Samples were then maintained at -80°C until RNA extraction. For samples under AG conditions, sterilized seeds of two contrasting lines were germinated in 8 cm deep water in a glass beaker for 4 days. For the control condition, sterilized seeds were germinated for 4 days on moist filter paper placed inside Petri dishes. Coleoptiles from four biological replicates of each line grown under AG stress (4 DAS) and normal condition (4 DAS) were collected and immediately frozen in liquid nitrogen and then maintained at -80°C until RNA extraction. For preparing samples for the RNA extraction under combination of AG stress and cold stress, sterilized seeds were germinated in 8 cm deep water in glass beakers at 13°C temperature for 7 days. For the control samples, both the contrasting lines were germinated for 2 days in a growth chamber maintained at 30°C in dark conditions. Germinating embryos from four biological replicates of each line grown under a combination of AG and cold stress (7 DAS) and under control (2 DAS) were collected and immediately frozen in liquid nitrogen, and the samples were then preserved at -80°C until RNA extraction. Any remaining embryos were included in the samples.
Total RNA was extracted using TRIzol reagent (Invitrogen, MA, USA) and QIAgen RNeasy Plant Mini Kit (Qiagen). RNA quality and concentration were checked using spectrophotometer NanoDrop ND-1000, and the integrity was checked using 1% agarose gel electrophoresis. Poly-A RNA containing mRNA was purified using poly-T oligo-attached magnetic beads and fragmented, and complementary DNA (cDNA) was synthesized using random hexamer primers, followed by purification, end-repairing, poly-A tailing, and adapter ligation. TruSeq Stranded RNA-seq libraries were prepared at Texas A&M AgriLife Genomics and Bioinformatics Service (TxGen; College Station, TX, USA) as per Standard Operating Protocol (Illumina). The quality of RNA libraries was further evaluated using Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) using the kit Agilent DNA 1000. The libraries were run on multiple lanes of a HiSeq 4000 to generate 75 nt pair-end reads.
RNA-seq analysis and differential expressed genes (DEGs)
FastQC was used for the analysis of the read quality and its visualization . The low-quality bases and library adapters were removed from each library using Trimmomatic version 0.36 . Reads were mapped against the reference genome of Oryza sativa cv. Nipponbare (IRGSP build 1.0 RAP-DB) [182, 183]. The mapping and alignment of the reads on the rice genome were done using HISAT2 version 2.1.0 . Expression levels of each gene were quantified by normalizing total gene counts with the effective library size. DESeq2 v1.26 was used to test the pairwise differential expression analysis . Genes with a p-value less than 0.05 were considered to be differentially expressed. We used the multi-factorial linear modeling and tested two null hypotheses of effects on gene expression: (1) whether stress treatment has a significant effect on the expression of each gene, and (2) whether gene expression was affected by stress treatment in a genotype-dependent manner. We fitted for models with our experimental data: (1) FMtrt: Y = τ + ε; (2) FMadd: Y = τ + γ + ε, and; (2) FMfull: Y = τ + γ + τ:γ + ε. In each model, Y is the expression value of each gene, τ is the effect of stress treatment, γ is the effect of different genotype, and ε is the random error. Comparing FMtrt to FMadd separately, we tested whether the expression of each gene was regulated by stress, and whether there was a significant genotypic effect under stress. Comparisons of FMfull and FMadd allowed us to test whether gene expression was affected by stress in a genotype-dependent manner. In all cases, expression values for genes were standardized using the expression z=(x−x̅)/stdev for cross-genotype comparisons.
Gene ontology and enrichment analysis
We used the web-based tool AGRIGO v2.0 to perform Gene Ontology (GO) enrichment analysis on DEGs , while the Fisher’s exact test and Benjamini-Hochberg’s false discovery rate (FDR) adjustment were used to control for multiple comparisons . The GO term with FDR < 0.05 was regarded as significantly enriched. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was performed using ShinyGO v0.75 . MapMan 3.5.1 software was applied to map DEGs to O. sativa pathways to determine genes involved in specific pathways .
Protein-protein interactions network of DEGs
Common DEGs detected in AG stress, cold stress, and a combination of AG and cold stress were used to construct protein-protein interactions (PPIs) network. The STRING v11.5 database was used to detect both validated and predicted protein-protein interactions . The resulting interactions were used to construct the PPIs network that was further analyzed and visualized using Cytoscape v3.7.0 . PPIs network was constructed to identify hub genes (the nodes with the highest degree) involved in different stresses.
Availability of data and materials
All the data supporting the results of this article are provided within the article or in the additional files.
Ismail AM. Flooding and submergence tolerance. In: Kole C, editor. Genomics and Breeding for Climate-Resilient Crops: Vol. 2 Target Traits. Berlin, Heidelberg: Springer; 2013. p. 269–90.
Septiningsih EM, Pamplona AM, Sanchez DL, Neeraja CN, Vergara GV, Heuer S, Ismail AM, Mackill DJ. Development of submergence-tolerant rice cultivars: the Sub1 locus and beyond. Ann Bot. 2008;103(2):151–60.
Gonzaga ZJC, Carandang J, Singh A, Collard BC, Thomson MJ, Septiningsih EM. Mapping QTLs for submergence tolerance in rice using a population fixed for SUB1A tolerant allele. Mol Breed. 2017;37(4):1–10.
Gonzaga ZJC, Carandang J, Sanchez DL, Mackill DJ, Septiningsih EM. Mapping additional QTLs from FR13A to increase submergence tolerance in rice beyond SUB1. Euphytica. 2016;209(3):627–36.
Singh A, Carandang J, Gonzaga ZJC, Collard BC, Ismail AM, Septiningsih EM. Identification of QTLs for yield and agronomic traits in rice under stagnant flooding conditions. Rice. 2017;10(1):1–18.
Septiningsih EM, Mackill DJ. Genetics and breeding of flooding tolerance in rice. In: Sasaki T, Ashikari M, editors. Rice genomics, genetics and breeding. Singapore: Springer; 2018. p. 275–95.
Kato Y, Collard BC, Septiningsih EM, Ismail AM. Increasing flooding tolerance in rice: combining tolerance of submergence and of stagnant flooding. Ann Bot. 2019;124(7):1199–209.
Magneschi L, Perata P. Rice germination and seedling growth in the absence of oxygen. Ann Bot. 2008;103(2):181–96.
Narsai R, Edwards JM, Roberts TH, Whelan J, Joss GH, Atwell BJ. Mechanisms of growth and patterns of gene expression in oxygen-deprived rice coleoptiles. Plant J. 2015;82(1):25–40.
Septiningsih EM, Ignacio JCI, Sendon PM, Sanchez DL, Ismail AM, Mackill DJ. QTL mapping and confirmation for tolerance of anaerobic conditions during germination derived from the rice landrace Ma-Zhan Red. Theor Appl Genet. 2013;126(5):1357–66.
Hsu S-K, Tung C-W. Genetic mapping of anaerobic germination-associated QTLs controlling coleoptile elongation in rice. Rice. 2015;8:1–12.
Zhang M, Lu Q, Wu W, Niu X, Wang C, Feng Y, Xu Q, Wang S, Yuan X, Yu H. Association mapping reveals novel genetic loci contributing to flooding tolerance during germination in Indica rice. Front Plant Sci. 2017;8:678.
Kim S-M, Reinke RF. Identification of QTLs for tolerance to hypoxia during germination in rice. Euphytica. 2018;214(9):160.
Baltazar MD, Ignacio JCI, Thomson MJ, Ismail AM, Mendioro MS, Septiningsih EM. QTL mapping for tolerance to anaerobic germination in rice from IR64 and the aus landrace Kharsu 80A. Breed Sci. 2019;69(2):227–33.
Ghosal S, Quilloy FA, Casal C, Septiningsih EM, Mendioro MS, Dixit S. Trait-based mapping to identify the genetic factors underlying anaerobic germination of rice: Phenotyping, GXE, and QTL mapping. BMC Genet. 2020;21(1):1–13.
Ghosal S, Casal C, Quilloy FA, Septiningsih EM, Mendioro MS, Dixit S. Deciphering genetics underlying stable anaerobic germination in rice: phenotyping, QTL identification, and interaction analysis. Rice. 2019;12(1):1–15.
Yang J, Sun K, Li D, Luo L, Liu Y, Huang M, Yang G, Liu H, Wang H, Chen Z. Identification of stable QTLs and candidate genes involved in anaerobic germination tolerance in rice via high-density genetic mapping and RNA-Seq. BMC Genom. 2019;20(1):1–15.
Jeong JM, Cho YC, Jeong JU, Mo YJ, Kim CS, Kim WJ, Baek MK, Kim SM. QTL mapping and effect confirmation for anaerobic germination tolerance derived from the japonica weedy rice landrace PBR. Plant Breed. 2020;139(1):83–92.
Ignacio JCI, Zaidem M, Casal C Jr, Dixit S, Kretzschmar T, Samaniego JM, Mendioro MS, Weigel D, Septiningsih EM. Genetic mapping by sequencing more precisely detects loci responsible for anaerobic germination tolerance in rice. Plants. 2021;10(4):705.
Tnani H, Chebotarov D, Thapa R, Ignacio JCI, Israel WK, Quilloy FA, Dixit S, Septiningsih EM, Kretzschmar T. Enriched-GWAS and Transcriptome Analysis to Refine and Characterize a Major QTL for Anaerobic Germination Tolerance in Rice. Int J Mol Sci. 2021;22(9):4445.
Su L, Yang J, Li D, Peng Z, Xia A, Yang M, Luo L, Huang C, Wang J, Wang H. Dynamic genome-wide association analysis and identification of candidate genes involved in anaerobic germination tolerance in rice. Rice. 2021;14:1–22.
Thapa R, Tabien RE, Thomson MJ, Septiningsih EM: Genetic factors underlying anaerobic germination in rice: Genome‐wide association study and transcriptomic analysis. Plant Genome 2022:e20261.
Kretzschmar T, Pelayo MAF, Trijatmiko KR, Gabunada LFM, Alam R, Jimenez R, Mendioro MS, Slamet-Loedin IH, Sreenivasulu N, Bailey-Serres J. A trehalose-6-phosphate phosphatase enhances anaerobic germination tolerance in rice. Nat Plants. 2015;1(9):15124.
Guglielminetti L, Perata P, Alpi A. Effect of anoxia on carbohydrate metabolism in rice seedlings. Plant Physiol. 1995;108(2):735–41.
Perata P, Guglielminetti L, Alpi A. Mobilization of endosperm reserves in cereal seeds under anoxia. Ann Bot. 1997;79(suppl_1):49–56.
Ismail AM, Ella ES, Vergara GV, Mackill DJ. Mechanisms associated with tolerance to flooding during germination and early seedling growth in rice (Oryza sativa). Ann Bot. 2008;103(2):197–209.
Pompeiano A, Fanucchi F, Guglielminetti L. Amylolytic activity and carbohydrate levels in relation to coleoptile anoxic elongation in Oryza sativa genotypes. J Plant Res. 2013;126(6):787–94.
Setter T, Ella E, Valdez A. Relationship between coleoptile elongation and alcoholic fermentation in rice exposed to anoxia. II. Cultivar differences. Ann Bot. 1994;74(3):273–9.
Gibbs J, Morrell S, Valdez A, Setter T, Greenway H. Regulation of alcoholic fermentation in coleoptiles of two rice cultivars differing in tolerance to anoxia. J Exp Bot. 2000;51(345):785–96.
Magneschi L, Kudahettige R, Alpi A, Perata P. Comparative analysis of anoxic coleoptile elongation in rice varieties: relationship between coleoptile length and carbohydrate levels, fermentative metabolism and anaerobic gene expression. Plant Biol. 2009;11(4):561–73.
Atwell BJ, Greenway H, Colmer TD. Efficient use of energy in anoxia-tolerant plants with focus on germinating rice seedlings. New Phytologist. 2015;206(1):36–56.
Edwards JM, Roberts TH, Atwell BJ. Quantifying ATP turnover in anoxic coleoptiles of rice (Oryza sativa) demonstrates preferential allocation of energy to protein synthesis. J Exp Bot. 2012;63(12):4389–402.
Lasanthi-Kudahettige R, Magneschi L, Loreti E, Gonzali S, Licausi F, Novi G, Beretta O, Vitulli F, Alpi A, Perata P. Transcript profiling of the anoxic rice coleoptile. Plant Physiol. 2007;144(1):218–31.
Huang S, Taylor NL, Narsai R, Eubel H, Whelan J, Millar AH. Experimental analysis of the rice mitochondrial proteome, its biogenesis, and heterogeneity. Plant Physiol. 2009;149(2):719–34.
Shingaki-Wells RN, Huang S, Taylor NL, Carroll AJ, Zhou W, Millar AH. Differential molecular responses of rice and wheat coleoptiles to anoxia reveal novel metabolic adaptations in amino acid metabolism for tissue tolerance. Plant Physiol. 2011;156(4):1706–24.
Jackson MB. Ethylene and responses of plants to soil waterlogging and submergence. Ann Rev Plant Physiol. 1985;36(1):145–74.
Crawford R. Oxygen availability as an ecological limit to plant distribution. In: Begon M, Fitter AH, editors. Advances in ecological research. Vol. 23. Academic Press; 1992. p. 93–185.
Perata P, Alpi A. Plant responses to anaerobiosis. Plant Sci. 1993;93(1–2):1–17.
Armstrong W, Brändle R, Jackson MB. Mechanisms of flood tolerance in plants. Acta Bot Neerl. 1994;43(4):307–58.
Drew MC, He C-J, Morgan PW. Programmed cell death and aerenchyma formation in roots. Trends Plant Sci. 2000;5(3):123–7.
Sauter M. Rice in deep water:" How to take heed against a sea of troubles". Naturwissenschaften. 2000;87(7):289–303.
Gibbs J, Greenway H. Mechanisms of anoxia tolerance in plants. I. Growth, survival and anaerobic catabolism. Funct Plant Biol. 2003;30(1):1–47.
Voesenek L, Colmer T, Pierik R, Millenaar F, Peeters A. How plants cope with complete submergence. New Phytol. 2006;170(2):213–26.
Baxter-Burrell A, Yang Z, Springer PS, Bailey-Serres J. RopGAP4-dependent Rop GTPase rheostat control of Arabidopsis oxygen deprivation tolerance. Science. 2002;296(5575):2026–8.
Hoeren FU, Dolferus R, Wu Y, Peacock WJ, Dennis ES. Evidence for a role for AtMYB2 in the induction of the Arabidopsis alcohol dehydrogenase gene (ADH1) by low oxygen. Genetics. 1998;149(2):479–90.
Liu F, VanToai T, Moy LP, Bock G, Linford LD, Quackenbush J. Global transcription profiling reveals comprehensive insights into hypoxic response in Arabidopsis. Plant Physiol. 2005;137(3):1115–29.
Dordas C, Hasinoff BB, Rivoal J, Hill RD. Class-1 hemoglobins, nitrate and NO levels in anoxic maize cell-suspension cultures. Planta. 2004;219(1):66–72.
Vriezen WH, Hulzink R, Mariani C, Voesenek LA. 1-Aminocyclopropane-1-carboxylate oxidase activity limits ethylene biosynthesis in Rumex palustrisduring submergence. Plant Physiol. 1999;121(1):189–96.
Mattana M, Coraggio I, Bertani A, Reggiani R. Expression of the enzymes of nitrate reduction during the anaerobic germination of rice. Plant Physiol. 1994;106(4):1605–8.
Saab IN, Sachs MM. A flooding-induced xyloglucan endo-transglycosylase homolog in maize is responsive to ethylene and associated with aerenchyma. Plant Physiol. 1996;112(1):385–91.
Sales E, Viruel J, Domingo C, Marqués L. Genome wide association analysis of cold tolerance at germination in temperate japonica rice (Oryza sativa L.) varieties. PloS one. 2017;12(8):e0183416.
Shakiba E, Edwards JD, Jodari F, Duke SE, Baldo AM, Korniliev P, McCouch SR, Eizenga GC. Genetic architecture of cold tolerance in rice (Oryza sativa) determined through high resolution genome-wide analysis. PloS one. 2017;12(3):e0172133.
Thapa R, Tabien RE, Thomson MJ, Septiningsih EM. Genome-Wide Association Mapping to Identify Genetic Loci for Cold Tolerance and Cold Recovery During Germination in Rice. Front Genet. 2020;11:22.
Ma Y, Dai X, Xu Y, Luo W, Zheng X, Zeng D, Pan Y, Lin X, Liu H, Zhang D. COLD1 confers chilling tolerance in rice. Cell. 2015;160(6):1209–21.
Liu C, Ou S, Mao B, Tang J, Wang W, Wang H, Cao S, Schläppi MR, Zhao B, Xiao G. Early selection of bZIP73 facilitated adaptation of japonica rice to cold climates. Nat Commun. 2018;9(1):1–12.
Zhang Z, Li J, Pan Y, Li J, Shi H, Zeng Y, Guo H, Yang S, Zheng W, Yu J. Natural variation in CTB4a enhances rice adaptation to cold habitats. Nature Commun. 2017;8(1):1–13.
Rabbani MA, Maruyama K, Abe H, Khan MA, Katsura K, Ito Y, Yoshiwara K, Seki M, Shinozaki K, Yamaguchi-Shinozaki K. Monitoring expression profiles of rice genes under cold, drought, and high-salinity stresses and abscisic acid application using cDNA microarray and RNA gel-blot analyses. Plant Physiol. 2003;133(4):1755–67.
Maruyama K, Todaka D, Mizoi J, Yoshida T, Kidokoro S, Matsukura S, Takasaki H, Sakurai T, Yamamoto YY, Yoshiwara K. Identification of cis-acting promoter elements in cold-and dehydration-induced transcriptional pathways in Arabidopsis, rice, and soybean. DNA Res. 2011;19(1):37–49.
Takasaki H, Maruyama K, Kidokoro S, Ito Y, Fujita Y, Shinozaki K, Yamaguchi-Shinozaki K, Nakashima K. The abiotic stress-responsive NAC-type transcription factor OsNAC5 regulates stress-inducible genes and stress tolerance in rice. Mol Genet Genom. 2010;284(3):173–83.
Lv D-K, Bai X, Li Y, Ding X-D, Ge Y, Cai H, Ji W, Wu N, Zhu Y-M. Profiling of cold-stress-responsive miRNAs in rice by microarrays. Gene. 2010;459(1–2):39–47.
Yun K-Y, Park MR, Mohanty B, Herath V, Xu F, Mauleon R, Wijaya E, Bajic VB, Bruskiewich R. de los Reyes BG: Transcriptional regulatory network triggered by oxidative signals configures the early response mechanisms of japonica rice to chilling stress. BMC Plant Biol. 2010;10(1):16.
Modi A, Vai S, Caramelli D, Lari M. The Illumina Sequencing Protocol and the NovaSeq 6000 System. Methods Mol Biol (Clifton, NJ). 2021;2242:15–42.
Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.
Alpi A, Beevers H. Effects of O2 concentration on rice seedlings. Plant Physiol. 1983;71(1):30–4.
Kawai M, Uchimiya H. Coleoptile senescence in rice (Oryza sativa L.). Ann Bot. 2000;86(2):405–14.
Cruz RPd, Milach SCK. Cold tolerance at the germination stage of rice: methods of evaluation and characterization of genotypes. Sci Agric. 2004;61(1):1–8.
Miro B, Longkumer T, Entila FD, Kohli A, Ismail AM. Rice seed germination underwater: Morpho-physiological responses and the bases of differential expression of alcoholic fermentation enzymes. Front Plant Sci. 1857;2017:8.
Miura K, Kuroki M, Shimizu H, Ando I. Introduction of the long-coleoptile trait to improve the establishment of direct-seeded rice in submerged fields in cool climates. Plant Prod Sci. 2002;5(3):219–23.
Hsu S-K, Tung C-W. Genetic Mapping of Anaerobic Germination-Associated QTLs Controlling Coleoptile Elongation in Rice. Rice. 2015;8(1):38.
da Maia LC, Cadore PR, Benitez LC, Danielowski R, Braga EJ, Fagundes PR, Magalhães AM, de Oliveira AC. Transcriptome profiling of rice seedlings under cold stress. Func Plant Biol. 2016;44(4):419–29.
Pan Y, Liang H, Gao L, Dai G, Chen W, Yang X, Qing D, Gao J, Wu H, Huang J. Transcriptomic profiling of germinating seeds under cold stress and characterization of the cold-tolerant gene LTG5 in rice. BMC Plant Biol. 2020;20(1):1–17.
Drew MC. Oxygen deficiency and root metabolism: injury and acclimation under hypoxia and anoxia. Ann Rev Plant Biol. 1997;48(1):223–50.
Felle HH. Apoplastic pH during low-oxygen stress in barley. Ann Bot. 2006;98(5):1085–93.
Limami AM, Diab H, Lothier J. Nitrogen metabolism in plants under low oxygen stress. Planta. 2014;239(3):531–41.
Bailey-Serres J, Fukao T, Gibbs DJ, Holdsworth MJ, Lee SC, Licausi F, Perata P, Voesenek LA, van Dongen JT. Making sense of low oxygen sensing. Trends Plant Sci. 2012;17(3):129–38.
Bailey-Serres J, Voesenek L. Flooding stress: acclimations and genetic diversity. Annu Rev Plant Biol. 2008;59:313–39.
Allègre A, Silvestre J, Morard P, Kallerhoff J, Pinelli E. Nitrate reductase regulation in tomato roots by exogenous nitrate: a possible role in tolerance to long-term root anoxia. J Exp Bot. 2004;55(408):2625–34.
Morard P, Silvestre J, Lacoste L, Caumes E, Lamaze T. Nitrate uptake and nitrite release by tomato roots in response to anoxia. J Plant Physiol. 2004;161(7):855–65.
Horchani F, Aschi-Smiti S, Brouquisse R. Involvement of nitrate reduction in the tolerance of tomato (Solanum lycopersicum L.) plants to prolonged root hypoxia. Acta Physiol Plant. 2010;32(6):1113–23.
Zhang G-Z, Jin S-H, Jiang X-Y, Dong R-R, Li P, Li Y-J, Hou B-K. Ectopic expression of UGT75D1, a glycosyltransferase preferring indole-3-butyric acid, modulates cotyledon development and stress tolerance in seed germination of Arabidopsis thaliana. Plant Mol Biol. 2016;90(1–2):77–93.
Guglielminetti L, Morita A, Yamaguchi J, Loreti E, Perata P, Alpi A. Differential expression of two fructokinases in Oryza sativa seedlings grown under aerobic and anaerobic conditions. J Plant Res. 2006;119(4):351–6.
Kennedy R, Fox T, Everard J, Rumpho M. Biochemical adaptations to anoxia: potential role of mitochondrial metabolism to flood tolerance in Echinochloa phyllopogon (barnyard grass). In: Plant life under oxygen deprivation Ecology, physiology and biochemistry The Hague: SPB Academic. 1991. p. 217–27.
Fox TC, Kennedy RA, Rumpho ME. Energetics of plant growth under anoxia: metabolic adaptations of Oryza sativa and Echinochloa phyllopogon. Ann Bot. 1994;74(5):445–55.
Vartapetian B, Mazliak P, Lance C. Lipid biosynthesis in rice coleoptiles grown in the presence or in the absence of oxygen. Plant Sci Lett. 1978;13(4):321–8.
Panda D, Barik J. Flooding tolerance in rice: Focus on mechanisms and approaches. Rice Sci. 2021;28(1):43–57.
Nanjo Y, Maruyama K, Yasue H, Yamaguchi-Shinozaki K, Shinozaki K, Komatsu S. Transcriptional responses to flooding stress in roots including hypocotyl of soybean seedlings. Plant molecular biology. 2011;77(1–2):129–44.
Kreuzwieser J, Hauberg J, Howell KA, Carroll A, Rennenberg H, Millar AH, Whelan J. Differential response of gray poplar leaves and roots underpins stress adaptation during hypoxia. Plant Physiol. 2009;149(1):461–73.
Lee Y-H, Kim K-S, Jang Y-S, Hwang J-H, Lee D-H, Choi I-H. Global gene expression responses to waterlogging in leaves of rape seedlings. Plant Cell Rep. 2014;33(2):289–99.
Klok EJ, Wilson IW, Wilson D, Chapman SC, Ewing RM, Somerville SC, Peacock WJ, Dolferus R, Dennis ES. Expression profile analysis of the low-oxygen response in Arabidopsis root cultures. Plant Cell. 2002;14(10):2481–94.
Clement M, Lambert A, Herouart D, Boncompagni E. Identification of new up-regulated genes under drought stress in soybean nodules. Gene. 2008;426(1–2):15–22.
Wrzaczek M, Vainonen J, Gauthier A, Overmyer K, Kangasjärvi J. Reactive oxygen in abiotic stress perception-from genes to proteins: InTech. 2011.
Trevaskis B, Watts RA, Andersson CR, Llewellyn DJ, Hargrove MS, Olson JS, Dennis ES, Peacock WJ. Two hemoglobin genes in Arabidopsis thaliana: the evolutionary origins of leghemoglobins. Proc Natl Acad Sci. 1997;94(22):12230–4.
Hunt P, Klok E, Trevaskis B, Watts R, Ellis M, Peacock W, Dennis E. Increased level of hemoglobin 1 enhances survival of hypoxic stress and promotes early growth in Arabidopsis thaliana. Proc Natl Acad Sci. 2002;99(26):17197–202.
Perazzolli M, Dominici P, Romero-Puertas MC, Zago E. Zeier Jr, Sonoda M, Lamb C, Delledonne M: Arabidopsis nonsymbiotic hemoglobin AHb1 modulates nitric oxide bioactivity. Plant Cell. 2004;16(10):2785–94.
Zheng W, Ma L, Zhao J, Li Z, Sun F, Lu X. Comparative transcriptome analysis of two rice varieties in response to rice stripe virus and small brown planthoppers during early interaction. PLoS One. 2013;8(12):e82126.
Liu CW, Hsu YK, Cheng YH, Yen HC, Wu YP, Wang CS, Lai CC. Proteomic analysis of salt-responsive ubiquitin-related proteins in rice roots. Rapid Commun Mass Spectrom. 2012;26(15):1649–60.
Endler A, Kesten C, Schneider R, Zhang Y, Ivakov A, Froehlich A, Funke N, Persson S. A mechanism for sustained cellulose synthesis during salt stress. Cell. 2015;162(6):1353–64.
Román Á, Andreu V, Hernández ML, Lagunas B, Picorel R, Martínez-Rivas JM, Alfonso M. Contribution of the different omega-3 fatty acid desaturase genes to the cold response in soybean. J Exp Bot. 2012;63(13):4973–82.
Iba K. Acclimative response to temperature stress in higher plants: approaches of gene engineering for temperature tolerance. Ann Rev Plant Biol. 2002;53(1):225–45.
Murata N, Ishizaki-Nishizawa O, Higashi S, Hayashi H, Tasaka Y, Nishida I. Genetically engineered alteration in the chilling sensitivity of plants. Nature. 1992;356(6371):710–3.
Cai X, Lytton J. The cation/Ca2+ exchanger superfamily: phylogenetic analysis and structural implications. Mol Biol Evol. 2004;21(9):1692–703.
Wang J, Wang J, Wang X, Li R, Chen B. Proteomic response of hybrid wild rice to cold stress at the seedling stage. PloS One. 2018;13(6):e0198675.
Thomashow MF. Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Ann Rev Plant Biol. 1999;50(1):571–99.
Ruelland E, Vaultier M-N, Zachowski A, Hurry V. Cold signalling and cold acclimation in plants. Adv Bot Res. 2009;49:35–150.
Hennessey TL, Field CB. Circadian rhythms in photosynthesis: oscillations in carbon assimilation and stomatal conductance under constant conditions. Plant Physiol. 1991;96(3):831–6.
Hennessey TL, Freeden AL, Field CB. Environmental effects on circadian rhythms in photosynthesis and stomatal opening. Planta. 1993;189(3):369–76.
Mansfield C, Carabasi R, Wells W, Borman K. Circadian rhythm in the skin temperature of normal and cancerous breasts. Int J Chronobiol. 1973;1(3):235.
Allen DJ, Ort DR. Impacts of chilling temperatures on photosynthesis in warm-climate plants. Trends Plant Sci. 2001;6(1):36–42.
Jalil SU, Ansari MI. Isoprenoids in Plant Protection Against Abiotic Stress. In: Protective Chemical Agents in the Amelioration of Plant Abiotic Stress: Biochemical and Molecular Perspectives. 2020. p. 424–36.
Brunetti C, Guidi L, Sebastiani F, Tattini M. Isoprenoids and phenylpropanoids are key components of the antioxidant defense system of plants facing severe excess light stress. Environ Exp Bot. 2015;119:54–62.
Kan C-C, Chung T-Y, Wu H-Y, Juo Y-A, Hsieh M-H. Exogenous glutamate rapidly induces the expression of genes involved in metabolism and defense responses in rice roots. BMC Genom. 2017;18(1):1–17.
Decottignies A, Goffeau A. Complete inventory of the yeast ABC proteins. Nat Genet. 1997;15(2):137.
Sánchez-Fernández RO, Davies TE, Coleman JO, Rea PA. The Arabidopsis thaliana ABC protein superfamily, a complete inventory. J Biol Chem. 2001;276(32):30231–44.
Martinoia E, Klein M, Geisler M, Bovet L, Forestier C, Kolukisaoglu Ü, Müller-Röber B, Schulz B. Multifunctionality of plant ABC transporters–more than just detoxifiers. Planta. 2002;214(3):345–55.
Moons A. Ospdr9, which encodes a PDR-type ABC transporter, is induced by heavy metals, hypoxic stress and redox perturbations in rice roots 1. Febs Letters. 2003;553(3):370–6.
Yazaki K, Shitan N, Sugiyama A, Takanashi K. Cell and molecular biology of ATP-binding cassette proteins in plants. Int Rev Cell Mol Biol. 2009;276:263–99.
Fukao T, Xu K, Ronald PC, Bailey-Serres J. A variable cluster of ethylene response factor–like genes regulates metabolic and developmental acclimation responses to submergence in rice. Plant Cell. 2006;18(8):2021–34.
Hattori Y, Nagai K, Furukawa S, Song X-J, Kawano R, Sakakibara H, Wu J, Matsumoto T, Yoshimura A, Kitano H. The ethylene response factors SNORKEL1 and SNORKEL2 allow rice to adapt to deep water. Nature. 2009;460(7258):1026.
Hinz M, Wilson IW, Yang J, Buerstenbinder K, Llewellyn D, Dennis ES, Sauter M, Dolferus R. Arabidopsis RAP2. 2: an ethylene response transcription factor that is important for hypoxia survival. Plant Physiol. 2010;153(2):757–72.
Butt HI, Yang Z, Gong Q, Chen E, Wang X, Zhao G, Ge X, Zhang X, Li F. GaMYB85, an R2R3 MYB gene, in transgenic Arabidopsis plays an important role in drought tolerance. BMC Plant Biol. 2017;17(1):142.
Tak H, Negi S, Ganapathi T. Overexpression of MusaMYB31, a R2R3 type MYB transcription factor gene indicate its role as a negative regulator of lignin biosynthesis in banana. PloS one. 2017;12(2):e0172695.
Xu Z-S, Feng K, Que F, Wang F, Xiong A-S. A MYB transcription factor, DcMYB6, is involved in regulating anthocyanin biosynthesis in purple carrot taproots. Sci Rep. 2017;7:45324.
Huang W, Lv H, Wang Y. Functional characterization of a novel R2R3-MYB transcription factor modulating the flavonoid biosynthetic pathway from Epimedium sagittatum. Front Plant Sci. 2017;8:1274.
Lv Y, Yang M, Hu D, Yang Z, Ma S, Li X, Xiong L. The OsMYB30 transcription factor suppresses cold tolerance by interacting with a JAZ protein and suppressing β-amylase expression. Plant Physiol. 2017;173(2):1475–91.
Wang F-Z, Chen M-X, Yu L-J, Xie L-J, Yuan L-B, Qi H, Xiao M, Guo W, Chen Z, Yi K. OsARM1, an R2R3 MYB transcription factor, is involved in regulation of the response to arsenic stress in rice. Front Plant Sci. 1868;2017:8.
Wei Q, Luo Q, Wang R, Zhang F, He Y, Zhang Y, Qiu D, Li K, Chang J, Yang G. A wheat R2R3-type MYB transcription factor TaODORANT1 positively regulates drought and salt stress responses in transgenic tobacco plants. Front Plant Sci. 2017;8:1374.
Zhou Y, Yang P, Cui F, Zhang F, Luo X, Xie J. Transcriptome analysis of salt stress responsiveness in the seedlings of Dongxiang wild rice (Oryza rufipogon Griff.). PLoS One. 2016;11(1):e0146242.
Garg R, Verma M, Agrawal S, Shankar R, Majee M, Jain M. Deep transcriptome sequencing of wild halophyte rice, Porteresia coarctata, provides novel insights into the salinity and submergence tolerance factors. DNA Res. 2013;21(1):69–84.
Agarwal M, Hao Y, Kapoor A, Dong C-H, Fujii H, Zheng X, Zhu J-K. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing tolerance. J Biol Chem. 2006;281(49):37636–45.
Guan S, Xu Q, Ma D, Zhang W, Xu Z, Zhao M, Guo Z. Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice. Gene. 2019;685:96–105.
Suh J, Jeung J, Lee J, Choi Y, Yea J, Virk P, Mackill D, Jena K. Identification and analysis of QTLs controlling cold tolerance at the reproductive stage and validation of effective QTLs in cold-tolerant genotypes of rice (Oryza sativa L.). Theor Appl Genet. 2010;120(5):985–95.
Sharoni AM, Nuruzzaman M, Satoh K, Shimizu T, Kondoh H, Sasaya T, Choi I-R, Omura T, Kikuchi S. Gene structures, classification and expression models of the AP2/EREBP transcription factor family in rice. Plant Cell Physiol. 2010;52(2):344–60.
Ramamoorthy R, Jiang S-Y, Kumar N, Venkatesh PN, Ramachandran S. A comprehensive transcriptional profiling of the WRKY gene family in rice under various abiotic and phytohormone treatments. Plant Cell Physiol. 2008;49(6):865–79.
Berri S, Abbruscato P, Faivre-Rampant O, Brasileiro AC, Fumasoni I, Satoh K, Kikuchi S, Mizzi L, Morandini P, Pè ME. Characterization of WRKY co-regulatory networks in rice and Arabidopsis. BMC Plant Biol. 2009;9(1):120.
Kotak S, Larkindale J, Lee U, von Koskull-Döring P, Vierling E, Scharf K-D. Complexity of the heat stress response in plants. Curr Opin Plant Biol. 2007;10(3):310–6.
Ciftci-Yilmaz S, Mittler R. The zinc finger network of plants. Cell Mol Life Sci. 2008;65(7–8):1150–60.
Agarwal P, Arora R, Ray S, Singh AK, Singh VP, Takatsuji H, Kapoor S, Tyagi AK. Genome-wide identification of C 2 H 2 zinc-finger gene family in rice and their phylogeny and expression analysis. Plant Mol Biol. 2007;65(4):467–85.
Miller G, Shulaev V, Mittler R. Reactive oxygen signaling and abiotic stress. Physiol Plant. 2008;133(3):481–9.
Kiełbowicz-Matuk A. Involvement of plant C2H2-type zinc finger transcription factors in stress responses. Plant Sci. 2012;185:78–85.
Huang J, Sun S, Xu D, Lan H, Sun H, Wang Z, Bao Y, Wang J, Tang H, Zhang H. A TFIIIA-type zinc finger protein confers multiple abiotic stress tolerances in transgenic rice (Oryza sativa L.). Plant Mol Biol. 2012;80(3):337–50.
Huang J, Sun S-J, Xu D-Q, Yang X, Bao Y-M, Wang Z-F, Tang H-J, Zhang H. Increased tolerance of rice to cold, drought and oxidative stresses mediated by the overexpression of a gene that encodes the zinc finger protein ZFP245. Biochem Biophys Res Commun. 2009;389(3):556–61.
Huang J, Yang X, Wang M-M, Tang H-J, Ding L-Y, Shen Y, Zhang H-S. A novel rice C2H2-type zinc finger protein lacking DLN-box/EAR-motif plays a role in salt tolerance. Biochim Biophys Acta (BBA) Gene Struct Express. 2007;1769(4):220–7.
Xu D-Q, Huang J, Guo S-Q, Yang X, Bao Y-M, Tang H-J, Zhang H-S. Overexpression of a TFIIIA-type zinc finger protein gene ZFP252 enhances drought and salt tolerance in rice (Oryza sativa L.). FEBS Lett. 2008;582(7):1037–43.
Sun S-J, Guo S-Q, Yang X, Bao Y-M, Tang H-J, Sun H, Huang J, Zhang H-S. Functional analysis of a novel Cys2/His2-type zinc finger protein involved in salt tolerance in rice. J Exp Bot. 2010;61(10):2807–18.
Kanneganti V, Gupta AK. Overexpression of OsiSAP8, a member of stress associated protein (SAP) gene family of rice confers tolerance to salt, drought and cold stress in transgenic tobacco and rice. Plant Mol Biol. 2008;66(5):445–62.
Kong Z, Li M, Yang W, Xu W, Xue Y. A novel nuclear-localized CCCH-type zinc finger protein, OsDOS, is involved in delaying leaf senescence in rice. Plant Physiol. 2006;141(4):1376–88.
Zheng X, Chen B, Lu G, Han B. Overexpression of a NAC transcription factor enhances rice drought and salt tolerance. Biochem Biophys Res Commun. 2009;379(4):985–9.
Lin R, Zhao W, Meng X, Wang M, Peng Y. Rice gene OsNAC19 encodes a novel NAC-domain transcription factor and responds to infection by Magnaporthe grisea. Plant Sci. 2007;172(1):120–30.
Nakashima K, Tran LSP, Van Nguyen D, Fujita M, Maruyama K, Todaka D, Ito Y, Hayashi N, Shinozaki K, Yamaguchi-Shinozaki K. Functional analysis of a NAC-type transcription factor OsNAC6 involved in abiotic and biotic stress-responsive gene expression in rice. Plant J. 2007;51(4):617–30.
Zhang T, Zhao X, Wang W, Pan Y, Huang L, Liu X, Zong Y, Zhu L, Yang D, Fu B. Comparative transcriptome profiling of chilling stress responsiveness in two contrasting rice genotypes. PloS one. 2012;7(8):e43274.
Chinnusamy V, Ohta M, Kanrar S. Lee B-h, Hong X, Agarwal M, Zhu J-K: ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev. 2003;17(8):1043–54.
Bourne Y, Henrissat B. Glycoside hydrolases and glycosyltransferases: families and functional modules. Curr Opin Struct Biol. 2001;11(5):593–600.
Henrissat B. A classification of glycosyl hydrolases based on amino acid sequence similarities. Biochem J. 1991;280(2):309–16.
Chandran AKN, Jeong HY, Jung K-H, Lee C. Development of functional modules based on co-expression patterns for cell-wall biosynthesis related genes in rice. J Plant Biol. 2016;59(1):1–15.
Dametto A, Sperotto RA, Adamski JM, Blasi ÉA, Cargnelutti D, de Oliveira LF, Ricachenevsky FK, Fregonezi JN, Mariath JE, da Cruz RP. Cold tolerance in rice germinating seeds revealed by deep RNAseq analysis of contrasting indica genotypes. Plant Sci. 2015;238:1–12.
Kuntothom T, Luang S, Harvey AJ, Fincher GB, Opassiri R, Hrmova M, Cairns JRK. Rice family GH1 glycoside hydrolases with β-D-glucosidase and β-D-mannosidase activities. Arch Biochem Biophys. 2009;491(1–2):85–95.
Wang E, Xu X, Zhang L, Zhang H, Lin L, Wang Q, Li Q, Ge S, Lu B-R, Wang W. Duplication and independent selection of cell-wall invertase genes GIF1 and OsCIN1 during rice evolution and domestication. BMC Evol Biol. 2010;10(1):1–13.
Guo Z, Liu C, Xiao W, Wang R, Zhang L, Guan S, Zhang S, Cai L, Liu H, Huang X. Comparative transcriptome profile analysis of anther development in reproductive stage of rice in cold region under cold stress. Plant Mol Biology Rep. 2019;37(3):129–45.
Sperotto RA, de Araújo Junior AT, Adamski JM, Cargnelutti D, Ricachenevsky FK, de Oliveira B-HN, da Cruz RP, Dos Santos RP, da Silva LP, Fett JP. Deep RNAseq indicates protective mechanisms of cold-tolerant indica rice plants during early vegetative stage. Plant Cell Rep. 2018;37(2):347–75.
Kong F-J, Oyanagi A, Komatsu S. Cell wall proteome of wheat roots under flooding stress using gel-based and LC MS/MS-based proteomics approaches. Biochim Biophys Acta (BBA) Proteins Proteomics. 2010;1804(1):124–36.
Calzadilla PI, Maiale SJ, Ruiz OA, Escaray FJ. Transcriptome response mediated by cold stress in Lotus japonicus. Front Plant Sci. 2016;7:374.
Usadel B, Blaesing OE, Gibon Y, Poree F, Hoehne M, Guenter M, Trethewey R, Kamlage B, Poorter H, Stitt M. Multilevel genomic analysis of the response of transcripts, enzyme activities and metabolites in Arabidopsis rosettes to a progressive decrease of temperature in the non-freezing range. Plant Cell Environ. 2008;31(4):518–47.
Wang X, Li W, Li M, Welti R. Profiling lipid changes in plant response to low temperatures. Physiol Plant. 2006;126(1):90–6.
Li K-H, Yu Y-H, Dong H-J, Zhang W-B, Ma J-C, Wang H-H. Biological functions of ilvC in branched-chain fatty acid synthesis and diffusible signal factor family production in Xanthomonas campestris. Front Microbiol. 2017;8:2486.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51(D1):D587–92.
Bray EA. Genes commonly regulated by water-deficit stress in Arabidopsis thaliana. J Exp Bot. 2004;55(407):2331–41.
Hu R, Zhu X, Xiang S, Zhan Y, Zhu M, Yin H, Zhou Q, Zhu L, Zhang X, Liu Z. Comparative transcriptome analysis revealed the genotype specific cold response mechanism in tobacco. Biochem Biophys Res Commun. 2016;469(3):535–41.
Muthurajan R, Rahman H, Manoharan M, Ramanathan V, Nallathambi J. Drought responsive transcriptome profiling in roots of contrasting rice genotypes. Indian J Plant Physiol. 2018;23(3):393–407.
Suwabe K, Suzuki G, Takahashi H, Shiono K, Endo M, Yano K, Fujita M, Masuko H, Saito H, Fujioka T. Separated transcriptomes of male gametophyte and tapetum in rice: validity of a laser microdissection (LM) microarray. Plant Cell Physiol. 2008;49(10):1407–16.
Magwanga RO, Lu P, Kirungu JN, Lu H, Wang X, Cai X, Zhou Z, Zhang Z, Salih H, Wang K. Characterization of the late embryogenesis abundant (LEA) proteins family and their role in drought stress tolerance in upland cotton. BMC Genet. 2018;19(1):1–31.
İbrahime M, Kibar U, Kazan K, Özmen CY, Mutaf F, Aşçı SD, Aydemir BÇ, Ergül A. Genome-wide identification of the LEA protein gene family in grapevine (Vitis vinifera L.). Tree Genet Genom. 2019;15(4):1–14.
Gao J, Lan T. Functional characterization of the late embryogenesis abundant (LEA) protein gene family from Pinus tabuliformis (Pinaceae) in Escherichia coli. Sci Rep. 2016;6(1):1–10.
Reddy PS, Reddy GM, Pandey P, Chandrasekhar K, Reddy MK. Cloning and molecular characterization of a gene encoding late embryogenesis abundant protein from Pennisetum glaucum: protection against abiotic stresses. Mol Biol Rep. 2012;39(6):7163–74.
Mizoi J, Yamaguchi-Shinozaki K. Molecular approaches to improve rice abiotic stress tolerance. In: Rice protocols. 2013. p. 269–83.
Wang X-S, Zhu H-B, Jin G-L, Liu H-L, Wu W-R, Zhu J. Genome-scale identification and analysis of LEA genes in rice (Oryza sativa L.). Plant Science. 2007;172(2):414–20.
Liu WX, Liu HL, Qu LQ. Embryo-specific expression of soybean oleosin altered oil body morphogenesis and increased lipid content in transgenic rice seeds. Theor Appl Genet. 2013;126(9):2289–97.
Ganguly M, Datta K, Roychoudhury A, Gayen D, Sengupta DN, Datta SK. Overexpression of Rab16A gene in indica rice variety for generating enhanced salt tolerance. Plant Signal Behav. 2012;7(4):502–9.
Fu C, Wang F, Liu W, Liu D, Li J, Zhu M, Liao Y, Liu Z, Huang H, Zeng X. Transcriptomic analysis reveals new insights into high-temperature-dependent glume-unclosing in an elite rice male sterile line. Front Plant Sci. 2017;8:112.
Bioinformatics B. FastQC: a quality control tool for high throughput sequence data. Cambridge: Babraham Institute; 2011.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.
Sasaki T, Burr B. International Rice Genome Sequencing Project: the effort to completely sequence the rice genome. Curr Opin Plant Biol. 2000;3(2):138–42.
Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013;6(1):1–10.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
Du Z, Zhou X, Ling Y, Zhang Z, Su Z. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 2010;38(suppl_2):W64–70.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Stat Soc B (Methodological). 1995;57(1):289–300.
Ge SX, Jung D, Yao R. ShinyGO: a graphical gene-set enrichment tool for animals and plants. Bioinformatics. 2020;36(8):2628–9.
Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, Müller LA, Rhee SY, Stitt M. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004;37(6):914–39.
Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P. STRING v11: protein–protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607–13.
Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
We thank the Genomic and Bioinformatics Service at Texas A&M Agrilife for the sequencing data.
This project was supported by Texas A&M AgriLife Research and USDA NIFA # 2022-67013-36210 to Endang M. Septiningsih.
Ethics approval and consent to participate
The rice germplasm used in this study were obtained from the United States Department of Agriculture (USDA) National Plant Germplasm System (NPGS) and from the inbred rice breeding program of the Texas A&M AgriLife Research at Beaumont. The collection and use of rice germplasm, as well as the methods carried out in this study comply with relevant institutional, national, and international guidelines and legislation.
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
PCA plot generated using regularized log transformation of total gene counts with DESeq2 under (a) flooding during germination, (b) cold stress, and (c) combined of flooding and cold stress during germination conditions. Supplementary Figure 2. Differentially expressed genes (DEGs) under flooding during germination conditions. (a) DEGs of ABC transporter family; (b) DEGs of AP2 domain containing protein family; (c) DEGs of helix-loop-helix DNA-binding domain containing protein family; (d) DEGs of pentatricopeptide (PPR) gene family; (e) DEGs of no apical meristem (NAM) protein family; (f) DEGs of RNA binding family; (g) DEGs of MYB transcription factor family; (h) DEGs of WRKY family; and (i) DEGs of C2H2 zinc finger protein family. Supplementary Figure 3. DEGs of Darij under a) flooding during germination, (b) cold stress, and (c) combined of flooding and cold stress during germination conditions were binned to MapMan metabolism bin. Up-regulated and down-regulated transcripts are shown in blue and red, respectively. Supplementary Figure 4. DEGs of 4610 under a) flooding during germination, (b) cold stress, and (c) combined of flooding and cold stress during germination conditions were binned to MapMan metabolism bin. Up-regulated and down-regulated transcripts are shown in blue and red, respectively. Supplementary Figure 5. DEGs of Darij under a) flooding during germination, (b) cold stress, and (c) combined of flooding and cold stress during germination conditions associated with secondary metabolism were binned to MapMan functional categories. Up-regulated and down-regulated transcripts are shown in blue and red, respectively. Supplementary Figure 6. DEGs of 4610 under a) flooding during germination, (b) cold stress, and (c) combined of flooding and cold stress during germination conditions associated with secondary metabolism were binned to MapMan functional categories. Up-regulated and down-regulated transcripts are shown in blue and red, respectively.
Summary of RNA seq data and sequence assembly of control samples for flooding during germination stress. Supplementary Table 2. Summary of RNA seq data and sequence assembly of samples under flooding during germination stress. Supplementary Table 3. Summary of RNA seq data and sequence assembly of control samples for cold stress. Supplementary Table 4. Summary of RNA seq data and sequence assembly of samples under cold stress. Supplementary Table 5. Summary of RNA seq data and sequence assembly of samples under combined flooding and cold stress during germination conditions.
Differentially expressed genes (DEGs) identified independently for each genotype (Darij (AG1) and line 4610 (AG2)) and each stress condition (germinating seeds under flooding (anaerobic germination, AG).
Differentially expressed genes (DEGs) identified by considering both genotype and treatment effects.
GO enrichment analysis in each stress condition.
About this article
Cite this article
Thapa, R., Tabien, R.E., Johnson, C.D. et al. Comparative transcriptomic analysis of germinating rice seedlings to individual and combined anaerobic and cold stress. BMC Genomics 24, 185 (2023). https://doi.org/10.1186/s12864-023-09262-z
- Rice (Oryza sativa L.)
- Anaerobic germination (AG) stress
- Cold stress during germination
- Combined stress
- Differentially expressed genes (DEGs)
- Hub genes