Phenotypic changes of rice seedlings under low-temperature stress
Effect of low-temperature stress on osmoregulatory substances in rice
As illustrated in (Fig. 1B), the PRO content accumulated with cold treatment time and peaked at 48 h. The PRO accumulation rate and PRO content in JL were significantly higher than in JH. Compared with the control, after 4 h, 12 h, 24 h, and 48 h of cold treatment, the PRO content in JH was 1.98, 2.46, 4.05, and 4.69 times higher than that at the stress initiation, respectively. At this time, JL was 2.18, 3.19, 5.54, and 6.08 times higher than at the stress initiation. This indicates that the cold-tolerant variety, JL, can rapidly respond to low-temperature stress and maintain a relatively high PRO content (Fig. 1B).
During the 0–24 h of cold stress, the SS and SP contents of both varieties increased significantly, with a significant increase in JL. At 24 h, the SS and SP contents in the cold-sensitive variety, JH, were 1.61 and 1.93 times, respectively, higher than those at the beginning of the stress. At that time, JL was 1.76 and 2.44 times, respectively, higher than those at the beginning of the stress. At 24–48 h, the SS and SP contents in JL increased and remained at an elevated level; however, SS and SP contents in JH decreased sharply with the extension of chilling time. At 48 h, the contents of SS and SP in JH were 1.49 and 1.62 times the initial stress, and the contents of JL were 2.14 and 2.92 times the initial stress. At 24–48 h, SS and SP contents in JL still increased and remained at an elevated level, but the contents of SS and SP in JH decreased sharply with the extension of chilling time. At 48 h, the contents of SS and SP in JH were 1.49 and 1.62 times the initial stress, and those of JL were 2.14 and 2.92 times the initial stress. The cold-tolerant variety—JL—could enhance its cold resistance by accumulating SS and SP to maintain them high during low-temperature stress. However, the SS and SP contents of the cold-sensitive variety—JH—gradually decreased after 24 h as cold stress continued to cause damage to the plants (Fig. 1B).
Effect of low-temperature stress on antioxidant enzymes in rice
The CAT and POD of both species increased continuously with cold treatment time and peaked at 48 h. However, the accumulation rates and contents of CAT and POD in JL were significantly higher than those in JH. Compared with the control, after 4 h, 12 h, 24 h, and 48 h of cold treatment, the CAT content in JH was 1.09, 1.19, 1.37, and 1.41 times that at the beginning of the stress, respectively, when JL was 1.07, 1.26, 1.42 and 1.72 times that at the stress initiation, respectively. POD in JH was 1.12, 1.24, 1.49, and 1.60 times at this time, and 1.23, 1.40, 1.83, and 2.16 times at this time in JL than at the stress beginning. These results indicate that the cold-tolerant variety JL can accumulate CAT and POD rapidly in response to low-temperature stress and maintain relatively high levels to remove H2O2 from leaf tissues and inhibit peroxides damage to plants (Fig. 1B).
During 0–24 h of cold stress, T-AOC and SOD contents of both varieties increased significantly, with a greater increase in the JL. At 24 h, T-AOC and SOD contents in the cold-sensitive variety—JH—were 3.14 and 2.87 times higher than at the beginning of the stress, while those in the cold-tolerant variety—JL—were 4.11 and 3.24 times higher than at the beginning of the stress. Furthermore, at 24–48 h, the T-AOC and SOD contents in JL increased and remained at high levels. Nonetheless, the T-AOC and SOD times in JH decreased sharply with prolonged cold damage. At 48 h, T-AOC and SOD in JH were 2.63 and 2.65 times higher than at the beginning of the stress and 5.16 and 4.34 times higher than at the beginning of the stress in JL. This indicates that the cold-tolerant variety—JL—can remove excessive reactive oxygen radicals and enhance the cold resistance of plants by maintaining high SOD levels (Fig. 1B).
Effects of low-temperature stress on oxidative stress of rice
During cold stress, the OH- and MDA contents in both species increased continuously, and the accumulation rates and contents of OH- and MDA contents in JH were significantly higher than those in JL. Compared with the control, the OH- content in JH was 1.14, 1.21, 1.57, and 1.76 times higher than at the beginning of the stress after 4 h, 12 h, 24 h, and 48 h of cold treatment. JL was 1.06, 1.15, 1.31, and 1.47 times higher than at the beginning of the stress, respectively. The MDA content in JH was 2.22, 3.25, 5.16, and 6.51 times higher than at the beginning of the stress. On the other hand, JL was 1.41, 2.43, 3.04, and 3.49 times higher than at the beginning of the stress, indicating that the cold-tolerant variety—JL—could scavenge OH- and MDA during the-low temperature stress; thus, maintaining OH- and MDA at a low level (Fig. 1B).
Identification of important antioxidant phenotypes
Variable importance in projection (VIP) values are orthogonal partial least squares discriminant analysis variable projection importance values that can measure the contribution of individual phenotypic changes to the classification, with larger VIP values contributing more. Usually, a VIP value >1 is a screening criterion (Lamine et al. 2018). Orthogonal partial least squares discriminant analysis was performed at different time points of cold treatment based on the nine major metabolites in response to low temperature, according to VIP scores. The VIP scores of POD and SOD were >1 in JH. The VIP scores of POD, SOD, and OH- were >1 in JL, where SOD and POD were in both varieties and identified as the most essential antioxidant phenotypes among them (Fig. 1C, D).
High throughput sequencing
To systematically identify the expression profiles of lncRNAs and mRNAs associated with cold tolerance in Japonica rice seedlings, we performed strand-specific whole transcriptome RNA-seq on control and cold-treated rice seedling leaves of two Japonica varieties using Illumina machines. All samples (2 varieties ×5 time points ×3 replicates) were analyzed in three independent biological replicates. After removing low-quality and contaminated reads, all samples from JH and JL yielded nearly 3.3 billion high-quality reads (approximately 100 million per sample on average). The mapping of clean reads in each library to the Oryza sativa reference genome was 95.97–97.86%. In addition, all samples had Q20 values >97.8%, Q30 values >93.64%, and GC content of 46.55–48.21% (Additional file 2: Table S1), indicating the reliability of the RNA-seq data. These results provide a wealth of data for further analysis of the expression profiles and metabolic pathways of mRNAs and lncRNAs.
Comparative analysis of lncRNAs and mRNAs
These identified mRNAs and lncRNAs are widely distributed on all chromosomes, among which the number of LncRNAs and mRNAs on chromosome 1 is the largest. In contrast, LncRNAs and mRNAs on chromosomes 9, 10, and 12 are relatively few. The number of mRNAs not on a chromosome is more than lncRNA (Fig. 2A). Moreso, we conducted a conservative analysis of lncRNAs, with E value <1e-5 as the standard. The results revealed that 4, 1, 12, 62, and 127 lncRNAs were sequence conserved with Arabidopsis thaliana, soybean, Terrestris alfalfa, sorghum, and maize, respectively. Therefore, this result further indicated that lncRNAs were less conserved among different organisms. In addition, we further analyzed the basic characteristics of identified lncRNAs and mRNAs in rice, as presented in (Fig. 2B). The length of lncRNAs were 200–200538 nt, with an average of 1362 nt. Most lncRNAs were <2000 nt in length, with the most lncRNAs <1000 nt in length, in decreasing order. On the other hand, mRNA lengths were 30–57648 nt, with an average of 3038 nt, and their lengths were mainly concentrated within 3000 nt. Most lncRNAs had only one exon, while the proportion of mRNAs with over two exons was significantly higher than that of lncRNAs (Fig. 2C). These results suggest that mRNAs and lncRNAs play different roles in biological pathways.
Identification of differentially expressed mRNAs
The inter-sample agreement was examined using the Pearson correlation coefficient (R-value) between both varieties. The minimum R-value for each variety was 91.49%, and the maximum was 99.31% for three biological replicates of each variety compared. Replicate samples revealed close correlations, proving the reproducibility of the samples and their continued use for subsequent analysis (Additional file 1: Fig. S1). We compared the clean reads of each sample with the reference genome and obtained 36,850 mRNAs. According to the screening criteria of Fold Change ≥ 2 and FDR ≤0.01, 12,817 and 15,806 mRNAs were identified from JH and JL, respectively, compared with the control 0 h, among which 8,088 genes in JH were significantly up-regulated, 4,729 were significantly down-regulated. Additionally, 10,446 genes in JL were significantly up-regulated, and 5,363 were significantly down-regulated. The number of up-regulated differential mRNAs was much larger than that of down-regulated differential mRNAs. With prolonged cold treatment time, the number of differential mRNAs gradually increased, indicating the changed gene expression of rice seedlings with prolonged cold treatment time (Fig. 3A).
We also compared both varieties at the same time point. We observed that 3,102 genes (including 70 TF-encoding genes) were highly expressed, and 3,028 genes (including 82 TF-encoding genes) displayed significantly low expression than JH. Similarly, the number of differentially expressed genes in the varieties at 48 h was the highest (Fig. 3B). Some TF families displayed significant differences in the number of significantly high and low expressed members. For instance, the expression of TF family members involved in abiotic stresses in plants, such as MYB, HB, and NAC, was significantly elevated. On the other hand, transcription factors with lower expression were mainly from the WRKY, bZIP, and MADS families (Fig. 3C). The expression profiles of these TF families are presented in (Additional file 1: Fig. S2). The volcano plot (Additional file 1: Fig. S3) indicates the gene expression profiles of DEmRNAs in the 13 comparative combinations, with a clear distribution of up/down-regulated genes, indicating different expression patterns in the different groups.
Identification of differentially expressed lncRNA
To characterize putative lncRNAs of rice, first, the transcripts with length <200 nt were filtered. Second, the coding potential of the remaining transcripts was evaluated using Rfam, Pfam, and CPC2, and the intersection of non-protein-coding transcripts predicted using the three software was candidate lncRNAs for further analysis. Finally, 7,087 lncRNAs were identified in JH and JL. Based on the position of lncRNAs in relation to protein-coding genes, we classified them into four types: (3,035, 42.82%) antisense lncRNA; (2,700, 38.1%) intergenic lncRNAs; (548, 7.74%) intronic lncRNAs; (804. 11.36%) sense lncRNAs (Fig. 4A).
For each transcribed region, transcribed fragments per kilobase per million reads (FPKM) values can be estimated using StringTie software to quantify its expression abundance and variation. The FPKM distribution of the 10 treatments revealed their basic expression patterns. As illustrated in (Additional file 5: Fig. S4), the density-FPKM curves revealed the same trend in all samples, indicating that lncRNA expression was comparable among the groups.
According to the screening criteria of Fold Change ≥ 2 and FDR ≤0.01, we first identified differential genes at each cold treatment time point of each variety. Compared with the control 0 h, 1,306 and 2,414 differential lncRNAs were identified in JH and JL, respectively. Additionally, 75.1% of the genes in JH were significantly up-regulated, and 24.9% were significantly down-regulated. In JL, 59.4% of the genes were significantly up-regulated, and 40.6% were significantly down-regulated (Fig. 4B). In comparisons, the total number of significant regulatory genes differed. However, the number of up-regulated genes was higher than that of down-regulated genes in all groups. The number of differential lncRNAs in the late cold-treatment period was significantly higher than that in the early period, suggesting to some extent that these DElncRNAs may be involved in the response of rice to low temperature.
At the same time, we also compared both varieties at the same time point. Compared with JH, 726 genes were significantly highly expressed, and 951 genes had significantly lower expression. The differential JL and JH genes were the highest at 4 h and lowest at 12 h (Fig. 4B).
Identification of differentially expressed miRNAs
We also performed small RNA sequencing on all samples and obtained 409055593 clean reads. As a result, 542 existing miRNAs were identified using Bowtie (version 1.1.2) software, of which 409 were known miRNAs, and 113 were novel miRNAs. Most of them were 21 nt (49.82%) and 22 nt (20.30%) in length (Fig. 5A). The nucleotide composition of mature miRNAs indicated that two clusters of highly abundant miRNAs (21 nt and 24 nt) had a significant sequence bias, in which 5'U was dominant (Fig. 5B). Using edgeR software, with Fold Change ≥ 1.5 and P-value ≤ 0.05 as the screening criterion, there were (41, 32, 55, 24), (33, 34, 27, 36) differential miRNAs in JH and JL at 0 h vs. 4 h, 0 h vs. 12 h, 0 h vs. 24 h, and 0 h vs. 48 h, respectively. Also, comparing both species (78, 86, 37, 78, 64) differential miRNAs in JH vs. JL at 0 h, 4 h, 12 h, 24 h, and 48 h, respectively (Fig. 5C), and based on these data, scatter plot analysis visualized the miRNA differences between the compared groups.
Identification of weighted gene co-expression network analysis modules associated with premature aging
To comprehensively analyze the relationship between differential genes and cold tolerance phenotypes, we performed a weighted gene co-expression network analysis (WGCNA). A gene clustering scheme with a power of 7 was constructed. Genes with similar expression patterns after filtering were divided into 17 modules (Fig. 6A). Each module’s gene expression profile was correlated with all samples to generate the heat map of the mode sample matrix (Fig. 6B). For instance, the gene expression profile of the black module (r=-0.90, P=4e×10-4) positively correlated with SOD. On the other hand, the gene expression profile of the green module (r=-0.85, P=0.005) was significantly negatively correlated with POD (Fig. 6B), indicating that the genes in the green module may negatively regulate the cold tolerance of rice seedlings. SOD and POD were antioxidant phenotypes with VIP scores >1 (Fig. 1B, C), further proving the reliability of selecting both modules.
Identification and analysis of key genes for cold tolerance
There were 1,130 expressed genes in the black module and 2,253 expressed genes in the green. According to TopGO enrichment analysis, genes in key co-expression modules are significantly enriched in numerous pathways that can be classified into three major categories: biological processes, cellular components, and molecular functions.
Genes in the black module were significantly enriched in terpene synthase activity, oxidoreductase activity, acting on the CH-NH group of donors, oxygen as acceptor, abscisic acid binding, protein phosphatase inhibitor activity, glutathione transferase activity, magnesium ion binding, anthocyanin-containing compound biosynthetic process, regulation of protein serine/threonine phosphatase activity, glutathione metabolic process, toxin catabolic process (Additional file 1: Fig. S5). Genes in the green module were significantly enriched in ADP binding, polynucleotide adenylyltransferase activity, RNA-directed DNA polymerase activity, iron ion binding, mitochondrion, plastid, DNA integration, RNA-dependent DNA replication, apoptotic process, mRNA polyadenylation (Additional file 1: Fig. S6). We also performed KEGG enrichment analysis, selecting the top 10 pathways in terms of the p-value. The genes in the black module were mainly concentrated in flavonoid biosynthesis, galactose metabolism, stilbenoid, diarylheptanoid, gingerol biosynthesis, endocytosis, arginine, and PRO metabolism, phosphatidylinositol signaling system, beta-Alanine metabolism, other types of O-glycan biosynthesis, selenocompound metabolism, phenylalanine metabolism. The genes in the green module are mainly concentrated in the linoleic acid metabolism, valine, leucine, and isoleucine biosynthesis, homologous recombination, pyrimidine metabolism, purine metabolism, pantothenate and CoA biosynthesis, pyruvate metabolism, 2-Oxocarboxylic acid metabolism, nucleotide excision repair, mismatch repair (Fig. 7A, B).
To identify hub genes controlling cold tolerance in rice, linkage count analysis was performed on genes within the key module, and genes with high linkage counts were considered very important. The top five genes with the highest connectivity in each module were selected as hub genes. Hub genes in the black module included a UDP-glucosyltransferase family protein (Os05g0527000), a sesquiterpene synthase (Os08g0167800), an indole-3-glycerophosphatase gene (Os03g0797500), a gene encoding oxidoreductase (Os04g0339400) and an unknown gene (Os05g0212900). The hub genes in the green module include a WRKY transcription factor (Os11g0685700), a plasma membrane hydrogen-ATPase (Os02g0825600) related to the ABA signaling pathway, a leucine-rich repeat-like receptor kinase containing a Beta-Ig-H3 structural domain protein (Os02g0615800), and two unknown genes; Os03g0103950 and Os08g0288050.
Transcription factors play an important role in regulating gene expression in plant cells, for which we analyzed differential transcription factors in the module. A WRKY family transcription factor (Os05g0137500) is present in the black module, and a number of differentially expressed typical (including cold-tolerant) abiotic stress related genes were identified during transcriptome sequencing of rice seedlings of cold-tolerant varieties, including this WRKY transcription factor (Guan et al. 2019), a LOB family transcription factor (Os01g0511000), one of the top 50 genes down-regulated under low temperature conditions during booting stage in rice according to microarray analysis, A DBB family transcription factor (Os04g0540200), a B3 family transcription factor (Os06g0112300), and a bHLH family transcription factor (Os09g0410700).
In addition to the hubWRKY transcription factor (Os11g0685700) in the module, three WRKY family transcription factors (Os01g0656400, Os05g0571200, and Os11g0686250) are also present in the green module when the transcriptome of drought-sensitive rice was sequenced under drought, Os05g0571200 was the hub gene under drought conditions (Yu et al. 2020), three WYB transcription factors (Os05g0490600, Os06g0605750, and Os07g0137000), of which Os07g0137000 in rice at late stages of heat stress was induced to be differentially expressed (Hernández-Soto et al. 2021). Four MADS (Os01g0726400, Os08g0112700, Os12g0501700, and Os03g0752800), two B3 (Os03g0619600 and Os03g0621600), two bZIP transcription factors (Os06g0614100, Os12g0601800), one AP2 (Os04g0572200), and one EIL transcription factor (Os03g0324200). Os03g0324200 is an EIL family transcription factor that regulates the ethylene signal transduction pathway. Overexpression of the ethylene response factor—TERF2—enhanced cold tolerance in rice seedlings (Tian et al. 2011), GRAS transcription factor (Os03g0103400), a GRF transcription factor (Os03g0674700), an HB transcription factor (Os03g0673000), an HSF transcription factor (Os02g0232000), and a Trihelix transcription factor (Os11g0163500) was observed.
Construction and analysis of lncRNA-miRNA-mRNA ceRNA network
Using TargetFinder to compare the sequences of all types of RNAs obtained by sequencing, we observed relationships between 1,138 lncRNAs, 5,601 mRNAs, and 225 miRNAs (Additional file 2: Table S2). To accurately analyze the role of these RNAs during the cold rice response, we intersected the fraction of differences in these relationship pairs with genes in the significantly related module in WGCNA. In total, 85 intersecting DEmRNAs were obtained in the black module (Additional file 2: Table S3), and 461 intersecting DEmRNAs were obtained in the green module (Additional file 2: Table S4).In the DEmRNAs intersection of the black module, we selected some hub genes of the module screened above and the redox process differential genes mainly related to cold tolerance in the enrichment pathway. The network consisted of 5 mRNAs, 8 miRNAs, and 11 lncRNAs (Fig. 8A). Os08g0532600 was predicted as the target gene of osa-miR162b, osa-miR2905, and Novel.58. On the other hand, Os12g0448900 was predicted as the target gene of Novel.65 and osa-miR156f-3p. Furthermore, Os08g0167800 was the target gene of osa-miR396d, Os04g0339400 was the target gene of Novel.44, and Os09g0368200 was predicted as the target gene of osa-miR1865-5p. Among them, the gene encoding REDOX enzyme (Os04g0339400), peroxidase family protein (Os12g0448900), similar to peroxidase (Os08g0532600), polyamine oxidase (Os09g0368200) are essential in redox.
In the intersection DEmRNAs of the green module, we also selected the hub gene of the module screened above and the genes related to linoleic acid metabolism mainly related to cold tolerance in the enrichment pathway. The network consisted of 6 mRNAs, 9 miRNAs, and 11 lncRNAs (Fig. 8B). Os02g0825600 was the target gene of osa-miR2873b and osa-miR1861c, Os06g0215900 was predicted as the target gene of Novel.65, Os03g0699700 was predicted to be the target gene of osa-miR1861d and osa-miR812o-5p, OsGLR1.2 was predicted as the target gene of osa-miR5493 and osa-miR810b.1. Additionally, Os02g0615800 was the target gene of osa-miR390-5p and Os08g0508800 is predicted to be the target gene of Novel.50. Among them, 12-oxydienate reductase (Os06g0215900) and lipoxygenase gene (Os08g0508800 and Os03g0699700) are the key genes in linoleic acid metabolism.
We observed common lncRNAs (LNC2186, LNC0595) and miRNAs (Novel.65) in the green and black modules. Finally, we integrated the networks where the common lncRNAs and miRNAs were located in both modules and constructed the ceRNA network using Cytoscape3.3.0 software. As a result, there is a complex relationship between these ceRNAs, as presented in (Additional file 1: Fig. S7).
Validation of differentially expressed genes in RNA-Seq using qRT-PCR
We compared the log2fold change of six selected DEGs between RNA-Seq and qRT-PCR to validate the RNA-Seq results, demonstrating a correlation coefficient (R2) of 0.8881 (Additional file 1: Fig. S8). The qRT-PCR data were consistent with the RNA-Seq data, indicating the reliability of the RNA-Seq results.
Sequencing analysis of OsGLR1.2 for identification, subcellular localization, and mutant construction
To further confirm the function of OsGLR1.2 in affecting cold tolerance in rice, the seedling cold-tolerant Nipponbare genome was used as the reference sequence, and the promoter and CDS sequences of OsGLR1.2 in JH and JL with significantly different cold tolerance were amplified using PCR and sequenced separately. Compared with the Nipponbare sequence, there were four SNPs and one Indel on CDS; the CDS sequence of JL was completely consistent with the CDS sequence of Nipponbare. Four non-synonymous amino acid substitutions and one code-shifting mutation were found in the JH exon sequence compared with the CDS sequence of Nipponbare (Fig. 9). The sequence of the JL promoter region was completely consistent with that of the Nipponbare sequence, and there were several variations in the JH promoter region, including 18 SNPS and 4 Indels (Additional file 1: Fig. S9). These variation sites on the promoter and CDS may be the cause of differences in cold tolerance.
OsGLR1.2 is localized in the cell membrane and Golgi apparatus
To determine the subcellular localization of OsGLR1.2, we expressed the GFP-OSGLR1.2 fusion construct in rice protoplasts and found that the control vector (35S: GFP) had a bright GFP signal green distributed throughout the whole cell. However, the fluorescence of 35S: GFP: OsGLR1.2 fusion protein was only localized in the cell membrane and Golgi apparatus (Fig. 10), indicating that OsGLR1.2 may be localized in the cell membrane and Golgi apparatus.
Mutant construction and phenotype analysis
To investigate the effect of OsGLR1.2 knockdown on rice growth, we used the CRISPR-Cas9 method to construct OsGLR1.2 mutants. A CRISPR/Cas9 vector targeting OsGLR1.2 coding region guide RNA was first constructed (Fig. 11A) and transformed into rice (DN428). Then, PCR amplification and sequencing were performed for each generation of regenerated knockout plant leaves, and three mutation types were screened for T3 generation pure plants with 1 bp insertion, 1 bp deletion, and 5 bp at the first exon target site, respectively (Fig. 11B), resulting in the early termination of code shifting and translation. We also used online tools (http://skl.scau.edu.cn/offtarget/) in the rice genome search and OsGLR1.2 to target highly similar sequences to test potential off-target cutting, with miss cut not detected.
SOD, POD, CAT, SS, PRO, and MDA are the key indicators of the physiological state of cells under cold stress. Compared with the control (room temperature treatment), SOD, POD, CAT, SS, PRO, and MDA contents were increased after 4 h, 12 h, 24 h, and 48 h of cold treatment, but the knockout strain had significantly higher SOD, POD, CAT, SS, and PRO and significantly lower MDA contents than the wild type (WT) (Fig. 11D–I). Meanwhile, we observed that the knockout strain was phenotypically affected by cold compared with the WT (Fig. 11C). The results showed that the knockdown of OsGLR1.2 improved the cold tolerance of rice seedlings.