Integrated RNA-seq Analysis and Meta-QTLs Mapping Provide Insights into Cold Stress Response in Rice Seedling Roots

Rice (Oryza sativa L.) is a widely cultivated food crop around the world, especially in Asia. However, rice seedlings often suffer from cold stress, which affects their growth and yield. Here, RNA-seq analysis and Meta-QTLs mapping were performed to understand the molecular mechanisms underlying cold tolerance in the roots of 14-day-old seedlings of rice (RPY geng, cold-tolerant genotype). A total of 4779 of the differentially expressed genes (DEGs) were identified, including 2457 up-regulated and 2322 down-regulated DEGs. The GO, COG, KEEG, and Mapman enrichment results of DEGs revealed that DEGs are mainly involved in carbohydrate transport and metabolism, signal transduction mechanisms (plant hormone signal transduction), biosynthesis, transport and catabolism of secondary metabolites (phenylpropanoid biosynthesis), defense mechanisms, and large enzyme families mechanisms. Notably, the AP2/ERF-ERF, NAC, WRKY, MYB, C2H2, and bHLH transcription factors participated in rice’s cold–stress response and tolerance. On the other hand, we mapped the identified DEGs to 44 published cold–stress-related genes and 41 cold-tolerant Meta-QTLs regions. Of them, 12 DEGs were the published cold–stress-related genes and 418 DEGs fell into the cold-tolerant Meta-QTLs regions. In this study, the identified DEGs and the putative molecular regulatory network can provide insights for understanding the mechanism of cold stress tolerance in rice. In addition, DEGs in KEGG term-enriched terms or cold-tolerant Meta-QTLs will help to secure key candidate genes for further functional studies on the molecular mechanism of cold stress response in rice.


Introduction
Low temperature is a natural disaster that is often encountered in the entire process of rice growth and development, which is one of the main environmental factors limiting its growth and development, geographical distribution, yield, and quality formation [1,2]. The reduction of rice yield caused by cold stress has been a common problem in the world [2][3][4], especially in China, Japan, Australia, Korea, India, etc. [5][6][7].
The inheritance of cold tolerance in rice is complex and it is difficult to explain the mechanism of cold tolerance with single or several genes [6,7]. According to incomplete statistics, more than 300 QTLs were reported to be involved in cold tolerance using gene mapping studies of various populations at

Transcriptome Sequencing and DEGs Analysis
A total of 38.29 Gb of clean data was obtained from six samples ( Table 1). The percentage of Q30 (an error rate of sequencing lower than 1% o ) reached over 92.86% and GC contents were between 52.12% and 53.41%. The mapped read ratios of all samples ranged from 78.87% to 83.60%. Pearson's correlation coefficients of the three biological replicates for each treatment were all greater than 0.92 ( Figure 1A). These results indicated that sequencing data can be used in subsequent RNA-seq analysis. According to FDR< 0.01 and |log 2 fold change | > 2, we identified 2457 up-regulated and 2322 down-regulated DEGs in response to cold stress ( Figure 1B, 1C, Tables S1 and S2). Of them, 1481, 3858, 886, 1874, 4684, 3570, 3213, and 3905 DEGs were annotated in COG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and eggNOG databases.

Transcriptome Sequencing and DEGs Analysis
A total of 38.29 Gb of clean data was obtained from six samples ( Table 1). The percentage of Q30 (an error rate of sequencing lower than 1%o) reached over 92.86% and GC contents were between 52.12% and 53.41%. The mapped read ratios of all samples ranged from 78.87% to 83.60%. Pearson's correlation coefficients of the three biological replicates for each treatment were all greater than 0.92 ( Figure 1A). These results indicated that sequencing data can be used in subsequent RNA-seq analysis. According to FDR< 0.01 and |log2 fold change | > 2, we identified 2457 up-regulated and 2322 down-regulated DEGs in response to cold stress ( Figure 1B, 1C, Table S1, and S2). Of them, 1481, 3858, 886, 1874, 4684, 3570, 3213, and 3905 DEGs were annotated in COG, GO, KEGG, KOG, NR, Pfam, Swiss-Prot, and eggNOG databases.

Cold Stress-Related TFs and Salt-Stress-Affected Pathways
To get cold-responsive TFs in the cold stress regulatory network, various types of TFs and TRs were identified by iTAK software. The top six up-regulated TFs were mainly AP2/ERF-ERF, NAC, WRKY, MYB, C2H2, and bHLH family members ( Figure 3). However, the number and types of down-regulated TFs were less than those in up-regulated TFs. The top three down-regulated TFs belonged to MYB, NAC, and bHLH families ( Figure 3).
number of DEGs is distinguished by the size of the circle and the circle from blue to red represents the q-value from large to small.

Cold Stress-Related TFs and Salt-Stress-Affected Pathways
To get cold-responsive TFs in the cold stress regulatory network, various types of TFs and TRs were identified by iTAK software. The top six up-regulated TFs were mainly AP2/ERF-ERF, NAC, WRKY, MYB, C2H2, and bHLH family members ( Figure 3). However, the number and types of down-regulated TFs were less than those in up-regulated TFs. The top three down-regulated TFs belonged to MYB, NAC, and bHLH families ( Figure 3). We mapped DEGs to different functional categories in the MapMan tool to reveal the salt-stress-affected pathways. All DEGs were associate to 35 pathways ( Figure 4 and Table S5). They were significantly enriched on misc (8.62%), stress (5.96%), signaling (5.34%), secondary metabolism (3.06%), and hormone metabolism (2.70%) except for RNA (7.72%), protein (7.89%), and transport (4.65%) ( Figure 4). Therefore, these process genes were listed in detail in Figure 5. Mapman's identification of TFs was similar to that of iTAK software. Most AP2-EREBP, WRKY, and C2H2 TFs genes were up-regulated, while half of the bHLH and MYB TFs genes were down-regulated ( Figure  5A). Interestingly, four bZIP TFs genes showed obvious up-regulation under cold stress ( Figure 5A). A large number of large enzyme families are also differentially expressed under cold stress ( Figure  5B). This result indicated that cold stress severely affected the enzyme homeostasis in the rice roots. Multiple IAA, ABA, BA, Ethylene, Cytokinin, Jasmonate, SA, and GA pathways genes were DEGs ( Figure 5C), which also fully supported the KEGG enrichment result (ko04075, Figure S3). Many DEGs were related to protein modification and protein degradation ( Figure 5C). As expected, redox and signaling DEGs including G-proteins, MAP Kinases, Calcium regulation, etc., were found ( Figure 5C). In addition, we noticed a lot of secondary metabolism-related DEGs ( Figure 5D). These results proved that rice cold stress tolerance was a very complex regulatory process involving multiple pathways and transcription factors. We mapped DEGs to different functional categories in the MapMan tool to reveal the salt-stress-affected pathways. All DEGs were associate to 35 pathways ( Figure 4 and Table S5). They were significantly enriched on misc (8.62%), stress (5.96%), signaling (5.34%), secondary metabolism (3.06%), and hormone metabolism (2.70%) except for RNA (7.72%), protein (7.89%), and transport (4.65%) ( Figure 4). Therefore, these process genes were listed in detail in Figure 5. Mapman's identification of TFs was similar to that of iTAK software. Most AP2-EREBP, WRKY, and C2H2 TFs genes were up-regulated, while half of the bHLH and MYB TFs genes were down-regulated ( Figure 5A). Interestingly, four bZIP TFs genes showed obvious up-regulation under cold stress ( Figure 5A). A large number of large enzyme families are also differentially expressed under cold stress ( Figure 5B). This result indicated that cold stress severely affected the enzyme homeostasis in the rice roots. Multiple IAA, ABA, BA, Ethylene, Cytokinin, Jasmonate, SA, and GA pathways genes were DEGs ( Figure 5C), which also fully supported the KEGG enrichment result (ko04075, Figure S3). Many DEGs were related to protein modification and protein degradation ( Figure 5C). As expected, redox and signaling DEGs including G-proteins, MAP Kinases, Calcium regulation, etc., were found ( Figure 5C). In addition, we noticed a lot of secondary metabolism-related DEGs ( Figure 5D). These results proved that rice cold stress tolerance was a very complex regulatory process involving multiple pathways and transcription factors.  Table S4 have been provided in Table S5.   Table S4 have been provided in Table S5.   Table S4 have been provided in Table S5.  Each square in A-D represents a separate DEG, red indicates that gene expression was induced and blue indicates that gene expression was repressed compared with the control.

Key DEGs locking Involved in the Response to Cold Stress
Cold stress has a great influence on the growth and yield of rice, therefore, for many years, a large number of cold stress functional genes have been cloned and demonstrated through different rice materials, genetic population, and cloning methods, such as COLD1 [10], OsCDPK7 [29], OsWRKY45 [30], etc. We found that 12 DEGs in this study are genes that have been cloned (Table 2 and Figure 6). A total of 417 DEGs were dropped into 41 Meta-QTLs ( Figure 6, Table S6, and S7). Interestingly, EVM0001397 (SNAC2) and EVM0031981 (OsSPX1) were in MCqtl1-4 and MCqtl6-3, respectively. This result indicated that DEGs analysis combined with Meta-QTLs mapping can effectively reduce the number of candidate genes and accurately locate the key genes of cold stress. Forty-one Meta-QTLs contained 0 to 48 DEGs, and the number of DEGs was much smaller than the number of genes in Meta-QTLs. Of them, the number of DEGs in 36 Meta-QTLs was less than 20 ( Figure 6, Tables S6 and S7), especially, MCqtl9-3 (1 DEGs), MCqtl4-4 (2 DEGs), MCqtl12-4 (2 DEGs), MCqtl1-2 (3 DEGs), MCqtl7-2 (3 DEGs), MCqtl8-4 (3 DEGs), MCqtl2-2 (4 DEGs), MCqtl5-3 (4 DEGs), MCqtl10-1 (5 DEGs), and MCqtl11-4 (5 DEGs). EVM0009341 was the only DEGs in MCqtl9-3, annotated as Myb-related protein. This gene and the previously cloned cold-stress corresponding genes (Osmyb4 [31], MYBS3 [16], and OsMYB3R-2 [32,33]) belonged to the same family and showed a 3.63-fold upregulation ( Figure 6 and Table S7). MCqtl4-4 included EVM0032414 related to amino acid transport and metabolism and EVM0018896 (Function unknown). Two DEGs in MCqtl12-4 were annotated as Auxin-responsive protein IAA30 (EVM0025632, KEGG: ko04075) and CBL-interacting protein kinase 4 (EVM0013862). MCqtl1-2 also had a DEG in signal transduction mechanisms. Similarly, other Meta-QTLs also contained one or several potential cold-stress-responsive DEGs.  The green genes were cold-stress-related genes that have been cloned previously, but not the DEGs of this study. The red genes were the DEGs in this study and also belong to cold-stress-related genes previously reported. Black boxes in each chromosome (Chr) represent Meta-QTLs (e.g., MCqtl1-2, MCqtl1-3) and the number after Meta-QTL represents the number of our DEGs in this Meta-QTL segment. The green genes were cold-stress-related genes that have been cloned previously, but not the DEGs of this study. The red genes were the DEGs in this study and also belong to cold-stress-related genes previously reported. Black boxes in each chromosome (Chr) represent Meta-QTLs (e.g., MCqtl1-2, MCqtl1-3) and the number after Meta-QTL represents the number of our DEGs in this Meta-QTL segment.

Cold-Stress-Responsive DEGs in Rice Seedling Roots Associated with Plant Stress/Hormone Signal Transduction, Secondary Metabolites Mechanisms, TFs, and Function Proteins
Cold stress tolerance is a complicated quantitative trait that is controlled by various quantitative trait loci [3,15]. RNA-seq technology can identify a large number of cold-stress-responsive DEGs [24,34]. The obtained DEGs were systematically classified by various enrichment methods such as GO, COG, KEGG, and Mapman, which can fully analyze the cold stress regulation mechanism in rice roots. In this study, we obtained 4779 DEG-associated multiple biological processes. Lots of DEGs involved redox regulation, peroxidases, alcohol dehydrogenases, and glutathione S transferases, which implied that rice roots have suffered a certain degree of damage under cold stress. Previous studies proposed that cross-talk of plant hormones play important roles in the plant responses and resistance to multiple abiotic stresses [3,26]. Different plant hormones usually perform different functions, but their integration or coordinated expression is very important for plant growth and even survival under multiple stress environments [35,36]. Our KEGG enrichment term (ko04075, Figure 2C,D) and Mapman enrichment ( Figure 5C) revealed that most of plant hormones DEGs including IAA, ABA, BA, Ethylene, Cytokinin, Jasmonate (JA), SA, and GA, participated in the cold stress response of rice roots. In addition, DEGs in Ca 2+ signal transduction and MAPK cascade pathway were differently expressed ( Figure 5C), which supported earlier discovery that Ca 2+ signal transduction and MAPK cascade pathways were two major categories mediating stress signal transduction [24,25,34]. Interestingly, the expression levels of all JA DEGs were up-regulated. In Arabidopsis thaliana, JA enhanced the freezing tolerance by positively regulating the ICE-CBF/DREB1 pathway [37]. Therefore, we speculated that rice also enhances cold stress tolerance through this pathway. Wu et al. (2014) showed that cold acclimation is related to disease resistance in Amur grape (Vitis amurensis) [38]. Our DEGs were significantly enriched in plant-pathogen interactions (ko04626, Figure S2) and 34 pathogenesis-related proteins (PR-proteins) DEGs were found. These results suggested that the plant-pathogen interaction pathway may also play a role in cold resistance.
On the other hand, KEGG, GO enrichment, and COG function classification all pointed out that secondary metabolites play an important role in rice cold tolerance, especially phenylpropanoid biosynthesis (GO: 0009699, ko00940). Phenylpropanoids play a key role in the growth and development of plants and response to stress [39,40]. Our analysis result of DEGs showed that many secondary metabolites were widely involved in cold stress tolerance in rice. Earlier studies cloned some TFs that can enhance rice cold stress tolerance, namely, OsMYB2 [41], OsWRKY45 [30], Osmyb4 [42], etc. According to TF's identification result from iTAK and Mapman software, various AP2/ERF-ERF, NAC, WRKY, MYB, C2H2, and bHLH TF genes were up-/down-regulated. These DEGs, especially those that fall on Meta-QTLs, may be good candidates for transgenic genes to enhance cold stress tolerance in rice. Finally, all regulatory DEGs were summarized and presented through Mapman (Figure 7).

Meta-QTLs Mapping of DEGs Can Effectively Lock Main Effect Candidate Genes Related to Cold Stress Tolerance in Rice
Hu et al. (2008) reported that rice SNAC2 conferred cold and salt tolerance in rice [43]. Increased expression of OsSPX1 caused high sensitivity to cold and oxidative stresses in rice seedlings [44] and enhanced cold/subfreezing tolerance in tobacco and A. thaliana [45]. Yang et al. (2018) effectively reduced the physical distance MCqtl1-4 and MCqtl6-3 to 1.0 and 2.7 Mb based on Meta-QTLs analysis, respectively. However, most of the Meta-QTLs still contained many genes. Fortunately, we found that the possible genes in Meta-QTLs were greatly reduced after DEGs mapping of their Meta-QTLs. It should be noted that two previously cloned genes, SNAC2 and OsSPX1, were highlighted in MCqtl1-4 and MCqtl6-3 after our DEGs Meta-QTLs mapping. This result indicated that DEGs Meta-QTLs mapping is an effective method to discover the major genes in Meta-QTLs regarding cold stress tolerance. As expected, in Meta-QTLs containing very few DEGs, we have found possible candidate genes and these genes were supported by multiple annotation results. For example, EVM0013862 encoded CBL-interacting protein kinase 4 belonging to Ca 2+ signal transduction pathway in MCqtl12-4 (only two DEGs) and EVM0009341 encoding Myb-related protein was only a DEGs in MCqtl9-3. In this study, DEGs Meta-QTLs mapping made all 28 Meta-QTLs containing less than 15 genes (Tables S6 and S7).
However, MCqtl3-1, MCqtl12-3, MCqtl1-4, MCqtl10-2, MCqtl6-3, MCqtl7-4, and MCqtl10-4 contained 16-48 DEGs. Therefore, we provided comparison results of multiple databases to facilitate the further screening of candidate genes through annotations (Table S7). Of course, our screened genes in Meta-QTLs can be further fine mapped by developing SNP markers in flanking of these genes in the previously reported population, observing the mutant phenotypes of these genes from the published mutant library, or building CRISPR/Cas9 plants for cold stress experiments.
Int. J. Mol. Sci. 2020, 21, 4615 9 of 14 biosynthesis (GO: 0009699, ko00940). Phenylpropanoids play a key role in the growth and development of plants and response to stress [39,40]. Our analysis result of DEGs showed that many secondary metabolites were widely involved in cold stress tolerance in rice. Earlier studies cloned some TFs that can enhance rice cold stress tolerance, namely, OsMYB2 [41], OsWRKY45 [30], Osmyb4 [42], etc. According to TF's identification result from iTAK and Mapman software, various AP2/ERF-ERF, NAC, WRKY, MYB, C2H2, and bHLH TF genes were up-/down-regulated. These DEGs, especially those that fall on Meta-QTLs, may be good candidates for transgenic genes to enhance cold stress tolerance in rice. Finally, all regulatory DEGs were summarized and presented through Mapman (Figure 7).  Table S5.

Meta-QTLs Mapping of DEGs Can Effectively Lock Main Effect Candidate Genes Related to Cold Stress Tolerance in Rice
Hu et al. (2008) reported that rice SNAC2 conferred cold and salt tolerance in rice [43]. Increased expression of OsSPX1 caused high sensitivity to cold and oxidative stresses in rice seedlings [44] and enhanced cold/subfreezing tolerance in tobacco and A. thaliana [45]. Yang et al. (2018) effectively reduced the physical distance MCqtl1-4 and MCqtl6-3 to 1.0 and 2.7 Mb based on Meta-QTLs analysis, respectively. However, most of the Meta-QTLs still contained many genes. Fortunately, we found that the possible genes in Meta-QTLs were greatly reduced after DEGs mapping of their Meta-QTLs. It should be noted that two previously cloned genes, SNAC2 and OsSPX1, were highlighted in MCqtl1-4 and MCqtl6-3 after our DEGs Meta-QTLs mapping. This result indicated that DEGs Meta-QTLs mapping is an effective method to discover the major genes in Meta-QTLs regarding cold stress tolerance. As expected, in Meta-QTLs containing very few DEGs, we have found possible candidate genes and these genes were supported by multiple annotation results. For example, EVM0013862 encoded CBL-interacting protein kinase 4 belonging to Ca 2+ signal transduction pathway in MCqtl12-4 (only two DEGs) and EVM0009341 encoding Myb-related  Table S5.

Plant Materials and Treatments
Fourteen-day-old seedlings with uniform growth of rice cultivar 'RPY geng' (Oryza sativa ssp. japonica) were grown in 96-well PCR plate as supporting materials in Yoshida solution (Coolaber, Beijing, China) and changed to 5 • C (cold treatment)/26 • C (control), with a 16/8 h light/dark photoperiod, 60% relative humidity in plant growth incubators (ZSX1500GS, Jingshen Instrument, Shanghai, China) for three days, respectively. Here, three biological replicates were set using three separate cold/control treatments [20,46]. Only 1-2 roots from each seedling were collected immediately with liquid nitrogen for RNA extraction in order to include maximum number of plants (> 30 rice seedling) in single biological replicate [46,47].

Transcriptome Sequencing and RNA-seq Analysis
RNA extraction and cDNA library construction of the total samples were conducted by Biomarker Technologies (Beijing, China) as described previously [48][49][50]. Then, these cDNA libraries were sequenced on an Illumina sequencing platform (HiSeq 2500) [48,50]. All raw reads after sequencing were filtered using Trimmomatic and raw reads-containing adapter, ploy-N, and low-quality reads were removed. The clean data with high quality were used for the RNA-seq analyses. Q30, GC-content, and sequence duplication level of the clean data were calculated, respectively. HISAT2 software was used to map all clean reads to the reference genome of RPY geng (unpublished data) [51,52], and the mapped reads were assembled and quantified using StringTie software [53]. Gene expression levels were estimated by fragments per kilobase of transcript per million fragments mapped reads (FPKM) method [54,55]. Pearson's correlation coefficient of all tested samples was calculated based on all FPKMs. DESeq2 was used to identify differential expression genes (DEGs) with a False Discovery Rate (FDR) < 0.01 and |log 2 fold change | > 2 [34].

Functional Annotation and GO, KEGG Enrichment Analysis of the DEGs
Gene functions of DEGs were annotated by BlastP the following databases (E-value of e −5 ): Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (KEGG Ortholog database); GO (Gene Ontology).
Gene Ontology (GO) enrichment analysis of the DEGs was implemented by the GOseq R packages based Wallenius non-central hyper-geometric distribution [56], which can adjust for gene length bias in DEGs. GO terms of DEGs with KS < 0.01 were considered significantly enriched. We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways [57]. KEGG terms of DEGs with corrected P-value < 0.05 were considered significantly enriched [26].
The full-length nucleotide sequences of DEGs were submitted to the online annotation tool of Mapman (http://www.plabipd.de/portal/mercator-sequence-annotation) for functional annotation using default parameters [58]. Mapman annotations and gene expression levels of DEGs were visualized by Mapman [58].

Transcription Factors (TFs) and Transcriptional Regulators (TRs) Identification and Meta-QTLs Mapping.
We used iTAK software (http://itak.feilab.net/cgi-bin/itak/index.cgi) to identify TFs and TRs of DEGs with default parameters [59]. According to statistics from the National Rice Data Center (http://www.ricedata.cn/gene/) and the published papers, 44 cold-stress-related genes have been cloned. To determine how many of our DEGs are cloned genes, gene IDs of IRGSP-1.0 and nucleotide sequences of these genes were downloaded. Orthologs of our DEGs matched to genes in IRGSP-1.0 by blastn search using the following parameters: max_target_seqs 1, evalue e −5 , perc_identity 95 [34]. Then, DEGs matching the cloned genes were listed.
Meta-QTLs were obtained from Yang et al. study [13]. The physical location in IRGSP-1.0 of Meta-QTLs determined by online blast tool (https://rapdb.dna.affrc.go.jp/tools/blast) in the RAP-DB website (https://rapdb.dna.affrc.go.jp/download/irgsp1.html); the sequence of SSR markers must be identical with that of the whole genome sequence (100% identical, unique location). Then, the DEGs falling into the physical interval of Meta-QTLs were summarized, including our DEGs' ID, IRGSP-1.0, FPKM, function annotations, and other information.

Conclusions
In this study, we have identified numerous DEGs involving the putative molecular regulatory network of rice cold stress tolerance through the systematic RNA-seq analysis of the roots of fourteen-day rice seedlings. Our DEGs mapping to cold-tolerant Meta-QTLs previously reported provides effective target genes for future cold stress resistant rice breeding.  Table S6: Detailed information of Meta-QTLs; Table S7: Detailed information of DEGs in Meta-QTLs. Figure S1: KEGG term of phenylpropanoid biosynthesis (ko00940); Figure S2: KEGG term of plant hormone signal transduction (ko04075); Figure S3: KEGG term of plant-pathogen interaction (ko04626); Figure S4: KEGG term of circadian rhythm-plant (ko04712).

Conflicts of Interest:
The authors declare that they have no competing interests for this research.