Cold shock induces differential gene expression in IR64 indica rice
This work aimed at studying the early cold stress response in rice (Oryza sativa L. ssp. indica) involved the treatment of 14 days old rice seedlings to cold shock (at 4°C) for 2h. cDNA libraries generated from IR64 seedlings grown under the control (28°C±1°C) and 2h cold shock (4ºC) conditions are denoted as CT replicates and CS replicates respectively. A total of 58.53 million (CT1), 33.61 million (CT2), 60.85 million (CS1), 39.82 million reads (CS2) were generated for the IR64 rice cultivars. Statistics of cleaned reads (approx. 63 nucleotides) were assessed with FastQC, which revealed that all reads were of fairly good quality and without adapters. The read mapping was carried out using HiSat2, resulting in 95.5% reads mapping to reference Japonica genome Oryza sativa japonica (Os-Nipponbare-Reference-IRGSP-1.0 (IRGSP-1.0) [Table:1]. Normalized expression profiling was done on the aligned reads resulting in the identification of 32161 transcripts expressed in at least one of the four samples profiled. As shown in Figure 1(A), 24988 transcripts were found in both the control data sets, whereas 24836 transcripts were common for both cold shock replicates. When compared between control and cold data sets, 72% (23232) of the transcripts were expressed in both, indicating a basal level of expression. Further, 539 transcripts were found exclusively expressed in control, whereas 931 transcripts were expressed only under cold shock condition.
We used Deseq2 package for differential expression analysis of the genes. A filter with a p-value cut-off of <0.05 and log2fold change ≥ 1.5 and ≤ -1.5 was set as the criteria to identify the differentially expressed genes (DEGs). These DEGs were visualized using volcano plot [Figure 1(B)] to understand the distribution of up and downregulated genes. For this analysis, FPKM of >=0.1 for a transcript was considered as expressed supported by a median read count of at least five reads per transcript covering 100% of the sequence. Among these DEGs, 380 genes were upregulated, whereas, 136 genes were downregulated in cold shock (CS) Vs control condition (CT) seedlings. [Additional file 1] Analysis of unsupervised hierarchical clustering of differentially expressed transcripts shows distinct gene expression patterns of up and down-regulation levels during cold shock treated (CS) [Figure 1(C)]. The differentially expressed genes could be categorized into five different clusters, based on their expression patterns. The largest group, Cluster V comprises of the differentially expressed genes, upregulated under cold shock, whereas, Clusters II and IV constitute the downregulated genes under cold shock. The clustering of all these replicates exhibited high sample reproducibility.
Functional annotation of cold shock-induced genes indicate a significant increase in cold-responsive TF and ROS activity
To understand the biological function of the cold induced differentially expressed genes (DEGs), GO enrichment was performed using an FDR adjusted p-value of ≤ 0.05 as the cut-off. The Blast2GO analysis for 516 DEGs featured 234 GO term annotation in biological process, 273 in molecular function (MF), and 262 for cellular component. Comparative analysis of the upregulated and downregulated GO terms indicates cell division (GO:0051301), proliferation (GO:0008283), developmental processes (GO:0032502) and growth (GO:0040007) were specific for downregulated genes. GO terms such as transport (GO:0006810), and homeostatic process (GO:0042592) were majorly associated with the upregulated genes [Figure 2(A)].
GO analysis identified oxidation-reduction process, processes related to water stress and lignin metabolism were significantly enriched during cold shock treatment, in addition to generic terms such as cellular biosynthetic processes and transcription regulation [Figure 2(B)]. Under stress response, significant enrichment for GO-terms such as response to alcohol (GO:0097305), response to temperature stimulus (GO:0009266), response to abscisic acid (GO:0009737), response to lipid (GO:0033993), response to acid chemical (GO:0001101), response to osmotic stress (GO:0006970), and cold acclimation (GO:0009631) was observed [Figure 2(C)]. For the oxidation-reduction process, response to Oxygen-containing compound (GO: 1901700) and lignin catabolic process (GO:0046274) were significantly enriched under cold shock. GO-molecular function (MF) terms comparison indicates that metal ion binding (GO:0005488), oxidoreductase activity (GO:0016491) and transcription regulator activity (GO:0140100) were enriched for upregulated genes [Figure 2(C)]. GO enrichment analysis was also performed for the downregulated genes, but no significantly enriched terms were detected.
The biological pathway associated (KEGG pathway) with cold shock response was analyzed using BLASTKOALA (24.9% of input sequences). The analysis revealed that most genes were assigned to metabolism (40) of carbohydrate, amino acid, lipid and secondary metabolites; environmental information processing (15), genetic information processing (5) like transcription, translation and protein processes [Additional file 6(C)]. Further, KEGG-BRITE reconstruction revealed that compared to control, a higher number of genes were assigned to ko01000 Enzymes (42), ko02000 Transporters (8), ko01003 Glycosyltransferases (6), ko03000 Transcription factors (4) and, ko04147 exosomes (4) in the cold shock treated sample [Additional file 2]. The Interpro domain search data indicated that DNA binding domain and cytochromes were most abundant in the upregulated genes [Additional file 6(B)].
Gene regulatory network induced during early cold stress
Analysis of the DEGs showed the presence of a milieu of stress responsive genes upregulated, which included heat shock protein genes (Os02g0758000, Os03g0266900, Os06g0253100) Terpene synthases (OsTPS1, OsTPS31), Dehydrins, LEA and RAB group of proteins coding genes (OsDHN1, OsLEA28, OsRAB16, OsLEA14), pathogen-related proteins and chitinase and glucanases. Among signaling molecules, calcium-calmodulin molecules and receptor kinases and protein phosphatases, along with several redox homeostasis proteins were induced under cold shock conditions. Besides hydrophilic proteins and signaling protein, our data set indicate the presence of 38 upregulated TFs (10% of total upregulated DEGs) and 9 downregulated TFs in cold shock transcriptome [Figure 3(A) and 3(B); Additional file 3]. Gene network of these upregulated DEGs shows three major clusters that are highly interconnected. Cluster I comprise Zinc finger and NAC transcription factors and signaling proteins such as calmodulin and kinases. Cluster II represents DREB/AP2 and MYB transcription factors as major nodes. Cluster III contain proteins that mostly belong to osmoprotectants activated in response to dehydration stress [Figure 3(C)].
Cell wall modification and ROS generation are crucial to stress perception during early cold shock in IR64
Differentially expressed gene set unique to this study [Additional file 4], has a significant number of genes responsible for cell wall modification and ROS generation [Figure 3(A) and 3(B)]. The genes coding for redox molecules comprised majorly of lignin catabolic laccase genes (OsLAC10: Os02g0749700, OsLAC17: Os10g0346300, OsLAC23: Os11g0641500, and OsLAC29: Os12g0258700), the germin-like oxalate oxidases, and other ROS generating enzymes. These findings suggest that the generation of ROS occurs during early cold shock and is essential for activating the redox signaling at the later stages of the stress response. Other stress responsive genes, unique to 2h cold shock include the terpene biosynthesis genes (OsTPS1, OsTPS10), salt stress responsive lectin proteins (Os01g0348800, Os01g0348900) and receptor-like kinases (Os11g0672200, Os04g0540900). Other genes that were induced within 2h include cell wall degrading enzymes, such as, chitinases (Os05g0399300, Os11g0701200, Os11g0702100), cellulases (Os01g0946600, Os01g0946700), pectin methylesterase (Os04g0458900), and a group of membrane transporter genes.
Cold shock-induced TFs upregulated under cold shock constitute major gene regulatory networks
Sequence analysis suggests that around 5-7% of coding sequences in plant genomes constitute transcription factors [37, 38]. In plants, the role of AP2/EREBP, bZIP, NAC/NAM (ATAF and CUC), MYC/MYB, and WRKY transcription factor families has been elucidated in the abiotic stress response that regulates stress-responsive gene expression via ABA-dependent or independent pathways [23, 39, 40]. The upregulated transcription factors in this study belong to various families; the majority being the AP2/ERF, DREB and MYB group of TFs [Figure 4(A)]. Downregulated genes included nine transcription factors needed for growth and development of the plant, which consisted primarily of bHLH TFs, and growth-related TFs such as OsGIF3 (Os03g0733600), OsGRAS1(Os01g0646300) and OsGRF7 (Os12g0484900) [Additional file 3]. A gene regulatory network analysis using STRING shows that these upregulated transcription factors constitute a major network consisting of 15 nodes and a second network with 3 nodes. Further, the search suggested MYB2, MYB4, DREB1B, ZFP37, DREB1E and, DREB1G were highly connected and formed the central cluster [Cluster I, Figure 4(B)]. The cluster II consists of NAC39, which connects other TF like DERF5 (ERF103), NAC077 and ERF71. The heat shock transcription factors, HSF21, HSF11, are connected among themselves and also connected to central cluster via DREB1E, DREB1G and MYB2. HOX transcription factors constitute a secondary network which may contribute to cold stress response [Figure 4(B)].
Differential expression of DREB1 regulon genes is integral to early cold stress response
The RNA-seq results from cold shock and control condition were validated using quantitative real-time PCR (qRT-PCR) for 25 genes. This included 15 genes from upregulated DEGs and 10 genes from downregulated DEGs of the cold shock transcriptome. Figure 5(A) shows the fold change obtained from three biological replicates for these 25 genes and their corresponding log2fold change obtained from RNA-seq data. As seen in the figure, 80% of validated genes exhibited differential expression with significant p-value (p< 0.05) that matches with the genes’ expression profiles from RNA-seq data.
In this study, we observed that the DREB1 group of genes viz, OsDREB1A, OsDREB1B, OsDREB1C, OsDREB1G, OsDREB1E, OsDREB1H were significantly upregulated under cold shock conditions. This indicates that the transcription factors which recognize the DRE motif are triggered early during cold shock, to induce downstream regulators for mounting the entire cold stress response in rice plants. Previous studies have indicated that DREB/CBF dependent regulation is considered as the major pathway in cold acclimation and is highly conserved in various plant species [41]. Among the 516 DEGs reported in this study, 27.7% of genes were identified to have at least one DRE-binding motif at the upstream region [Additional file 5]. The upregulated DREB1 regulon genes identified in this study may be grouped into four major categories: the hydrophilic proteins, LEA, DHN1 (COR410) proteins, stress associated proteins (SAP17), osmotin precursors, photoperiod sensitive and transporter proteins; signal transducing molecules, such as phosphatases, membrane kinases, Ca2+-CaM proteins, and various cold-responsive ubiquitin ligases; multiple classes of enzymes including catabolic chitinases, exocyst complex proteins and oxidoreductases, laccases; and the final group of other zinc fingers, HOX, and AP2/ERF group of transcription factors. 20 genes from the 107 upregulated DREB regulon genes were validated, among which expression of the genes coding for AP2/ERF transcription factor OsERF102 (Os09g0457900), Calcium-binding protein OsCCD1 (Os06g0683400), and the exocyst subunit EXO70 family protein, OsEXO70FX14 (Os01g0905300) exhibited high levels of significant upregulation [Figure 5(B)].
Expression of cold-responsive genes in different unexplored indica cultivars
We extended our study to ten different indica rice cultivars that were never characterized for cold stress response. These varieties are mostly high yielding varieties (HYVs), including both hybrids and field-selected varieties, cultivated all over India [Table: 2]. The expression profiles of 18 genes (14 upregulated and 4 downregulated, as validated in IR64) from the above-validated list were tested for this study. Our data shows upregulation of transcript levels of OsDREB1b, ONAC039, OsCML31, ERD15, OsLEA14, OsCIPK7 and the DREB regulon genes Quinone Oxidoreductase, OsFbox208, Transcription factor IIS, Pyridoxal phosphate dependent transferase, OsERF102, OsCCD1, Pectin-glucuronyltransferase, and OsEXO70FX14 genes post cold shock treatment in the majority of the varieties [Figure: 6(A)]. In the case of the downregulated genes, OsGRF7, OsPRP1, OsbHLH155 showed consistent down-regulation in all varieties in response to cold shock treatment [Figure: 6(B)]. Our data indicated that expression of majority of the cold-responsive genes was significantly upregulated in the CB1, and Heera cultivars under control conditions, compared to the other varieties. Interestingly, the DREB regulon genes OsEXO70FX14 (Os01g0905300), Transcription Factor IIS (Os06g0693800), Pyridoxal-phosphate dependent transferase (Os04g0614500), OsFbox208 (Os04g0414500) which exhibited highly significant upregulation under cold shock conditions in most varieties, was highly upregulated under control condition in CB1 rice variety [Figure: 6(C)]. This higher level of expression of cold-responsive genes under control condition in CB1 and Heera varieties prompted us to examine the physiological response of the indica rice varieties used in this study.
The germination experiments performed in this study revealed that the indica rice variety CB1 remains nearly unaffected by cold stress (72hrs) prior to germination, whereas, the variety IR36 shows the highest sensitivity to the same. As illustrated by Figure 6(D), the decrease in percent germination of CB1 with increasing cold stress is not significant, compared to the other varieties. The rice variety Heera also exhibits a slight reduction in the percent germination, which indicates that this variety may also have better resilience to the low-temperature stress. IR64 and IR36 varieties show a significant drop in their percent germination, which suggests their high level of sensitivity to cold stress conditions. All the other indica varieties used in this study exhibit low to moderate sensitivity to cold stress as shown by their percent germination [Figure 6(D)].
The electrolyte leakage of the indica rice varieties IR64, IR36, Hamsahamas, CB1, and Heera under control and cold stress conditions (4°C treatment for 2h, 4h, 24h, 48h, and 72h) was studied by determining the relative electrical conductivity (REC%) of leaf tissue samples of 14 days old seedlings. As shown in Figure 6(E), REC% shows a significant increase under cold stress in IR64, IR36, and Hamsahamas varieties. In contrast, a delayed increase of REC% was observed for CB1, and Heera varieties (48h and 72h cold stress treatments). There was no significant change in REC% in CB1 and Heera during initial cold stress time points.
The ROS generation of the five indica varieties, under control and cold stress conditions (4°C treatment for 2h, 4h, 8h, and 24h) was quantitated using H2DCFDA. Our data shows that the ROS generation is induced significantly under cold shock conditions in IR64, IR36 and Hamsahamas varieties, with higher levels of increase in IR64 and IR36 varieties. As exhibited in Figure 6(F), the ROS levels in these three varieties show a significant increase with increasing stress duration, up to 8 hours. However, the ROS generation decreases after 24h of cold treatment in IR64, IR36 and Hamsahamas, suggesting the induction of the ROS scavengers during later stages of cold stress. The tolerant varieties CB1, and Heera show no significant changes in the ROS levels during cold stress.