Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis
BMC Plant Biology volume 19, Article number: 352 (2019)
Rice plants show yellowing, stunting, withering, reduced tillering and utimately low productivity in susceptible varieties under low temperature stress. Comparative transcriptome analysis was performed to identify novel transcripts, gain new insights into different gene expression and pathways involved in cold tolerance in rice.
Comparative transcriptome analyses of 5 treatments based on chilling stress exposure revealed more down regulated genes in susceptible and higher up regulated genes in tolerant genotypes. A total of 13930 and 10599 differentially expressed genes (DEGs) were detected in cold susceptible variety (CSV) and cold tolerant variety (CTV), respectively. A continuous increase in DEGs at 6, 12, 24 and 48 h exposure of cold stress was detected in both the genotypes. Gene ontology (GO) analysis revealed 18 CSV and 28 CTV term significantly involved in molecular function, cellular component and biological process. GO classification showed a significant role of transcription regulation, oxygen, lipid binding, catalytic and hydrolase activity for tolerance response. Absence of photosynthesis related genes, storage products like starch and synthesis of other classes of molecules like fatty acids and terpenes during the stress were noticed in susceptible genotype. However, biological regulations, generation of precursor metabolites, signal transduction, photosynthesis, regulation of cellular process, energy and carbohydrate metabolism were seen in tolerant genotype during the stress. KEGG pathway annotation revealed more number of genes regulating different pathways resulting in more tolerant. During early response phase, 24 and 11 DEGs were enriched in CTV and CSV, respectively in energy metabolism pathways. Among the 1583 DEG transcription factors (TF) genes, 69 WRKY, 46 bZIP, 41 NAC, 40 ERF, 31/14 MYB/MYB-related, 22 bHLH, 17 Nin-like 7 HSF and 4C3H were involved during early response phase. Late response phase showed 30 bHLH, 65 NAC, 30 ERF, 26/20 MYB/MYB-related, 11 C3H, 12 HSF, 86 Nin-like, 41 AP2/ERF, 55 bZIP and 98 WRKY members TF genes. The recovery phase included 18 bHLH, 50 NAC, 31 ERF, 24/13 MYB/MYB-related, 4 C3H, 4 HSF, 14 Nin-like, 31 bZIP and 114 WRKY TF genes.
Transcriptome analysis of contrasting genotypes for cold tolerance detected the genes, pathways and transcription factors involved in the stress tolerance.
Dependence of more than 50% of world population on rice for carbohydrate source makes it one of the most important food crops. Rice is cultivated in a wide range of ecology starting from upland to deepwater ecology including varying altitudes, climate and diverse soil types. Abiotic stresses like cold, high temperature, drought, salinity and flooding affect adversely in rice production. Rice being chilling-sensitive plant shows yellowing, slow seedling growth, stunting, withering, reduced tillering and utimately low productivity in cold susceptible varieties under low temperature stress [1,2,3,4]. Globally, about 15 million hectares of rice fields in 24 countries are affected by cold weather. Cultivation of rice is not possible in approximately 7 million hectares of land in south and south east Asia due to cold stress . Now-a-days, the extent of unexpected cold weather has increased due to extreme weather events. The yield of rice in Australia, Korea, India, China, Japan and many other countries is affected by cold stress [2, 3, 6]. The dry season rice of India particularly boro rice faces cold stress during seedling stage. The stress covers around 4 million hectares of rice areas of the country targeting the seedling-stage causing delay in growth of the plant which ultimately coincide with high temperature stress during flowering stage [3, 7]. Therefore, developing high yielding varieties possessing cold stress tolerance may help in further augmenting rice production in these vulnerable regions. Besides, more insight into the molecular mechanisms of tolerance in response to chilling stress need to be revealed.
Plants adopt different strategies to escape from stresses based on their occurrence, severity and ecology . Stresses faced by plants in response to wounding, pathogen attack, UV exposure, flood, drought, salt, heat and cold stress and combination of these stresses enhances production of reactive oxygen species (ROS) levels [9,10,11,12]. ROS homeostasis is observed to regulate strongly the cold responsive genes in arabidopsis. Earlier research results revealed the role of cold stress binding factors (CBFs) in ROS detoxification [13, 14]. Besides, ROS-mediated regulatory module that functions as an early component of the chilling stress response pathway is reported in rice . Reactive oxygen species may be taken as a point of convergence for various gene networks in response to stress situation . Transcription factors (TFs) are the master regulators for controlling expression of many target genes in a net work mode of single TF of down-stream genes. Earlier results suggest the role of AREB/ABF, NAC regulons and DREB1/ CBF in regulating salinity, temperature and drought stress responsive gene expression in rice [16, 17]. Reports of heat shock factor (HSFs) genes have been suggested as an important node of cross-talk in rice . Research findings of a common set of TFs working in convergence with stress signaling networks and regulons of the upstream regulatory genes for controlling plant responses to multiple stresses . A network mode of functioning of target TFs for similar molecular, biochemical, genetical and physiological processes are performed in response to stresses. A structural basis for the generation of unique patterns of gene expression is provided by binding of the TFs in combinations with cis-element. In a recent cloning study of COLD1 quantitative trait locus revealed mechanism of chilling tolerance in which the QTL interact with alpha sub-unit of G-protein for activating the Ca2+ channel which enhance the GTPase activity of G-protein for expressing chilling tolerance in japonica rice . However, the mechanisms underlying rice chilling tolerance is still not fully known.
Many researchers previously studied rice transcriptome using microarrays [20,21,22,23,24,25]. In recent years, researchers used RNA-seq technology to profile and compare the transcriptome and understand molecular pathways in rice [26,27,28,29,30,31]. The RNA-seq technology has true advantage of representing the global transcriptome, which is limited in microarray since it depends on predefined probes. Rice genotype, Geetanjali was categorized as very highly tolerant and Sahabhagidhan as highly susceptible type to chilling stress reported by earlier publications [2,3,4]. The highly tolerant genotype Geetanjali could survive upto 25 days under chilling stress of 4 °C. Therefore, in the present investigation, we have compared the transcriptome of contrasting genotypes for chilling stress tolerance to identify novel transcripts, gain new insights into different gene expression and pathways involved in cold tolerance in rice.
Phenotyping of contrasting genotypes for chilling stress tolerance
Upon exposure to low temperature stress, the susceptible genotype, Shahabhagidhan turned pale green and rolled leaves exhibiting a score of 5 after 7 days of exposure at 15 °C; whereas the resistant genotype, Geetanjali remained normal in growth and leaf coloration. At 4 °C, Sahabhagidhan had rolled leaf with yellow to brown in color and reduced growth showing a score of 7 while tolerant genotype, Geetanjali was normal. After 7 days of chilling at 4 °C, susceptible genotype, Sahabhagidhan died showing a SES score of 9 while Geetanjali was highly tolerant with a score of 1. The tolerant genotype showed scores of 5, 7 and 9 at chilling stress of 4 °C on 15th, 20th and 25th day, respectively.
RNA sequencing and differentially expressed genes in response to chilling stress
RNA-seq generated 566 million raw reads in which 431 million high quality reads (QV > 25) with an average of 36 million reads per sample were obtained (Additional file 7: Table S1). High quality reads were mapped to the rice reference genome of cv. Nipponbare (Oryza sativa L. subsp. japonica) ranging from 92.4–95.1% per sample.
Differentially expressed genes for chilling stress and subsequent recovery conditions were identified using fold change (FC) greater than 0 and significant in ‘t’ test (p < 0.05). We detected a total of 24529 DEGs in both the genotypes of which 13930 and 10599 DEGs were observed in CSV and CTV, respectively (Fig. 1). Out of these, 1666 and 2206 DEGs were found in subsequent recovery condition of CSV and CTV, respectively (Fig. 1). Furthermore, there was significant increase in DEGs detected in both the genotypes at 6, 12, 24 and 48 h exposure of cold stress with FC value > 4. To assess similarities and differences between cold contrasting genotypes, hierarchical cluster analysis using CuumeRbund was performed. The down-regulated gene populations in CSV was more than up-regulated genes while reverse was observed in CTV (Table 2). Further 12, 24 and 48 h cold stressed samples of both genotypes cluster together, but at 6 h cold stressed samples of both genotypes formed separate cluster. In addition, 24 h recovery samples of these two genotypes were grouped into a single cluster (Additional file 1: Figure S1). The comparative profile of differentially expressed genes in both the genotypes showed that more numbers of DEGs were down-regulated as compared to up-regulated genes except 48 h recovery condition (Additional file 8: Table S2; Additional file 2: Figure S2). Genes up-regulation and down-regulation in CSV above FC value was 24.97 and 75.03%, respectively. But in CTV, 43.92 and 56.08% genes were detected to be up-regulated and down-regulated, respectively (Additional file 2: Figure S2; Additional file 9: Table S3; Additional file 10: Table S4; Additional file 11: Table S5).
Functional classification of DEGs
Out of total 24529 DEGs, 14632 DEGs were GO term enriched in which 8047 DEGs observed in CSV and 6585 DEGs in CTV of which at least one term of GO classification categorized as molecular function, cellular component and biological process. Among all, 18 and 28 GO terms were significantly involved in molecular function, cellular component and biological process of CSV and CTV, respectively (Figs. 2 and 3).
Under molecular function category, five GO terms (GO:0003700; GO: GO:0003677; GO:0003824; GO:0016787; GO:0008289) were enriched and played a significant role in transcription regulation, oxygen, DNA and lipid binding, catalytic and hydrolase activity at any of three phases of CSV while nine GO terms (GO:0019825; GO:0030246; GO:0016740; GO:0003700; GO:0016301; GO:0016772; GO:0003824; GO:0005215; GO:0016787) were enriched in CTV. Both the genotypes showed 3 common GO terms enrichment during the cold stress. Catalytic activity, hydrolase activity and lipid binding molecular function were not enriched during the stress period in susceptible genotype. The cellular component category was enriched with six GO terms (GO: 0030312: GO: 0019825; 0005618; GO: 0005576; GO: 0009579; GO: 0009536) of which 3 GO terms (GO: 0030312; GO: 0005618; GO: 0005576) were significantly involved in all three phase of CSV genotype while in CTV, these were enriched in 24 h recovery condition. The GO term GO: 0009536 was enriched in 24 h recovery in CSV. The biological process category was enriched with 7 and 15 GO terms in at least any of the phases of CSV and CTV, respectively. Amongst these, 6 GO terms namely secondary metabolic process, response to stimulus, response to endogenous stimulus, response to stress, response to abiotic stimulus and response to biotic stimulus were enriched in both the genotypes. Examination of lipid metabolic process showed the same biological function in both susceptible and tolerant genotypes while signal transduction, regulation of cellular process, pollen-pistil interaction, regulation of biological process, photosynthesis, generation of precursor metabolites and energy, carbohydrate metabolic process were detected and enriched in early response phase of CTV (GO: 0007165; GO: 0050794; GO: 0065007; GO: 0009875; GO: 0050789; GO: 0015979; GO: 0006091; GO: 0005975). Further, DEGs in all the three response phases were described on the basis of GO term enrichment. In early response phase, the DEGs in biological process, molecular function and cellular component were detected in response to cold stress (Additional file 11: Table S5; Additional file 17: Table S11). It was highly surprising that no DEGs of CTV enriched for cellular component, whereas 102 DEGs of CSV were represented mainly for cell wall and all of them were down regulated. Similarly, the DEGs obtained for 3 classes of GO categorization in late response phase enriched due to cold stress (Additional file 12: Table S6; Additional file 15: Table S9; Additional file 18: Table S12). In recovery phase, CSV showed DEGs in 6 biological categories (secondary metabolic process, response to stimulus, response to abiotic stimulus, lipid metabolic process, response to endogenous stimulus, and response to stress) and no DEGs represented for transcription factor activity/regulation. In CTV, biological process DEGs were represented in 8 functional categories, out of which response to stimulus and stress, metabolic process, lipid and carbohydrate metabolism were affected (Additional file 13: Table S7; Additional file 16: Table S10).
Kyoto encyclopedia of genes and genomes (KEGG) pathway annotation
KEGG pathway database was used to identify cold stress related genes involved in the pathways using DEGs. It was found that cold treatment continuously increased the number of up-regulated and down-regulated genes in both CSV and CTV genotypes, while it decreased after 24 h recovery, except exclusively expressed genes of CTV during recovery (Table 2; Additional file 2: Figure S2). Comparative analysis of two genotypes under control condition revealed 29216 and 30867 genes expressed in susceptible and tolerant genotypes, respectively.
The commonly expressed genes (26009) in both genotypes may be involved in plant growth and development. In early response phase, out of 2173 DEGs, 297 DEGs were enriched in KEGG pathway in CSV while 243 DEGs were enriched in CTV from 1655 DEGs (Fig. 3: Additional file 5: Figure S5). Comparison of early response phase in both the genotypes indicated that around 46% of DEGs were up-regulated and 57% DEGs were down regulated in CTV while 20% DEGs were up-regulated and 80% DEGs were down-regulated in CSV. Furthermore, comparison of CSV and CTV, DEGs involved in KEGG pathway indicates that around 49 DEGs were commonly involved in early response phase while 155 and 142 DEGs exclusively present in the contrasting two varieties revealed cold susceptible and cold tolerant DEGs, respectively. Furthermore, comparison at late response phase (T4) and recovery response phase (T5) revealed that 324 genes were up-regulated exclusively in T4 and 195 were down-regulated in T5 in which 4 genes (Os02g44230/Os02g44235, trehalose 6-phosphate phosphatase; Os07g42960, 3-deoxy-7-phosphoheptulonate synthase; Os08g06344, ATP-dependent RNA helicase DDX46/PRP5) were significantly down-regulated and up-regulated in T5 and T4, respectively (Fig. 3; Additional file 19: Table S13 and Additional file 20: Table S14). Interestingly, out of four genes, three genes (Os02g44230, Os02g44235, Os08g06344) were up-regulated in S4 (48h late response phase of susceptible genotype) and one gene (Os08g06344) was down-regulated in S5 (24 h recovery condition of susceptible genotype). More interestingly, out of 324 up-regulated genes of T4, 88 genes were also up-regulated in S4 genotype. In the susceptible genotype Sahabhagidhan, our results showed that 702 genes were significantly down-regulated in late response phase in which 11 genes (Os03g09020, S- hydroxymethyl glutathione dehydrogenase/alcohol dehydrogenase; Os04g38540, aldose 1-epimerase; Os12g03470, alpha-N-arabinofuranosidase; Os07g34520, isocitrate lyase; Os07g05365, photosystem II 10 kDa protein; Os12g08260, 2-oxoisovalerate dehydrogenase E1 component; Os05g43830, proline iminopeptidase; Os10g02070, peroxidase; Os03g15960, HSP20 family protein; Os04g53210, (S)-2-hydroxy-acid oxidase; Os03g22680, RING finger and CHY zinc finger domain-containing protein) were significantly up-regulated in tolerant genotype. Out of 11, 9 genes (except Os10g02070, Os03g15960) were commonly expressed in tolerant genotype at 24 h recovery phase. Furthermore, out of these two exclusive genes, one gene (Os03g15960) was also down-regulated in toerant genotype after 48 h late phase and another gene (Os10g02070) was exclusively down-regulated and up-regulated in CSV at 48 h late response phase and 24 h recovery phases, respectively (Additional file 4: Figure S4 and Additional file 5: Figure S5; Additional file 21: Table S15). These results indicated the fact that all 11 genes were significantly involved in different pathways and act as cold responsive genes.
Expression of genes (Os02g44230 and Os02g44235) encoding enzyme trehalose-6- phosphate phosphatase (TPP) were induced by cold response in CSV while reduced expression was observed in CTV at early response phase and 48 h late response phase. Down-regulation of these two genes was found after 24 h recovery condition of both in CSV and CTV genotypes. Genes (Os07g48160 and Os10g35110) encoding α-Gal were down-regulated at 48 h late response phase while up-regulated after 48 h recovery in CTV genotype. Further, gene encoding α-glucosidase (Os06g46284) and β-glucosidase (Os03g53800) were found down-regulated and up-regulated at 48 h late response phase and 24 h recovery in CTV genotype. In energy metabolism pathways, 24 and 11 DEGs were enriched in CTV and CSV at early response phase, respectively in which only 3 genes (Os04g56160, Os01g36720 and Os07g32570) were commonly down-regulated in both varieties. Furthermore, genes (Os04g16732, Os10g21394) encoding NAD (P) H-quiNone oxidoreductase subunit H (ndhH) were consistently up-regulated at each time interval of stress condition and down-regulated during 24 h recovery condition in CTV genotype. In addition, gene (Os07g42960) encoding 3-deoxy-7-phosphoheptulonate synthase, a chloroplast precursor involved in photosynthesis was significantly up-regulated in cold stress condition while down-regulated after 24 h recovery condition in CTV. Moreover, in our study, the rice gene (Os03g18070) encoding omega-3 fatty acid desaturase (desB) was significantly up-regulated during 48 h stress condition while down-regulated after 24 h recovery condition in CTV genotype.
Role of transcription factors under cold stress
Of the 2438 known or annotated TF genes in rice genome , 1583 differentially expressed TF genes were identified in which 453, 704 and 444 were present in both genotypes at 6 h early response phase, 48 h late response phase and 24 h recovery condition, respectively (Fig. 4; Additional file 6: Figure S6; Additional file 22: Table S16). In early response phase, a total of 192, 154 and 25 TF genes were detected as CSV specific, CTV specific, cold responsive genes, respectively. During late response phase of under the stress, a total of 239 CSV specific; 231 CTV specific and 233 commonly regulated TF genes were detected. While during recovery phase, a total of 158 CSV specific; 189 CTV specific and 93 commonly regulated TF genes were detected after the stress condition. Among the transcription factors, 69 WRKY, 46 bZIP, 41 NAC, 40 ERF, 31/14 MYB/MYB-related, 22 bHLH, 17 Nin-like 7 HSF and 4C3H were involved during 6 h early response phase cold stress. In 48 h late response phase, TF genes include 30 bHLH, 65 NAC, 30 ERF, 26/20 MYB/MYB-related, 11 C3H, 12 HSF, 86 Nin-like, 41 AP2/ERF, 55 bZIP and 98 WRKY members (Additional file 23: Table S17). Furthermore, 18 bHLH, 50 NAC, 31 ERF, 24/13 MYB/MYB-related, 4 C3H, 4 HSF, 14 Nin-like, 31 bZIP and 114 WRKY TF genes involved after 24 h recovery. Interestingly, 4 TF genes (LOC_Os01g72680, LOC_Os04g10010, LOC_Os07g46920 and LOC_Os09g25070) were down-regulated in susceptible genotype while up-regulated in tolerant genotype at 6 h early response phase (Additional file 24: Table S18). Moreover, 1 TF gene (Os09g10054) was up-regulated in tolerant while down-regulated in susceptible genotype at 48 h late response phase (Additional file 24: Table S18). Besides, 4 TF genes (Os09g25070; Os03g03164; Os07g46920; Os11g02540) were up-regulated in tolerant while down-regulated in susceptible genotype after 24 h recovery condition (Additional file 25: Table S19).
The results revealed 41 differentially-regulated AP2/ERF genes at low temperature stress, amongst them 15 were commonly detected in both genotypes. There were 8 (Os01g04020; Os06g47590; Os06g10780; Os02g43820; Os02g45420; Os02g48184; Os02g45450; Os07g13170) and 7 (Os03g60430; Os10g42150; Os04g46400; Os10g22600; Os02g13710; Os06g11860; Os06g09688) AP2/ERF genes exclusively induced and repressed in CSV genotype by chilling stress, respectively while 11 (Os05g41780; Os04g48350; Os02g42585; Os04g44670; Os09g20350; Os10g41130; Os05g41760; Os02g43970; Os04g52090; Os02g51670; Os04g57340) genes were repressed in CTV genotype. Also, a set of WRKY genes were found to be chilling-regulated in both genotypes, of which 12 were commonly induced, 22 were exclusively up-regulated in CSV and 18 were exclusively up-regulated in CTV genotype. In the rice genome, there are many MYB-encoding genes with diverse roles in developmental processes and defense responses. In our study, we identified 46 MYB and MYB-related TF genes showing differential expression to low temperature stress in both the genotypes. Out of the total MYB genes detected, 29 were commonly regulated while 17 and 11 were exclusively regulated in CSV and CTV genotype, respectively. Further, we identified three MYB TF genes in which two genes (Os02g38980; Os04g49450) were induced and one gene (Os02g46030) repressed in CTV genotype under cold stress.
Further, TF genes in response to cold stress were analyzed from transcription factor database (STIFDB V2.0; http://caps.ncbs.res.in/stifdb2/) for early response, late response and recovery phases. Results showed that 6 ERF (Os01g21120; Os01g73770; Os02g43790; Os06g03670; Os08g36920; Os09g28440) and 1 MYB (Os02g41510), 3 C2H2 (Os03g41390; Os03g60560; Os03g60570) gens were up-regulated in CSV at 6 h early response phase. Moreover, 2 ERF (Os09g35010; Os09g35030), 2 MYB (Os03g20090; Os04g43680), 2 Nin-like (Os01g68490; Os10g35640), 2 WRKY (Os02g08440, Os05g25770), 1 C2H2 (Os03g32230), 1 NAC (Os07g34280) and 1 MYB-related (Os04g49450) TF genes were commonly up-regulated in CSV and CTV, revealing a positive regulatory role in response to low temperature stress in rice. Besides this, 2 TFs (ERF, Os03g64260; WRKY, Os02g32230) genes were up-regulated in CTV only. At 48 h late response phase, one TF gene (Os01g14870) C3H and two Nin-like TF genes (Os01g68490; Os10g35640) were induced while down-regulated WRKY (Os01g60640), CO-like (Os02g49870), HSF (Os08g43334) and NAC (Os11g05614) in response to cold stress in CTV. In addition, comparative analysis of CSV genotype showed TF gene C3H (Os01g14870), Nin-like (Os01g68490), C2H2 (Os03g32230; Os03g60560; Os03g60570), CO-like (Os08g08120) were down-regulated while CO-like (Os02g49870) and HSF (Os08g43334) were down-regulated. After 24 h recovery conditions in CTV genotype, TF NF-YA (Os06g46799), CO-like (Os08g08120), NAC (Os11g05614) were down-regulated in response to cold stress while WRKY (Os02g32230) was up-regulated.
Quantitative real-time reverse transcription– PCR for validation of novel transcripts
The qRT-PCR analysis was performed for validating the results of RNA-seq sequencing as well as to test the differential expression of novel transcripts identified from the analysis. For the purpose the DEGs were selected carefully following certain criteria that (1) they must be highly up regulated or down regulated ones having at least 5-fold difference between the susceptible and tolerant genotype, (2) they should be differentially expressed being stage specific and/or genotype specific. On the basis of these criteria nine transcripts were considered for validation. The results of qRT-PCR for the genes LOC_Os02g44230.1, LOC_Os09g25070.1, LOC_Os09g10054.1, LOC_Os03g53800.1, LOC_Os03g18070.1, LOC_Os03g03164.1, LOC_Os05g42250.1, LOC_Os02g33030.1 and LOC_Os02g08440.1 normalized with housekeeping gene b-tubulin in response to different chilling treatments were validated. A positive correlation (R2 = 0.985) was observed between RNA-seq and qRT-PCR expression analyses (Fig. 5).
Phenotyping for chilling stress tolerance
There are few reports for classifying rice genotypes based on seedling stage cold tolerance under low temperature regimes and duration [1, 32]. Various classes of genotypes to cold tolerance at the seedling stage were evaluated at 10 °C for 10 and 13 days . Cold tolerance at seedling stage at 9 °C for 8, 14, 16 and 18 days were reported earlier . However, rice genotypes classification with respect to seedling stage cold tolerance by exposing the seedlings to a temperature regime of 25 to 4 °C and further extending the exposure up to 15, 20 and 25 days at 4 °C were reported in previous study where the rice genotypes were grouped into six classes [2, 3]. Rice genotype, Geetanjali was categorized as very highly tolerant and Sahabhagidhan as highly susceptible type to chilling stress reported by earlier publications [2, 3]. Thus, both the studied genotypes for RNA-seq analyses were confirmed to be very highly contrasting for seedling stage chilling stress tolerance.
RNA sequencing and DEGs in response to chilling stress
Previous study on seedling stage cold tolerance to identify underlying molecular genetics pathways imparting cold stress adaptation using microarray on affymetrix has been reported [25, 33]. The RNA-seq technology with an advantage of global representation and precise measurement of expression level of each gene in a sample by mapping short DNA sequences on a reference . A total of 24529 DEGs in both the contrasting genotypes were observed (Fig. 1). Previous findings on cold and salinity tolerance also revealed relatively large number of DEGs in contrasting genotypes [34,5,36]. On the basis of CuumeRbund hierarchical cluster analysis of transcripts clearly divided the response to cold stress into three phases: an early response phase (6 h cold stress), late response phase (12, 24 and 48 h cold stress) and recovery phase (after 24 h). These results were consistent with the previous reports on cold stress gene profiling [25, 33, 34]. Further, comparative analysis of DEGs of CSV and CTV revealed 219 and 116 DEGs commonly involved at each time interval, respectively in which 198 and 95 DEGs were exclusively present in CSV and CTV genotypes and 21 DEGs were commonly involved in both varieties (Fig. 1). When we compared the expression profile of DEGs of both the genotypes, significant numbers of DEGs were down-regulated as compared to up-regulated genes except 48 h recovery condition of CTV (Additional file 2: Figure S2; Additional file 9: Table S3).
Comparative gene expression level on the basis of FC revealed 24.97 and 75.03% genes up-regulated and down-regulated in CSV while 43.92 and 56.08% genes were up-regulated and down-regulated in CTV, respectively (Additional file 9: Table S3; Additional file 10: Table S4). This signifies that CTV harbors potential to specifically maintain the higher percentage of genes up-regulated during the course of cold stress. This suggests that presence of these genes may confer tolerance to CTV’s during exposure to cold stress. Similar result was observed for K354 and C418, during cold stress, the favorable alleles from indica donor Bg300 improved the cold tolerance of K354 over C418 . Our results also showed that more number of genes participated to recover CTV from cold stress. Further, in CTV, both increase and decrease in exclusively expressed genes were observed during various response phases, but in 48 h recovery condition a sharp decline was observed, indicating that plant regained the expression pattern similar to control condition during this short recovery period, while still maximum number of genes were involved in CSV during this phase also, indicating large number of genes were involved in its survival (Fig. 3; Additional file 2: Figure S2). For normal growth and development, 2207 and 4858 genes were exclusively expressed in CSV and CTV (Fig. 1). Thus, it may be concluded that more number of genes were involved in becoming the genotype as tolerant. Earlier studies on gene expression have suggested that the highly constitutive gene expression before exposure to abiotic stress might represent a constitutive tolerance in tolerant genotypes [37,38,39,40]. Many research findings have suggested that the major group comprises genes unresponsive to chilling stress in either genotype as like 154 genes that were heavily detected in CTV, such as those encoding glutathione stransferase, oxidoreductase and thioredoxin, those increased chilling tolerance by maintaining cell redox homeostasis [41,42,43,44].
Functional characterization of DEGs
In molecular function category, 3 GO common terms were enriched in both the genotypes. This indicated few common molecular functions in both CSV and CTV. Under chilling stress, the susceptible genotype ‘Sahabhagidhan’ was not enriched with catalytic activity, hydrolase activity and lipid binding molecular function. Hence, under stress, genotype that lacks these functions may be susceptible to the stress. The cellular component characterization showed six GO terms enriched of which 3 GO terms were significantly involved in all three phase of CSV while in CTV, these were enriched in 24 h recovery condition. This indicates the absence of photosynthesis related genes, storage products like starch and synthesis of other classes of molecules like fatty acids and terpenes used for energy production and as raw material for the synthesis of other molecules during cold stress condition in CSV. In biological process category, lipid metabolic process showed the same biological function in both susceptible and tolerant genotypes while signal transduction, regulation of cellular process, pollen-pistil interaction, regulation of biological process, photosynthesis, generation of precursor metabolites and energy, carbohydrate metabolic process were detected and enriched in early response phase of CTV. This indicates that significant numbers of DEGs are involved in genotype as tolerant during cold stress (Table 1; Additional file 12: Table S6).
In early response phase, a very high number of DEGs were enriched in biological and molecular process in CTV as compared to CSV. Therefore, Getanjali showed tolerance response to early phase cold stress. However, no DEGs were enriched in CTV for cellular component, while 272 DEGs were enriched in 3 categories. This indicated that in CTV, cold stress though changed the expression profiles of genes involved in molecular and biological process, the genes related to cellular component remained unchanged indicating inherent strong morphological features to with stand cold stress. In both CSV and CTV, 798 and 688 DEGs were represented for response to stress and stimulus, respectively. The DEGs effects were different in both the genotypes; in CTV it induced many DEGs under several categories except for cellular components, on contrary in CSV, no DEGs were represented for many of the important biological components: metabolic processes and signal transduction; molecular functions: catalytic activity, kinase activity, and protein binding. In cellular component, 102 DEGs of CSV were represented mainly for cell wall and all of them were down regulated (Fig. 2; Table 2; Additional file 3: Figure S3; Additional file 13: Table S7; Additional file 16: Table S10). During 12-48 h of cold stress (late response phase), similar to early response DEGs, more DEGs were enriched in response to stress and stimulus. In this response phase, cellular component showed enrichment in susceptible genotype. While, tolerant genotype was enriched in photosynthesis, generation of precursor metabolites and energy, thylakoid, transporter activity, transcription regulator activity, oxygen binding and transcription factor activity, along with new DEGs representing in photosynthesis, generation of precursor metabolites and energy, catalytic activity, kinase activity, biological regulation, signal transduction, transferase activity-transferring phosphorus-containing groups, regulation of cellular process, secondary metabolic process, regulation of biological process, lipid metabolic process, carbohydrate binding, and pollen-pistil interaction got stabilized during the late cold stress (Fig. 2; Additional file 14: Table S8; Additional file 17: Table S11). Thus, genotype Getanjali remain tolerant in response to chilling stress. Previous research finding reveals that DEGs are genotype specific and observed to be regulated in the late phase of chilling stress treatment . In recovery phase, the tolerant genotype showed more functional category enrichment than susceptible genotype particularly in response to stimulus and stress, metabolic process, lipid and carbohydrate metabolism and transcription factor/regulator activity, whereas in CSV, DEGs related to lipid binding, thylakoid, and plastid were exclusive (Additional file 15: Table S9; Additional file 18: Table S12). Hence, the genotype Getanjali is tolerant to chilling stress. Previous report revealed that commonly up-regulated genes during recovery conditions enriched in metabolic, oxidation-reduction, and stress-response processes, indicates their general role in recovery mechanisms of rice plants after chilling stress .
KEGG pathway annotation
It was found that cold treatment continuously increased the number of up-regulated and down-regulated genes in both CSV and CTV genotypes, while it decreased after 24 h recovery, except exclusively expressed genes of CTV during recovery (Table 2). These results indicate that prolonged cold stress (6-48 h) increased number of up-regulated and down-regulated genes that might be involved in tolerance/susceptibility response. Further, decrease in number of up-regulated and down-regulated genes after 24 h recovery indicated that plant adapted to normal environment and suppressed expression of stress related genes. Therefore, we analyzed early response phase, late response phase specifically at 48 h cold stress and 24 h recovery to reveal the role of genes during cold stress and after recovery, normalized with control condition in both the genotypes. Comparative analysis of two genotypes revealed 29216 and 30867 genes expressed in susceptible and tolerant genotypes, respectively under control condition. These results indicate that the CTV contain higher DEGs in different pathways by which the genotype tolerate better under the stress condition. The commonly expressed genes (26009) in both varieties may be involved in plant growth and development. Comparison of early response phase in both varieties indicated that around 46% of DEGs were up-regulated and 57% DEGs were down regulated in CTV while 20% DEGs were up-regulated and 80% DEGs were down-regulated in CSV. These results indicate that tolerant genotype contains high percentage of up-regulated genes than susceptible genotype making the genotype tolerant.
Earlier research findings have suggested that cold responsive genes encode proteins involved in enzymes for respiration and metabolism of carbohydrates, lipids, phenylpropanoids and antioxidants, molecular chaperones and antifreeze proteins [45, 46]. Trehalose acts against stress protection of metabolite and for storage of carbohydrate. In plants, trehalose biosynthesis is catalyzed by trehalose-6-phosphate synthase (TPS) and trehalose- 6- phosphate phosphatase (TPP). Earlier finding suggests that TPP genes express in rice and their expression is influenced by cold stress . Chilling stress accumulates trehalose rapidly and transiently due to transient induction of TPP activity in rice tissues. Over expression of TPS and TPP genes increased the level of trehalose accumulation and showed more tolerance to cold stress in transgenic tobacco and rice [48,49,50,51]. Our results showed that expression of genes (Os02g44230 and Os02g44235) encoding enzyme TPP were induced by cold response in CSV while reduced expression was observed in CTV at early response phase and 48 h late response phase. Down-regulation of these two genes were found after 24 h recovery condition of both CSV and CTV, suggesting the fact that cold tolerant plant synthesize increased amount of TPP/TPS for regulating cell shape and plant architecture. In addition, genes (Os02g51680 and Os09g20390) encoding enzyme TPP/TPS were down-regulated in CSV. In transgenic petunia, down-regulation of α-Gal (α-Galactosidase) exhibited more tolerance under chilling stress suggested for engineering of raffinose metabolism through transformation with α–Gal . Our results showed that genes (Os07g48160 and Os10g35110) encoding α-Gal were down-regulated at 48 h late response phase while up-regulated after 48 h recovery in CTV. Further, gene encoding α-glucosidase (Os06g46284) and β-glucosidase (Os03g53800) were found down-regulated and up-regulated at 48 h late response phase and 24 h recovery in CTV.
In energy metabolism pathways, 24 and 11 DEGs were enriched in CTV and CSV at early response phase, respectively in which only 3 genes (Os04g56160, Os01g36720 and Os07g32570) were commonly down-regulated in both varieties. Results indicated that a large proportion of genes were exclusively expressed and involved in different pathways of CTV, which were not participated in the pathway of CSV and made a substantial difference between tolerant and susceptible types in early response phase. Furthermore, genes (Os04g16732, Os10g21394) encoding NAD (P) H-quiNone oxidoreductase subunit H (ndhH) were consistently up-regulated at each time interval of stress condition and down-regulated during 24 h recovery condition in CTV. These results indicated that NDH-dependent cyclic pathway around PSI participates in ATP supply in conditions of high ATP demand in response to severe stress conditions. In addition, gene (Os07g42960) encoding 3-deoxy-7-phosphoheptulonate synthase, a chloroplast precursor involved in photosynthesis was significantly up-regulated in cold stress condition while down-regulated after 24 h recovery condition in CTV. Moreover, in our study, the rice gene (Os03g18070) encoding omega-3 fatty acid desaturase (desB) was significantly up-regulated during 48 h stress condition while down-regulated after 24 h recovery condition in CTV indicating the fact that increased linolenic acid levels are essential for the maintenance of membrane fluidity and chloroplast function under chilling stress exposure . Cold stress also leads to change in the expression profile of enzyme peroxidase (POD) and their activity increases under low temperature . Our results showed that gene encoding POD (Os01g15830) was significantly down-regulated at 6 h, 24 h, 48 h and 24 h recovery conditions in CSV while up-regulated at 24 h and after 24 h recovery conditions in CTV. In regard to catalase (CAT) activity, gene expression profile of both varieties was observed differently. Gene encoding catalase (Os06g51150 and Os02g02400) was significantly up-regulated after 24 h recovery condition in CTV while down-regulated during 48 h stress condition in CSV. Previous results also suggested different expression profile of gene encoding CAT in chilling stress microarray study of rice .
We identified gene, Os08g06344 encoding ATP-dependent RNA helicase, that acts on acid anhydrides to facilitate cellular and sub cellular movement, was up-regulated at 6 h and 48 h late response phase in CTV and may be responsible for cold tolerance in rice. Low expression of osmotically responsive gene (LOS4), one of the helicase encoded by the Arabidopsis, is responsible for tolerance to chilling and freezing stress . Defects in the nucleocytoplasmic transport of RNA seem to affect cold tolerance preferentially, because the LOS4 mutant plants do not have severe growth or developmental phenotypes, nor they are strongly altered in other abiotic stresses. Further, we indentified 12 differentially expressed genotype specific genes in which 4 DEGs (Os08g44360; Os02g54140; Os02g07110; Os05g42250) were exclusively present in CTV and 8 DEGs (Os05g03610; Os02g33030; Os10g34480; Os01g39860; Os09g04050; Os02g09490; Os06g14240; Os05g03610) were exclusively present in CSV. The roles of these genes need to be further examined.
Role of transcription factors under cold stress
A large group of transcription factors are involved in response to abiotic stresses comprising the AP2/EREBP family [31, 56]. In our study, we detected 41 differentially-regulated AP2/ERF genes at low temperature stress, amongst them 15 were commonly detected in both varieties. These results are in conformity with the previous reports on low temperature stress regulation due to presence of CBF/DREB subfamily [57,58,59]. More than half of the 103 WRKY TF genes have been identified in the rice varieties that are shown to be transcriptionally regulated under conditions of biotic and/or abiotic stresses [60, 61]. In our study, a set of WRKY genes were found to be chilling-regulated in both genotypes, of which 12 were commonly induced, 22 were exclusively up-regulated in CSV and 18 were exclusively up-regulated in CTV. Three chilling-induced WRKY genes (Os05g0343400, Os01g0246700, and Os01g0826400) in both varieties were also observed to be induced by chilling stress in a japonica rice genotype . The remaining WRKY TF genes were differentially regulated involved in differential response of the varieties.
In our study, we identified 46 MYB and MYB-related TF genes showing differential expression to low temperature stress in both the varieties. Out of the total MYB genes detected, 29 were commonly regulated while 17 and 11 were exclusively regulated in CSV and CTV, respectively. Eleven of these MYB/MYB-related genes have been reported in japonica rice in response to chilling stress . There are many evidences in arabidopsis and rice to support for the involvement of many MYB proteins in response to abiotic stress particularly to seedling stage cold stress [21, 62,63,64]. Therefore, our detected MYB TF genes are functionally involved in response to chilling stress. We found that OsMYB4 (Os04g43680) was induced in response to chilling stress in CSV. Earlier finding also suggest that Myb4 over expressing plants show increased cold tolerance . Further, we identified three MYB TF genes in which two genes (Os02g38980; Os04g49450) were induced and one gene (Os02g46030) repressed in CTV under cold stress.
Further, TF genes in response to cold stress were analyzed from transcription factor database (STIFDB V2.0; http://caps.ncbs.res.in/stifdb2/) for early response, late response and recovery phases. Our results showed that 6 ERF and 1 MYB, 3 C2H2 gens were up-regulated in CSV at 6 h early response phase. Moreover, 2 ERF, 2 MYB, 2 Nin-like, 2 WRKY, 1 C2H2, 1 NAC and 1 MYB-related TF genes were commonly up-regulated in CSV and CTV, revealing a positive regulatory roles in response to low temperature stress in rice. Besides this, ERF (Os03g64260) and WRKY (Os02g32230) genes were up-regulated in CTV only. Earlier report indicated that TFs (AP2/EREBP, MYB, HSF, and NAC) regulated most of the genes after 2 h chilling stress . At 48 h late response phase, one TF gene (Os01g14870) C3H and two Nin-like TF genes (Os01g68490; Os10g35640) were induced while down-regulated WRKY (Os01g60640), CO-like (Os02g49870), HSF (Os08g43334) and NAC (Os11g05614) in response to cold stress in CTV. In addition, comparative analysis of CSV showed TF gene C3H (Os01g14870), Nin-like (Os01g68490), C2H2 (Os03g32230; Os03g60560; Os03g60570), CO-like (Os08g08120) were down-regulated while CO-like (Os02g49870) and HSF (Os08g43334) were down-regulated. After 24 h recovery conditions in CTV, TF NF-YA (Os06g46799), CO-like (Os08g08120), NAC (Os11g05614) were down-regulated in response to cold stress while WRKY (Os02g32230) was up-regulated. In CSV, C2H2 (Os03g60560; Os03g60570) was up-regulated while MYB-related (Os02g46030) and NF-YA (Os06g46799) were down-regulated [65,66,67].
Quantitative real-time reverse transcription– PCR for validation of novel transcripts
For validating the results of RNA-seq sequencing, qRT-PCR was performed to test the differential expression of novel transcripts identified from analysis. The results of qRT-PCR for the nine genes normalized with housekeeping gene b-tubulin in response to different chilling treatments were validated. A positive correlation (R2 = 0.985) was observed between RNA-seq and qRT-PCR expression analyses (Fig. 5). This confirmed the reliability of expression of differentially expressed genes identified through RNA-seq.
The down regulated genes were more in susceptible genotype in all 3 phases while up regulated genes were more in tolerant genotype. The down-regulated genes declined during recovery in CSV and CTV. The numbers of up-regulated genes were more in CTV during recovery phase under cold stress. A significant involvement of transcription regulation, oxygen and lipid binding, catalytic and hydrolase activity was established in gene ontology classification. Absence of photosynthesis related genes, storage products like starch and synthesis of other classes of molecules like fatty acids and terpenes used for energy production and as raw material for the synthesis of other molecules during cold stress condition were observed in cold susceptible varieties. During recovery phase of cold stress, no genes for photosynthesis and storage products like starch and the synthesis of many classes of molecules are detected in the CSV. Under biological process category, response to biotic stimulus, response to abiotic stimulus, secondary metabolic process and lipid metabolic process were influenced for both the varieties. Catalytic activity, hydrolase activity and lipid binding molecular function are absent in susceptible genotype. More number of DEGs were involved in different pathways of tolerant genotype imparting the genotype as tolerant to chilling stress. A total of 41 differentially-regulated AP2/ERF transcription factors genes were detected at low temperature stress. A set of WRKY genes comprising 12 were commonly induced, 22 were exclusively up-regulated in CSV and 18 were exclusively up-regulated in CTV. Also, 46 MYB and MYB-related TF genes were identified showing differential expression to low temperature stress in both the varieties.
Plant materials and growth conditions
The highly tolerant genotype Geetanjali and susceptible cultivar, Sahabhagidhan to vegetative stage cold tolerance were collected from ICAR-National Rice Research Institute, Cuttack, Odisha, India and raised in pots till 21 days. The potted plants were exposed to cold stress condition in a walk-in cold growth chamber (Conviron PGV36). In the growth chamber, 75–85% relative humidity (RH), 800 μ moles s− 1 m− 2 light intensity above 60 cm from the floor and 12 h photoperiod were maintained. Phenotyping for seedling stage cold tolerance in the contrasting varieties was performed as per the published protocol [2, 3]. In the protocol, seedlings were grown in RGA-cum-Phytotron till three-leaf stage at a set temperature of 25 °C and photoperiod of 12 h. During this stage, the weak seedlings were removed and healthy seedlings were exposed to low temperature regime (LTR) treatments in the growth chamber with the pre setting of temperature and RH as in the protocol. The low temperature regime1 (LTR1) was applied by exposing the seedling to a temperature of 25 °C exposed for 3 days; low temperature regime 2 (LTR2) was achieved by a gradual decrease in temperature to 15 °C within 7 days from 25 °C and maintained at 15 °C; low temperature regime 3 (LTR3) exposed to gradual decrease in temperature to 8 °C from 15 °C within 7 days and maintained at 8 °C; low temperature regime 4 (LTR4) with a gradual decrease in temperature to 4 °C from 8 °C within 7 days and maintained at 4 °C; low temperature regime 5 (LTR5) exposed to temperature of 4 °C for 14 days and low temperature regime 6 (LTR6) to temperature of 4 °C for 21 days. The varieties were evaluated using a modified IRRI-SES score under control facility.
Leaf samples of both the varieties were collected at 25 °C, after 6 h, 12 h, 24 h, 48 h at 4 °C exposure and recovery was observed after 48 h at 4 °C and kept for 24 h at 25 °C. The leaf samples from three different biological replicates were collected in different vials containing RNA stabilizer solution (Xcelris, India) and immediately stored at − 80 °C. 1–2 leaves from each seedling were collected in order to include maximum number of plants in single biological replicate. The three biological replicates included three separate set of chilling stress treatment experiments. For total mRNA isolation 100 mg of tissue from each biological replicates were taken together and crushed with liquid nitrogen to make powder form. 100 mg of the this powdered tissue was taken for further mRNA isolation followed by RNAseq analysis. The treatments for RNA-seq analyses were S0 to S5 and T0 to T5 for CSV and CTV at 0, 6, 12, 24, 48 and 24 h recovery after 48 h cold stress conditions, respectively. Zero hour stress indicates normal/control condition (S0/T0) whereas S1 to S5 and T1 to T5 represent the treatments.
RNA isolation and sequencing
The leaf samples from three different biological replicates were collected in different vials containing RNA stabilizing solution. 1–2 leaves from each seedling were collected in order to include maximum number of plants in single biological replicate. The three biological replicates included three separate set of chilling stress treatment experiments. For mRNA isolation 100 mg of tissue from each biological replicates were taken together and crushed with liquid nitrogen to make powder form. 100 mg of the this powdered tissue was taken for further RNA isolation followed by RNAseq analysis. Total mRNA was isolated XcelGen total RNA isolation kit (Xcelris Genomics, India) as per manufacturer’s instructions. The yield and purity of RNA were evaluated by recording absorbance at 260 and 280 nm on Nanodrop 8000 Spectrophotometer (Thermoscientific). Further, the integrity of RNA was checked using RNA 6000 Nano LabChip using Agilent Bioanalyzer 2100 (Agilent Technologies, Germany). The pair end sequencing libraries for RNA-seq were prepared using illumina TruSeq® RNA Library Preparation Kit as per manufacturer’s protocol (illumina®, San Diego, CA). Library quality control and quantification were performed on Caliper LabChip GX using HT DNA High Sensitivity Assay Kit. RNA-Seq was performed using illumina HiSeq2000 (illumina®, San Diego, CA) as per manufacturer’s protocol to generate 30 million 2 × 100 bp reads for each sample. The processing of fluorescent images into sequences, base-calling and quality value calculations were performed using the illumina data processing pipeline (version 1.8).
Reads mapping and annotations
Rice genome and gene information for reference cultivar, Nipponbare (Oryza sativa L. subsp. japonica) was downloaded from ftp:// ftp.plantbiology.msu.edu/pub/data/EukaryoticProjects/osativa/annotationdbs/pseudomolecules/version_7.0/all.dir/. High quality (QV > 25) transcriptome libraries from all 12 samples were individually mapped to the rice genome using TopHat v1.3.3 and Bowtie v0.12.9 with default parameters. Transcript abundance and differential gene expression were estimated using the programme Cufflinks v1.3.0 alongwith FPKM (Fragments Per Kilobase of transcript per Million mapped reads) from reference-guided mapping . Genes expressed at very low levels (read counts < 10 across all 12 libraries) were excluded from differential gene expression analysis. Significance tests for differential expression were based on a modified exact test. Differentially expressed genes were identified at 0.05 false discovery rate (FDR) 0.05 and p-value ≤0.05. The CuumeRbund, a package of R was used to produce the heat maps of differentially expressed genes .
Gene ontology (GO), KEGG Orthology (KO), and enrichment analysis
Functional classification of the DEGs was performed using agriGo; GO analysis toolkit and database using singular enrichment analysis (SEA) against supported species Oryza sativa MSU build 7.0 and further annotations were plotted using Web Gene Ontology Annotation Plot (WEGO) software. Statistical methods like Hypergeometric was used to obtain significant GO-terms considering P-value < 0.05. Hypergeometric tests with Benjamini and Hochberg FDRs were used for analysis taking the default parameters to adjust the P-value. Functional annotation of differentially expressed genes was done by BLAST comparisons against KEGG GENES database KAAS (KEGG Automatic Annotation Server). The KO terms were assigned (representative gene data set for rice Oryza sativa Ref.Seq) by using the SBH (Single-directional best hit) option. Pathway mapping was done using the KEGG Orthology database (http://www.genome. jp/kegg/ko.html) .
Transcription factor identification
For the identification of transcription factor families represented in rice transcriptome, 10639 DEG were searched against 2438 transcription factors of Oryza sativa subsp. japonica, classified into 56 families at rice transcription factor database (PlantTFDB v2.0) using BLASTX with an E-value cut-off of 1E-05.
Quantitative real-time reverse transcription–PCR for validation of novel transcripts
The expression analysis for validation of novel transcripts involved in chilling stress, the transcript analysis for genes LOC_Os02g44230.1, LOC_Os09g25070.1, LOC_Os09g10054.1, LOC_Os03g53800.1, LOC_Os03g18070.1, LOC_Os03g03164.1, LOC_Os05g42250.1, LOC_Os02g33030.1 and LOC_Os02g08440.1 performed with real-time method. Twenty one days old rice seedlings were exposed to temperature at 25 °C, after 6 h, 12 h, 24 h, 48 h at 4 °C exposure and recovery after 48 h at 4 °C and kept for 24 h at 25 °C. Total mRNA was isolated from the stressed and non-stressed plants with biological replicates using Trizol LS reagent (Qiagen) as per the manufacturer’s instructions. cDNA was synthesized by using oligo-dT primers with SuperscriptII (Invitrogen) cDNA synthesis kit. Expression analysis of Os02g44230.1, Os09g25070.1, Os09g10054.1, Os03g53800.1, Os03g18070.1, Os03g03164.1, Os05g42250.1, Os02g33030.1 and Os02g08440.1 in response to duration of stress was performed by real-time analysis with real-time primers (Additional file 26: Table S20). β-tubulin gene from rice was taken for normalization. Fold change in expression of the genes in various hours of cold stress conditions as compared to controlled condition were calculated by using 2-ΔΔ Ct method . Final expression result was obtained from three independent biological replicates and three technical replicates.
Availability of data and materials
Extra data has been appended as supplementary Tables. The accession number for sequence data generated in this study is PRJNA288892 available at http://trace.ncbi.nlm.nih.gov/.
ABA responsive element
Cold stress binding factors
Cold susceptible variety
Cold tolerant variety
Differentially expressed genes
Dehydration responsive element binding proteins
Ethylene responsive element
False discovery rate
Fragments per kilo base of transcript per million mapped reads
Heat shock factor
Kyoto encyclopedia of genes and genomes
Low temperature regime
Nicotinamide adenine dinucleotide
Quantitative real time reverse transcription polymerase chain reaction
Reactive oxygen species
Standard evaluation system
Web gene ontology annotation plot
Andaya VC, Mackill DJ. Mapping of QTLs associated with cold tolerance during the vegetative stage in rice. J Exp Bot. 2003;54:2579–85.
Pradhan SK, Nayak DK, Guru M, Pandit E, Das S, Barik SR, Mohanty SP, Anandan A. Screening and classification of genotypes for seedling stage chillingstress tolerance in rice and validation of the trait using SSR markers. Plant Genet Resour; Charact Util. 2015;15(46). https://doi.org/10.1017/S1479262115000192.
Pandit E, Tasleem S, Barik SR, Mohanty DP, Nayak DK, Mohanty SP, Das S, Pradhan SK. Genome-wide association mapping reveals multiple qtls governing tolerance response for seedling stage chilling stress in indica rice. Front Plant Sci. 2017;8:552. https://doi.org/10.3389/fpls.2017.00552.
Das S, Pandit E, Guru M, Nayak DK, Tasleem S, Barik SR, Mohanty DP, Mohanty SP, Patra BC, Pradhan SK. Genetic diversity, population structure, marker validation and kinship analysis for seedling stage cold tolerance in indica rice. Oryza. 2018;55:396–405.
Sthapit BR, Witcombe JR. Inheritance of tolerance to chilling stress in rice during germination and plumule greening. Crop Sci. 1998;38:660–5.
Jiang L, Liu SJ, Hou MY, Tang JY, Chen LM, Zhai HQ, et al. Analysis of QTLs for seed low temperature germinability and anoxia germinability in rice (Oryza sativa L.). Field Crop Res. 2006;98:68–75. https://doi.org/10.1016/j.fcr.2005.12.015.
Pradhan SK, Barik SR, Sahoo A, Mohapatra S, Nayak DK, Mahender A, Meher J, Anandan A, Pandit E. Population structure, genetic diversity and molecular marker-trait association analysis for high temperature stress tolerance in rice. PLoS One. 2016;11(8):e0160027. https://doi.org/10.1371/journal.pone.0160027.
Mittler R. Abiotic stress, the field environment and stress combination. Trends Plant Sci. 2006;11:15–9.
Moller IM. Plant mitochondria and oxidative stress: Electron transport, NADPH turnover, and metabolism of reactive oxygen species. Annu Rev Plant Physiol Plant Mol Biol. 2001;52:561–91.
Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends Plant Sci. 2002;7:405–10.
Miller G, Shulaev V, Mittler R. Reactive oxygen signaling and abiotic stress. Physiol Plant. 2008;133:481–9.
Miller G, Shulaev V, Mittler R. Redox regulation in photosynthetic organisms: signaling, acclimation, and practical implications. Antioxid Redox Signal. 2008b;11:861–905.
Chinnusamy V, Zhu J, Zhu JK. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12:444–51.
Zhou MQ, Shen C, Wu LH, Tang KX, Lin J. CBF-dependent signaling pathway: a key responder to low temperature stress in plants. Crit Rev Biotechnol. 2011;31:186–92.
Cheng C, Yun KY, Ressom HW, Mohanty B, Bajic VB, Jia Y, Yun SJ, de los Reyes BG. An early response regulatory cluster induced by low temperature and hydrogen peroxide in seedlings of chilling-tolerant japonica rice. BMC Genomics. 2007;8:175.
Fujita M, Fujita Y, Noutoshi Y, Takahashi F, Narusaka Y, Yamaguchi-Shinozaki K, Shinozaki K. Crosstalk between abiotic and biotic stress responses: a current view from the points of convergence in the stress signaling networks. Curr Opin Plant Biol. 2006;9:436–42.
Nakashima K, Ito Y, Yamaguchi-Shinozaki K. Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 2009;149:88–95.
Mittal D, Chakrabarti S, Sarkar A, Singh A, Grover A. Heat shock factor genefamily in rice: genomic organization and transcript expression profiling in response to high temperature, low temperature and oxidative stresses. Plant Physiol Biochem. 2009;47:785–95.
Kant P, Gordon M, Kant S, Zolla G, Davydov O, Heimer YM, Chalifa-Caspi V, Shaked R, Barak S. Functional-genomics-based identification of genes that regulate Arabidopsis responses to multiple abioticstresses. Plant Cell Environ. 2008;31:697–714.
Ma TL, Wu WH, Wang Y. Transcriptome analysis of rice root responses to potassium deficiency. BMC Plant Biol. 2012;12:161.
Ishimaru Y, Suzuki M, Masuda H, Takahashi M, Nakanishi H, Mori S. OsZIP4, a novel zinc-regulated zinc transporter in rice. J Exp Bot. 2005;56:3207–14.
Wei G, Tao Y, Liu G, Chen C, Luo R, XiaH GQ, Zeng H, Lu Z, Han Y, Li X, Song G, Zhai H, Peng Y, Li D, Xu H, Wei X, Cao M, Deng H, Xin Y, Fu X, Yuan L, Yu J, Zhu Z, Zhu L. A transcriptomic analysis of super hybrid rice LYP9 and its parents. Proc Natl Acad Sci. 2009;106:7695–701.
Mittal D, Madhyastha DA, Grover A. Genome-wide transcriptional profiles duringtemperature and oxidative stress reveal coordinated expression patterns and overlapping regulons in rice. PLoS One. 2012;7:e40899.
Dubey S, Misra P, Dwivedi S, Chatterjee S, Bag SK, Mantri S, Mehar HA, Rai A, Kumar S, Shri M, Tripathi P, Tripathi RD, Trivedi PK, Chakrabarty D, Tuli R. Transcriptomic and metabolomic shifts in rice roots in response to Cr (VI) stress. BMC Genomics. 2010;11:648.
Zhang T, Zhao X, Wang W, Pan Y, Huang L, Liu X, Zong Y, Zhu L, Yang D, Fu B. Comparative transcriptome profiling of colding stress responsiveness in two contrasting rice genotypes. PloSone. 2012;7:e43274.
Lu T, Lu G, Fan D, Zhu C, Li W, Zhao Q, Feng Q, Zhao Y, Guo Y, Li W, Huang X, Han B. Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-seq. Genome Res. 2010;20:1238–49.
Zhai R, Feng Y, Wang H, Zhan X, Shen X, Wu W, Zhang Y, Chen D, Dai G, Yang Z, Cao L, Cheng S. Transcriptome analysis of rice root heterosis by RNA-Seq. BMC Genomics. 2013;14:19.
He G, Zhu X, Elling AA, Chen L, Wang X, Guo L, Liang M, He H, Zhang H, Chen F, Qi Y, Chen R, Deng XW. Global epigenetic and transcriptional trends among two rice subspecies and their reciprocal hybrids. Plant Cell. 2010;22:17–33.
Oono Y, Kawahara Y, Kanamori H, Mizuno H, Yamagata H, Yamamoto H, Hosokawa S, Ikawa H, Akahane I, Zhu Z, Wu J, Itoh T, Matsumoto T. mRNA-Seq reveals a comprehensive transcriptome profile of Rice under phosphate stress. Rice. 2011;4:50–65.
Hu Y, Zhang L, He S, Huang M, Tan J, Zhao L, Yan S, Li H, Zhou K, Liang Y, et al. Cold stress selectively unsilences tandem repeats in heterochromatin associated with accumulation of H3K9ac. Plant Cell Environ. 2012;35:2130–42.
Gao G, Zhong Y, Guo A, Zhu Q, Tang W, Zheng W, Gu X, Wei L, Luo J. DRTF: a database of rice transcription factors. Bioinformatics. 2006;22(10):1286–7.
Zhang ZH, Su L, Chen W, Li W, Zhu YG. A major QTL conferring cold tolerance at early seedling stage using recombinant inbred lines of rice (Oryza sativa L.). Plant Sci. 2005;168:527–53.
Zhan F, Huang L, Wang W, Zhao X, Zhu L, Fu B, Li Z. Genome-wide gene expression profiling of introgressed indica rice alleles associated with seedling cold tolerance improvement in a japonica rice background. BMC Genomics. 2012;13:461.
Socquet-Juglard D, Kamber T, Pothier JF, Christen D, Gessler C, Duffy B, Patocchi A. Comparative RNA-Seq Analysis of Early-Infected Peach Leaves by the Invasive Phytopathogen Xanthomonas arboricola pv. pruni. PLoS ONE. 2013;8(1):e54196. https://doi.org/10.1371/journal.pone 0054196.
Wang H, Wang H, Shao H, Tang X. Recent advances in utilizing transcription factors to improve plant abiotic stress tolerance by transgenic technology. Front Plant Sci. 2016;7:67. https://doi.org/10.3389/fpls.2016.00067.
Walia H, Wilson C, Condamine P, Liu X, Ismail AM, Zeng LH, Wanamaker SI, Mandal J, Xu J, Cui XP, Close TJ. Comparative transcriptional profiling of two contrasting rice genotypes under salinity stress during the vegetative growth stage. Plant Physiol. 2005;139(2):822–35.
Taji T, Seki M, Satou M, Sakurai T, Kobayashi M, Ishiyama K, Narusaka Y, Narusaka M, Zhu JK, Shinozaki K. Comparative genomics in salt tolerance between Arabidopsis and Arabidopsis-related halophyte salt cress using Arabidopsis microarray. Plant Physiol. 2004;135:1697–709.
Frank G, Pressman E, Ophir R, Althan L, Shaked R, FreedmanM SS, Firon N. Transcriptional profiling of maturing tomato (Solanum lycopersicum L.) microspores reveals the involvement of heat shock proteins, ROS scavengers, hormones, and sugars in the heat stress response. J Exp Bot. 2009;60:3891–908.
Kumar K, Kumar KM, HSR R, Cho YG. Insights into genomics of salt stress response in rice. Rice. 2013;6:27. https://doi.org/10.1186/1939-8433-6-27.
Bita CE, Zenoni S, Vriezen WH, Mariani C, Pezzotti M, Gerats T. Temperature stress differentially modulates transcription in meiotic anthers of heat-tolerant and heat-sensitive tomato plants. BMC Genomics. 2011;12:384.
Binh LT, Oono K. Molecular cloning and characterization of genes related to chilling tolerance in rice. Plant Physiol. 1992;99:1146–50.
Fowler S, Thomashow MF. Arabidopsis transcriptome profiling indicatesthat multiple regulatory pathways are activated during cold acclimation inaddition to the CBF cold response pathway. Plant Cell. 2002;14:1675–90.
Soranzo N, SariGorla M, Mizzi L, De-Toma G, Frova C. Organizationand structural evolution of the rice glutathione S-transferase gene family. Mol Gen Genomics. 2004;271:511–21.
Kavanagh KL, Jörnvall H, Persson B, Oppermann U. Medium- andshort-chaindehydrogenase/reductase gene and protein families: the SDRsuperfamily: functional and structural diversity within a family of metabolic and regulatory enzymes. Cell Mol Life Sci. 2008;65:3895–906.
Mahapatra SS, Wolfraim L, Poole RJ, Dhindsa RS. Molecular cloning and relationship to freezing tolerance of cold acclimation- specific genes of alfalfa. Plant Physiol. 1989;89:375–80.
Guy CL. Cold acclimation and freezing stress tolerance: role of protein metabolism. Annu Rev Plant Physiol Plant Mol Biol. 1990;41:187–223.
Pramanik MH, Imai R. Functional identification of a trehalose 6-phosphate phosphatasegene that is involved in transient induction of trehalose biosynthesis during chilling stress inrice. Plant Mol Biol. 2005;58:751–62.
Garg AK, Kim JK, Owens TG, Ranwala AP, Choi YD, Kochian LV, Wu RJ. Trehalose accumulation in rice plants confers high tolerance levels to different abiotic stresses. Proc Natl Acad Sci. 2002;99:15898–903.
Jang IC, Oh SJ, Seo JS, Choi WB, Song SI, Kim CH, Kim YS, Seo HS, Choi YD, Nahm BH. Expression of a bifunctional fusion of the Escherichia coli genes fortrehalose-6- phosphate synthase and trehalose- 6-phosphate phosphatase in transgenicriceplants increases trehalose accumulation and abiotic stress tolerance without stunting growth. Plant Physiol. 2003;131:516–24.
Ge LF, Chao DY, Shi M, Zhu MZ, Gao JP, Lin HX. Overexpression of the trehalose-6 phosphate phosphatase gene OsTPP1 confers stress tolerance in rice and results in the activation of stress responsive genes. Planta. 2008;228:191–201.
Iordachescu M, Imai R. Trehalose biosynthesis in response to abiotic stresses. J Integr Plant Biol. 2008;50:1223–9.
Pennycooke JC, Jones ML, Stushnoff C. Down-regulating α-galactosidase enhances freezing tolerance in transgenic Petunia. Plant Physiol. 2003;133:901–9.
Routaboul JM, Fischer SF, Browse J. Trienoic fatty acids are required to maintain chloroplast function at low temperatures. Plant Physiol. 2000;124(4):1697–705.
Wenying L, Kenming Y, Tengfei H, Feifei L, Dongxu Z, Jianxia L. The temperature induced physiological responses of Avene nuda L. A cold tolerant plant species. Sci World J. 2013; Article ID 658793. https://doi.org/10.1155/2013/658793.
Gong Z, Dong CH, Lee H, Zhu J, Xiong L, Gong D, Stevenson B, Zhu JK. ADEAD box RNA helicase is essential formRNA export and important for development and stress responses in Arabidopsis. Plant Cell. 2005;17:256–67.
Sharoni AM, Nuruzzaman M, Satoh K, Shimizu T, Kondoh H, Sasaya T, Choi R, Omura T, Kikuchi S. Gene structures, classification and expression models of the AP2/EREBP transcription factor family in rice. Plant Cell Physiol. 2011;52:344–60.
Yun KY, 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:16.
Dubouzet JG, Sakuma Y, Ito Y, Kasuga M, Dubouzet EG, Miura S, Seki M, Shinozaki K, Yamaguchi-Shinozaki K. OsDREB genes in rice, Oryzasativa L., encode transcription activators that function in drought-, high-salt- and cold-responsive gene expression. Plant J. 2003;33:751–63.
Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 1819;2012:86–96.
Ramamoorthy R, Jiang SY, 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:865–79.
Berri S, Abbruscato P, Faivre-Rampant O, Brasileiro AC, Fumasoni I, Satoh K, Kikuchi S, Mizzi L, Morandini P, Pè ME, Piffanelli P. Characterization of WRKY co-regulatory networks in rice and Arabidopsis. BMC Plant Biol. 2009;9:120.
Vannini C, Locatelli F, Bracale M, Magnani E, Marsoni M, Osnato M, Mattana M, Baldoni E, Coraggio I. Overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis thaliana plants. Plant J. 2004;37:115–27.
Dai X, Xu Y, Ma Q, Xu W, Wang T, Xue Y, Chong K. Overexpression of an R1R2R3 MYB gene, OsMYB3R-2, increases tolerance to freezing, drought, and salt stress in transgenic Arabidopsis. Plant Physiol. 2007;143:1739–51.
Su CF, Wang YC, Hsieh TH, Lu CA, Tseng TH, Yu SM. A novel MYBS3-dependent pathway confers cold tolerance in rice. Plant Physiol. 2010;153:145–58.
Mahantesha NK, Shameer OK, Mathew RG, Ramanathan S. STIFDB2: an updated version of plant stress-responsive transcription factor Data Base with additional stress signals, Stress-Responsive Transcription Factor Binding Sites and Stress-Responsive Genes in Arabidopsis and Rice. Plant and Cell Physiol. 2013;54:e8. https://doi.org/10.1093/pcp/pcs185.
Shameer K, Ambika S, Varghese SM, Karaba N, Udayakumar M, Sowdhamini R. STIFDB-Arabidopsis stress responsive transcription factor Data Base. Int J Plant Genomics. 2009;2009:583429. https://www.ncbi.nlm.nih.gov/pubmed/19841686.
Ambika S, Varghese SM, Shameer K, Udayakumar M, Sowdhamini R. STIF: hidden Markov model-based search algorithm for the recognition of binding sites of stress upregulated transcription factors and genes in Arabidopsis thaliana. Bioinformation. 2008;2(10):431–7.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.
Thumma BR, Sharma N, Southerton SG. Transcriptome sequencing of Eucalyptus camaldulensis seedlings subjected to water stress reveals functional single nucleotide polymorphisms and genes under selection. BMC Genomics. 2012;1(13):364.
Kyndt T, Denil S, Haegeman A, Trooskens G, De-Meyer T, Criekinge VW, Gheysen G. Transcriptome analysis of rice mature root tissue and root tips in early development by massive parallel sequencing. J Exp Bot. 2012;63:2141–57.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-∆∆CT method. Methods. 2001;25:402–8.
We thank Director, ICAR-NRRI for providing all facilities for the experiment. We thank Xcelaris Company, Ahmedabad for technical support.
The financial support for the work was provided by ICAR-National Rice Research Institute, Cuttack, Odisha, India. The funding Institute approved and allowed the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Hierarchical cluster analysis of transcripts in CSV and CTV. The median ratio (stressed/control) was log (base 2)-transformed and subjected to linkage hierarchical clustering. S1-S5 and T1-T5 denote CSV and CTV at 6, 12, 24, 48 and 24 h recovery after 48 h stress conditions, respectively. (TIF 724 kb)
Figure S2. Expression profiles of significant DEGs (up-regulated, down-regulated) and exclusively expressed genes in both genotypes i.e. CSV and CTV. Up and down-regulated genes were identified by comparing the fold change using FPKM values with control condition (S0). T denotes tolerant variety, Geetanjali (CTV) and S denotes susceptible variety, Sahabhagidhan (CSV). (JPG 44 kb)
Figure S3. GO analysis of DEGs in early, late and 24 h recovery stages of susceptible cultivar, ‘Sahabhagidhan’. (JPG 62 kb)
Figure S4. Venn diagram of DEGs under 6 h, 48 h and after 24 h recovery. A) Down-regulated genes of S1 and T1 B) Up-regulated genes of S1 and T1 C) Down-regulated genes of S4 and S5, D) Up-regulated genes of S4 and S5, E) Commonly down-regulated expressed genes of S4-T4 and S5-T5, F) Commonly up-regulated expressed genes of S4-T4 and S5-T5, G) Down-regulated genes of T4 and T5, H) Up-regulated genes of T4 and T5. S, T, 1, 4 and 5 denote CSV, CTV, 6 h, 48 h stress condition and 24 h recovery, respectively. (JPG 33 kb)
Figure S5. KEGG pathway analysis of DEGs of CSV genotype at 6 h early response phase (S1), 48 h late response phase (S4) and 24 h recovery (S5). (JPG 61 kb)
Figure S6. Differential expression of transcription factors genes involved during 48 h late response cold stress and 24 h recovery conditions of CSV genotype. a) up-regulated TF genes in S4, b) down-regulated TF genes in S4, c) up-regulated TF genes in S5, d) down-regulated TF genes in S5. (JPG 82 kb)
Table S1. Summary of RNA-Seq reads and their mapping on the rice genome. (DOCX 22 kb)
Table S2. Differentially expressed genes (DEGs) of CSV and CTV genotype at each time interval. (DOCX 13 kb)
Table S3. Comparison of DEGs of CSV and CTV genotypes at each time interval. (DOCX 14 kb)
Table S4. Differentially expressed genes (DEGs) of S1, S2, S3, S4 and S5 cold susceptible cultivar (CSV) against control sample (S0). (XLS 4297 kb)
Table S5. Differentially expressed genes (DEGs) of T1, T2, T3, T4 and T5 cold tolerant cultivar (CTV) against control sample (T0). (XLS 3788 kb)
Table S6. Top DEGs detected at different response phases of chilling stress in both the constrasting parents. (XLS 834 kb)
Table S7. Significant GO terms enriched DEGs of cold susceptible variety (CSV) on the basis of FDR corrected p-value of three different functional categories. (DOCX 23 kb)
Table S8. Significant GO terms of early response phase (S1) of CSV genotype. (DOCX 17 kb)
Table S9. Significant GO terms of late response phase (S2-S4) of CSV genotype. (DOCX 13 kb)
Table S10. Significant GO terms of 24 h recovery condition (S5) of CSV genotype. (DOCX 13 kb)
Table S11. Significant GO terms of early response phase (T1) of CTV genotype. (DOCX 14 kb)
Table S12. Significant GO terms of late response phase (T2-T4) of CTV genotype. (DOCX 13 kb)
Table S13. Significant GO terms of 24 h recovery condition (T5) of CTV genotype. (DOCX 13 kb)
Table S14. KEGG Pathway analysis of DEGs at 48 h late response phase and after 24 h recovery of CTV and their comparison with CSV genotype. (XLS 528 kb)
Table S15. KEGG Pathway analysis of DEGs at 48 h late response phase and after 24 h recovery of CSV and their comparison with CTV genotype. (XLS 489 kb)
Table S16. Distribution of differentially expressed transcription factors genes in 48 h late response phase and 24 h recovery condition in CSV and CTV genotypes. (XLS 35 kb)
Table S17. Differentially expressed TF involved at 48 h late response phase and 24 h recovery conditions of both CSV and CTV genotypes. (XLS 32 kb)
Table S18. Differentially expressed TF involved at 48 h late response phase of both CSV and CTV genotypes with annotation. (XLS 88 kb)
Table S19. Differentially expressed TF involved in after 24 h recovery condition of both CSV and CTV genotypes with annotation. (XLS 66 kb)
Table S20. List of primers used for qRT-PCR study. (XLSX 12 kb)
About this article
Cite this article
Pradhan, S.K., Pandit, E., Nayak, D.K. et al. Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis. BMC Plant Biol 19, 352 (2019). https://doi.org/10.1186/s12870-019-1922-8