The phenotypic changes of rice seedlings at different stress stages are shown in Fig. 1. After 1 day of cold treatment, the appearance of 9311 and DC907 did not change significantly, and their growth conditions were similar, indicating that 9311 and DC907 had the ability to withstand short-term low temperatures. After 3 days of cold treatment, the leaves of 9311 curled and wilted, whereas those of DC907 did not. After 5 days of cold treatment, the leaves of 9311 began to wilt and turn yellow, whereas those of DC907 had a lower wilting degree. These results showed that DC907 had better low-temperature resistance than 9311 under longer low-temperature conditions. In the process of recovery test, with prolonged recovery time, CK and T1 treatment grew normally. In the recovery process of T3 treatments, 9311 leaves gradually turned yellow with prolonged recovery time, and all plants died at the end of 5 days of recovery. Conversely, for DC907, most original curled leaves gradually recovered and only a few leaves finally turned yellow and died with prolonged recovery time. In the recovery process of T5 treatment, leaf curl appeared in most rice varieties. With prolonged recovery time, all 9311 plants gradually turned yellow and eventually died at 5 days of recovery, whereas DC907 also turned yellow and died at the same time. After 5 days of recovery, only a few plants could maintain normal growth. Moreover, the cold resistance of DC907 was stronger than that of 9311, and the recovery ability was stronger after cold injury.
Global transcriptome analysis of rice varieties
To avoid increasing the difficulty of data analysis and the analysis error caused by comparing multiple similar groups, some treatment groups were selected for transcriptome sequencing and subsequent analysis based on the above phenotypic observations. The materials and methods for transcriptome sequencing are described in detail. The transcriptomes of 10 experimental treatments (AT1, AT1H1, AT3, AT3H1, BT1, BT1H1, BT1H3, BT3, BT3H1, and BT3H3) and 2 blank controls (ACK and BCK) were sequenced.
RNA-seq analysis was performed using rice aboveground-tissue samples. Each treatment contained three biological replicates of 36 samples. The number of raw reads ranged from 39.72 million to 47.26 million, and the number of clean reads obtained by data filtering ranged from 38.96 million to 46.05 million. The percentage of Q30 base to total base ranged from 91.59% to 95.24% (Additional file 1). Rice reference genomes were compared using HISAT2 software [33], and results showed that the relative efficiency of 9311 and DC907 was 82.96%–95.15% and 89.89%–95.15% respectively. Amongst them, the relative efficiency of AT1_3 with genome was low (67.51%) (Additional file 2). However, we consider it as sequencing based on the previous quality-control results and sample correlation analysis. Considering the lack of depth, we did not believe it would adversely affect the follow-up analysis and decided to use it for such analysis.
A total of 43 255 genes were identified by Cufflinks and Cuffcompare, including 40 745 annotated genes and 2509 new genes. Overall, about 51.6% of the genes were expressed in at least one of the 12 samples (fragments per kilobase of transcript per million [FPKM] > 1). The number of expressed genes ranged from 46.7% (T1) to 51.0% (CK) in 9311 and from 48.4% (T3H3) to 51.9% (T1H1) in DC907 (Fig. 2a). In all treatments, the number of genes with low (1 < FPKM < 3), medium (3 < FPKM < 15), high (15 < FPKM < 60), and very high expression (FPKM < 60) were similar (Fig. 2b). The number of genes with low expression (CK and T1) had the highest gene expression in the two cultivars.
Identification of differentially expressed genes (DEGs) in different treatments
To investigate the transcriptional differences between two rice varieties at different treatment stages, we identified the DEGs of the two rice varieties at each treatment stage. Each treatment stage was compared with CK, and DEGs were filtered with padj < 0.5. Based on this formula, 4178 and 6054 of genes were significantly differentially expressed in 9311 and DC907. Amongst 9311 cold-sensitive rice cultivars, 2911 DEGs were identified in AT1, with the largest number of DEGs; 174 DEGs were identified in AT1H1, with the least number of DEGs. In the cold-insensitive rice cultivar DC907, 4380 DEGs were identified in BT3, with the largest number of DEGs; 285 DEGs were identified in BT1H3, with the least number of DEGs.
Interestingly, the proportion of upregulated genes in AT1 was 70% and that of upregulated genes in AT3 was 51%. The proportion of upregulated genes in all differential genes decreased with increased cold treatment time; the proportion of upregulated genes in BT1 was 74%, and that of upregulated genes in BT3 was 60%. The proportion of upregulated genes in DC907 decreased with increased cold treatment time. The proportion of upregulated genes in DC907 was slightly higher than 9311 on the first day of cold treatment, but the difference enlarged on the third day of cold treatment. The changes in number of these abnormal genes indicated that they may have had different molecular stress mechanisms at different stages of cold treatment. We also found that the number of DEGs in recovery treatment after 3 days of cold treatment was higher than that in recovery treatment after 1 day of cold treatment (Table 1). Results showed that with increased effect of cold treatment on rice, more genes were needed to participate in the recovery process to restore normal rice growth.
Table 1. Numbers of DEGs in each experimental-processing stage.
Sample
|
Up
|
Down
|
Total
|
AT1 vs ACK
|
2034
|
877
|
2911
|
AT1H1 vs ACK
|
143
|
31
|
174
|
AT3 vs ACK
|
576
|
528
|
1104
|
AT3H1 vs ACK
|
522
|
501
|
1023
|
BT1 vs BCK
|
1612
|
558
|
2170
|
BT1H1 vs BCK
|
73
|
159
|
232
|
BT1H3 vs BCK
|
182
|
103
|
285
|
BT3 vs BCK
|
2658
|
1722
|
4380
|
BT3H1 vs BCK
|
292
|
300
|
592
|
BT3H3 vs BCK
|
345
|
608
|
953
|
Functional-enrichment analysis of DEGs
A total of 7 506 DEGs were screened out after subcooling stress and recovery experiments. To understand the regulatory effects of these genes in different experimental approaches, Gene Ontology (GO) enrichment analysis was performed on DEGs in different experimental groups. At the T1 stage, DEGs were categorised in different GO terms, including RNA biosynthetic process, cellular metabolic process, and primary metabolic process. Results indicated that the transcriptional, translation, and metabolic activities of rice increased significantly within 1 day of cold treatment. At the T3 stage, the GO terms related to photosynthesis were significantly enriched. This finding showed that after a long period of cold treatment, the photosynthesis system of plants was damaged. In the process of recovery, the significantly enriched GO terms showed difference between the two rice varieties. In 9311 cultivars, many GO terms related to antioxidant activity were significantly enriched at the T1H1 and T3H1 stages, but a large number of genes at T3H1 were depressed. In DC907 cultivars, a significant difference existed in GO enrichment between recovery 1 and 3 days. At the T1H1 and T3H1 stages, many specific genes were enriched in flavonoid metabolic process and flavonoid biosynthetic process GO terms. At the T3H1 stage, GO terms related to ROS scavenging were enriched, such as detoxification, response to toxic substance, and oxidation-reduction process (Fig. 3). At different treatment stages, rice may resist low-temperature stress through various regulatory methods.
Identification of DEGs at different treatment stages of rice seedlings
We analysed the expression profiles of genes related to plant-hormone biosynthesis, plant-hormone signal transduction, ROS scavenging, and plant antifreeze biosynthesis. CDPKs are a class of protein kinases that exist only in plants and some protozoa. Calcium participates in the signal-transduction process of calcium ion as a second messenger [34] and plays an important role in regulating plant growth and development signals and stress signals [35, 36]. Six CDPK-coding genes were upregulated during cold stress (Fig. 4a). Notably, the overall expression of CDPK-coding genes in 1 day of cold treatment was higher than that in 3 days of cold treatment. CML family proteins containing predicted Ca2+-binding EF-hand motif are potential calcium sensors [37]. Under abiotic stresses such as high temperature, low temperature, hypoxia, and darkness, the expression of CML protein significantly increases [38]. A total of 12 CML coding genes were identified to be upregulated under cold stress. The total expression level in 9311 was higher than that in DC907. Calcium participates in the signal-transduction pathway of rice cold-stress response, which is widespread in sensitive and cold-insensitive rice lines. Calcium signalling pathway primarily plays a regulatory role in the early cold-stress process of rice.
Reactive oxygen species (ROS) are important intermediates in the plant abiotic-stress process. An appropriate amount of ROS can regulate abiotic-stress responses in plants through signal-transduction pathways. With prolonged stress time, ROS accumulation inflicts oxidative damage to plants. During the cold treatment and restoration of rice, numerous differentially expressed antioxidant synthesis genes were screened (Fig. 4b). Glutathione transferase is ubiquitous in plants. It plays an important role in plant stress response by preventing the oxidative damage of ROS [39, 40]. Six glutathione S-transferase coding genes were overexpressed in AT1. At the AT3 stage, three glutathione S-transferases, one glutathione peroxidase, and one monodehydroascorbate reductase increased after cold treatment. Five glutathione S-transferase coding genes were induced in AT1. Two glutathione S-transferase coding genes BGIOSGA033329 and BGIOSGA001631, one glutathione peroxidase coding gene BGIOSGA016904, and one single dehydroascorbate reductase coding gene BGIOSGA029182 were induced to be expressed in BT3. In Arabidopsis, regulating the redox state through the ascorbic acid–glutathione cycle can reduce oxidative stress under many abiotic stresses [41]. Results showed that the ASA-GSH cycle played a positive role in the regulation of cold stress in rice, and with prolonged cold-treatment time, more enzymes became involved in the induction of expression.
Sugar is an important regulator of plant growth, morphogenesis, and expression. It also plays an important role in maintaining basic plant-life activities and signal transduction. Trehalose accumulation in plants can improve plant tolerance to abiotic stress [42, 43]. In the present study, many genes involved in trehalose biosynthesis were induced and expressed. Two trehalose 6-phosphate phosphatase encoding trehalose 6-phosphate were significantly induced to be expressed in cold-treated day (AT1 and BT1) (Fig. 4c). It can remove phosphate from trehalose phosphate to produce free trehalose [42]. The expression of Trehalase gene BGIOSGA031535 increased in AT3H1 and BT3. Its main biological function is to hydrolyse trehalose into glucose. Four trehalose 6-phosphate phosphorylase coding genes were identified in BT3at, which were upregulated and differed from the above genes. These results indicated that trehalose played an important role in maintaining the normal life activities of rice. Under cold stress, different genes were expressed at different treatment stages, causing rice to have cold tolerance. The role of amylase in plants is to hydrolyse starch and ultimately decompose it into glucose. The expression of one α-amylase and two β-amylases were found to be induced during treatment. BGIOSGA028835 and BGIOSGA031385 were induced to be expressed during the first day of cold treatment, and BGIOSGA011829 was induced to be expressed in BT3. It can mediate maltose accumulation under freezing stress, thereby helping to protect the photosynthetic electron-transport chain [44, 45]. Two genes encoding ribulose bisphosphate carboxylase (ribulose bisphosphate carboxylase) inhibited during treatment were involved in the immobilisation of carbon dioxide and the oxidation of pentose substrates in plants. Results showed that the photosynthesis capacity of plants decreased during cold treatment and short-term recovery, consistent with previous experimental observations and functional-enrichment analysis.
We found that most plant-hormone biosynthesis and signal-transduction genes were upregulated during cold stress, suggesting that cold-stress response was the main process of cold tolerance in rice, and that the regulation of gene expression may play only an auxiliary role during recovery. Four NCED coding genes were significantly differentially expressed; amongst them, NCED was the rate-limiting enzyme in the process of ABA biosynthesis [46]. BGIOSGA025169 and BGIOSGA013244 were upregulated in BT3, whereas BGIOSGA083804 was induced during DC907 cold stress and in 9311 for 3 days. This finding indicated that the ABA biosynthesis pathway played an important role in the cold tolerance of rice. CYP707A7 participated in the oxidative degradation of abscisic acid. BGIOSGA029635, the coding gene of CYP707A7, was significantly induced during cold treatment. The expression of CYP707A7 in AT1 was higher than that in BT1H1 and AT1H1. We speculated that rice needed to accumulate more ABA to participate in cold signal transduction during 3 days of cold treatment. SAPK2 plays an important role in ABA signal transduction. We found that the SAPK2 coding gene BGIOSGA026164 was induced to be expressed in two cultivars during cold treatment, and that the expression level was higher than 3 days in 1 day. In rice, SAPK2 was strongly induced under salt and drought stress, and SAPK2 gene overexpression could reduce rice sensitivity to salt [47]. Jasmonic acids (JAs) are an important component of plant hormones. Lipoxygenase is a key enzyme in the biosynthesis of JA. In the current work, we found that four genes encoding lipoxygenase were induced to be expressed during cold stress. BGIOSGA037588 was highly expressed in BT3 and the other three in AT1. TIFY protein is a repressor of JA reaction [48]. Six TIFY coding genes were identified to be differentially expressed during the treatment. Notably, TIFY protein was very highly expressed in BT3 compared with the DC907 treatment group. In the 9311 treatment group, they were expressed in AT1, AT3, and AT3H1 to varying degrees. Cytokinin dehydrogenase 8 (CKX8) can catalyse cytokinin oxidation in plants. BGIOSGA016782, the coding gene of CKX8, was upregulated in BT3, and its expression was higher than that of other treatments. Flavin-containing monooxygenase may be involved in auxin biosynthesis in rice. We found that BGIOSGA015673, the coding gene of flavin-containing monooxygenase, was significantly induced to be expressed during cold treatment for 1 day. In this study, we also found that nine genes encoding auxin-response protein(ARF) were significantly differentially expressed during treatment. IAA4 was upregulated during the first day of cold stress, IAA3 was inhibited in each treatment, and other ARF genes were inhibited or induced to varying degrees during the recovery treatment of DC907. Results showed that auxin signal transduction played an important role in cold stress and recovery. Two genes encoding the ethylene-response sensor proteins ERS1 and ERS2 were induced to be expressed during cold treatment. BGIOSGA009308 and BGIOSGA009973 were significantly upregulated during cold treatment. Through the above analysis of plant-hormone biosynthesis and signal-transduction pathways, we found complex synergistic regulatory networks and signal-transduction pathways amongst ABA, JA, auxin, and ethylene, which enabled rice to withstand low temperatures.
In addition to the aforementioned pathways that had been extensively studied under abiotic stress, we screened several genes that were differentially expressed during cold treatment and recovery, and these genes may play important roles in regulating cold tolerance or recovery in rice. From literature, we can find that Phosphoinositide phospholipase C coding gene is significantly downregulated during cold treatment. It can mediate the production of the second messenger molecule diacylglycerol and inositol 1,4,5-triphosphate through the phosphatidylinositol signalling pathway [49]. Herein, we also found that one gene encoding Inositol polyphosphate multikinase (IPK2) and one gene encoding Inositol-tetrakisphosphate 1-kinase 2 (ITPK2) were upregulated during cold treatment. Notably, their expression level on 1 day after cold treatment was higher than that 3 days after cold treatment and BT3. Both genes were involved in inositol phosphate metabolism [50, 51]. CAM1-3 was induced to be expressed during cold treatment, and the highest expression was in AT1. The differential expression of these genes suggested that the activation of the phosphoinositol signalling pathway may play a potential role in maintaining cell homeostasis under cold stress. Polyamines are plant growth regulators that participate in important physiological processes, such as plant growth and development, sex differentiation, aging, and environmental adaptation [52, 53]. We found that two polyamine oxidase coding genes were significantly differentially expressed during treatment. BGIOSGA014319 was induced in BT3 and BGIOSGA014122 was significantly downregulated during the recovery of 9311 and DC907 after 3 days of cold treatment. Polyamine oxidase participates in secondary amino oxidation and plays an important role in regulating intracellular polyamine solubility [54]. After cold treatment, we found that the genes involved in spermidine biosynthesis showed an induced expression trend. Spermine may play an antifreezing role in rice. Gamma-aminobutyric acid (GABA) is a non-protein amino acid converted from l-glutamic acid. Aluminum activates malate transporter activity to regulate plant growth, development, and stress-response signals [55-57]. The expression of the glutamate decarboxylase gene BGIOSGA012193 was significantly upregulated after cold treatment and increased with prolonged cold-treatment time. Glutamic acid decarboxylase can catalyse the conversion of l-glutamic acid to GABA. According to the above experimental results, GABA and spermidine may play the role of antifreeze protector in rice and induce expression in cold treatment.
Numerous transcription factors are reportedly involved in abiotic stress [58-60]. Two transcription-factor encoding genes, namely, the dehydration-responsive element binding factors DREB1B and DREB1C, were significantly upregulated at the BT3 stage. Interestingly, DREB1B was also significantly induced in 9311 cultivars during 1 day of cold treatment (Fig. 5). The high expression of DREB1B may endow tobacco plants with resistance to high salt, cold, and drought stress, and inserting OsDREB1B gene into tobacco can make plants more tolerant to oxidative and freezing stress than the wild type [61]. In Arabidopsis, DREB1C negatively regulates the expression of CBF1/DREB1B and CBF3/DREB1A and plays a negative role in stress tolerance [62]. However, in rice cultivars, DREB1C may play an essential role in cold tolerance. In DC907 cultivars, stress-induced transcription factor NAC1 is induced after cold stress. It is also upregulated in AT1. In rice, OsPP18 is regulated by SNAC1 and regulates ROS homeostasis through the abscisic acid-independent pathway, thereby conferring oxygen stress and drought tolerance [63]. Ten putative WRKY transcription factors were filtered from all upregulated DEGs. WRKYs had different biological functions in plants and participate in response to biotic and abiotic stresses [64-66].
Analysis of gene-coexpression network
We used WGCNA analysis to further understand the gene expression of the two near-isogenic lines rice in different treatments and to screen out the characteristic genes related to cold tolerance. After filtering raw data in materials and methods, 19 894 genes were retained for WGCNA analysis. The coexpression network was constructed based on the gene expression of all samples. The module was a cluster of highly interconnected genes, agreeing that the genes in the cluster were highly correlated with one another. A total of 26 modules were identified for different colours, and the grey for genes meant no assignment to any module. Considering the specific expression of genes in different samples, the correlation coefficients amongst different samples of the 26 modules differed. Notably, 13 out of the 26 coexpression modules showed high correlation with single sample (r > 0.6). WGCNA analysis can also be used to construct gene-expression networks, where each node represents a gene and the line (edges) represents the correlation between two connecting genes (Fig. 6). Hub genes refer to genes with high connectivity in modules. Hub genes with high connectivity are usually regulatory factors (the upstream position of the regulatory network), whereas those with low connectivity are usually downstream genes in the regulatory network (e.g., transporters, catalytic enzymes, etc.). Hub genes are linked to more genes in the interaction network and have a high KME value.
In the yellow module, 1410 genes highly correlated with AT1 were identified. A total of 39 specific genes of AT1H1 were found in the orange module. A total of 1097, 1095, and 1114 genes were found to be highly correlated with AT3 in red, black, and green. A total of 1455 AT3H1-specific genes were identified in the brown module. In the dark-grey and cyan modules, 62 and 210 genes were found to be highly correlated with AT1. A total of 114 specific genes of BT1H3 were identified in the light-yellow module. In the dark-orange and magenta modules, 35 and 888 genes highly associated with BT3 were found. In the midnight-blue module, 206 genes were found to be highly correlated with BT3H3.
In the yellow module, 80 gene-coding transcription factors were identified, including 5 bHLH transcription factors, 3 bZIP transcription factors, 4 HOX transcription factors, 9 AP2-EREBP transcription factors, 2 MADS transcription factors, 5 MYB transcription factors, 8 NAC transcription factors, 1 TIFY protein, and 6 WRKY transcription factors. NAC family transcription factors (BGIOSGA001951) were identified as candidate Hub genes for this module. In the cyan module, 13 presumptive transcription-factor genes were identified, including 1 ARF, 1 bZIP, 1 bHLH, 1 NAC, and 2 Tubby-like F-box domain transcription factors. Calcium-transporting ATPase and auxin-response protein were identified as candidate central genes for this module.
As shown in the figure, the dark-orange and magenta modules had relatively high correlation, and both of them were presumably BT3-specific genes. In the dark-orange module, a gene encoding C3H transfer factor (BGIOSGA004750) was identified. Peptidylprolyl isomerase was identified as a candidate central gene for this module. In the magenta module, four transcription factors of AP2 domain, two bHLH, two bZIP, one MYB, two NACs, and three WRKY transcription factors were found in the magenta module. BGIOSGA002846 and membrane proteins associated with secretory vectors were identified as candidate central genes for this module.
Midnight-blue and turquoise modules were highly correlated, and both modules were positively correlated with BT3H3, so the genes in both modules may be BT3H3-specific genes. In the midnight-blue module, two presumptive transcription factors, one AP2 family transcription factor, and one bHLH transcription factor were found. BGIOSGA015677 is an ubiquitin-transferase coding gene identified as a candidate central gene for this module. In the turquoise module, 63 transcription factors were identified, including 5 bHLH, 5 bZIP, 9 MYB, 3 CH3, 5 HB, 9 NAC, and 4 WRKY transcription factors. Before screening, Hub genes were not annotated, and BGIOSGA021307 and BGIOSGA021806 genes were relatively high in expression. We speculated that these two genes were candidate central genes of the turquoise module.
In the cluster diagram, black and green are the specific module genes of AT3. In the black module, we identified 26 presumptive transcription-factor proteins, including 2 bHLH, 3 bZIP, 2 MYB, 5 NAC, and 3 WRKY transcription factors. A LOB transcription-factor coding gene (BGIOSGA010454) and phosphatidyl cytidyltransferase were identified as candidate central genes for this module. In the green module, 60S ribosomal proteins l44, L18a, L11, and SPT4 homologues were identified as candidate central genes for this module. Brown was an AT3H1 characteristic module gene. In this module, we identified 43 presumed transcription-factor proteins, including three AP2 domain transcription factors, one bHLH, five bZIP, two HBs, seven MYB, four NAC, three TIFY proteins, six WRKY, and one leucine zipper domain homologue protein. Aspartate aminotransferase, Peroxygenase, and Reticulon-like protein were identified as candidate central genes for this module.
Dark-grey was a characteristic module gene of BT1. In this module, we identified four presumed transcription-factor proteins, including one C2H2, one G2-like, one AP2 domain protein, and one HB transcription factor. DelLA protein SLR1 and xylan transglucanase/hydrolase were identified as candidate central genes for this module. Light yellow was a characteristic module gene of BT1H3, in which two predictive coding transcription-factor genes zf-HD and bHLH were identified. ITPK3 was identified as a candidate central gene for this module.