Cold induced differential gene expression analysis
This work was aimed at studying the early cold stress response in rice (Oryza sativa L. ssp. indica). For this, 14 days old rice seedlings were treated at 4°C for 2h. Libraries under the control condition are denoted as CT replicates and 2h cold shock (4ºC) are indicated as CS replicates. 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 were assessed with FastQC, which revealed that all reads were of fairly good quality and without adapters. After removal of adaptors, the clean read length was approximately 63 nucleotides on an average. 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 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). As indicated by the volcano plot [Figure 1(B)], 516 genes were differentially regulated. 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 and classification of cold induced differentially expressed genes
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)]
The GO terms show that beside cellular biosynthetic processes and transcription regulation, processes related to water stress, oxidation-reduction process and lignin metabolism were significantly enriched during cold shock treatment [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 4]
The Interpro domain search data indicated that DNA binding domain and cytochromes were most abundant in the upregulated genes. [Additional file 6(B)]
Expression profiling of early cold responsive genes in IR64 variety
Transcriptome studies in rice under cold stress have identified cold-responsive genes in a number of indica rice lines, and cultivars such as 93-11, and Pusa Basmati [37-40]. Upon comparison, we have identified 72 cold induced genes differentially regulated for both our data set and 93-11 studies [Additional file 7]. We further refined this comparison by including early cold stress studies in japonica rice, Jumli Marshi [41]. Considering all the data sets, we categorized our DEGs into three groups: (A) genes common to all studies; (B) genes unique to IR64 variety, and (C) genes showing opposite regulation in IR64 under 2h shock compared to later stress time points. Category A consists of 14 genes which include cold stress responsive transcription factors: OsDREB1B, ZFP182, OsDREB1A, ERD15; players of calcium signaling pathways: OsCML16, OsACA7; Phospholipase A2; Peptidase S54; and senescence-associated protein [Additional file 2].
The 194 unique genes that were regulated only in IR64 during cold shock condition (category B) included 25 stress responsive genes, 8 transcription factors coding genes, and 20 genes involved in redox signaling. These unique upregulated 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. The 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. Among the novel transcription factors specific for IR64 data set, OsWRKY27 (Os01g0586800) was upregulated, whereas, growth-regulating factors OsGIF3 (Os03g0733600) and OsGRF7 (Os12g0484900), which act as positive regulators of cell proliferation, were downregulated under cold shock condition. 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. Category (C) genes were identified mainly by comparing our data with 93-11 indica variety subjected to 48h and 72h cold treatment. The analysis showed 28 genes that were upregulated in IR64 during 2h cold shock was downregulated at later time points in 93-11 variety, whereas 17 genes showed the opposite expression profile. The early upregulated genes include metabolic genes, transporters, protein phosphatase, and early transcription activators.
Interestingly our study shows that the DREB1 group of genes viz, OsDREB1A, OsDREB1B, OsDREB1C, OsDREB1G, OsDREB1E, OsDREB1H were significantly upregulated under cold shock conditions. Studies indicate that DREB/CBF dependent regulation is considered as the major pathway in cold acclimation and is highly conserved in various plant species [42]. 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 3]. Upregulated DREB regulons (107 genes), obtained from this analysis include abiotic stress responsive TFs, calcium-binding proteins, RLCKs, redox- signaling molecules, and other hydrophilic water-deficit responsive proteins.
Gene regulatory network induced during early cold stress
Sequence analysis suggests that around 5-7% of coding sequences in plant genomes constitute transcription factors [43, 44]. In plants, the role of AP2/EREBP, bZIP, NAC/NAM (ATAF and CUC), MYC/MYB, and WRKY transcription factor families has been elucidated in abiotic stress response that regulates stress-responsive gene expression via ABA-dependent or independent pathways [23, 45, 46]. Our data indicate the presence of 38 upregulated TFs (10% of total upregulated DEGs) and 9 downregulated TFs in cold shock transcriptome. [Figure 3(A)] The upregulated transcription factors constitute members of various families, the majority being the AP2/ERF, DREB and MYB group of TFs. Downregulated genes included nine transcription factors needed for growth and development of the plant, which consisted primarily of bHLH TFs, and growth-related TFs (Os03g0733600, Os12g0484900, Os01g0646300) [Additional file 5].
Gene regulatory network analysis of these upregulated transcription factors using STRING shows the presence of a major network consisting of 15 nodes and a second network with 3 nodes. Further, the search suggested MYB2, MYB4, DREB1B, ZFP37 and, DREB1G were highly connected and formed the central cluster. DREB1E TF is also linked with this main cluster, and it works in tandem to regulate gene expression in response to cold stress. NAC39, which connects other TF like DERF5 (ERF103), NAC077 and ERF71 also played an important role to hold these proteins together with the central cluster. Other TFs like ZFP182, ERF24, ERF26 were connected via DREB1B. 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 3(B)].
Our results show that a large number of DREB1 subfamily of TFs were significantly upregulated under cold shock, which is consistent with other cold stress transcriptomic studies [37, 38]. This indicates that the genes which regulate through DRE binding site are triggered early during cold shock, to induce downstream regulators for mounting the entire cold stress response in rice plants. The large number of other AP2/ERF transcription factors genes, such as OsERF141, OsDERF5, EREBP139, OsERF102, and OsDERF8 were also upregulated for early cold stress response in indica rice [47]. The upregulation of the cold specific MYB factor genes, OsMYB2 and OsMYB4 indicate their important role in early signaling events. ZFP182 (Os03g0820300) gene coding for a TFIIIA-type zinc finger protein type transcription factor known to be involved in multiple abiotic stress tolerance mechanisms in rice [48, 49] was significantly upregulated under cold 2h shock in this study.
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 a number of redox homeostasis proteins were induced under cold shock conditions [Figure 4(A); Additional file 5] 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 are mostly belonging to osmoprotectants activated due to dehydration stress [Figure 4(C)].
Expression of cold responsive genes in different unexplored indica cultivars
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.
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 profile of ten genes (six upregulated and four downregulated, as validated in IR64) from the above-validated list were tested for this study [Figure: 5(B)]. Our data shows the expression of OsDREB1b, ONAC039, OsCML31, ERD15, genes were upregulated in all ten varieties. However, upregulation of OsLEA14 and OsCIPK7 genes were restricted to some of the varieties. In the case of the downregulated genes, OsGRF7, OsPRP1, OsbHLH155 showed consistent downregulation in all varieties. Our data also indicated that expression of some of the cold-responsive genes was significantly upregulated in the Heera and GB1 cultivars under control conditions, compared to the other varieties. This higher level of expression of cold-responsive genes under control condition may be related to their better adaptability to cold shock.