An association between the nitrogen metabolism pathway and cold tolerance in rice was identified using comparative transcriptome and proteome profiling

Rice ( Oryza sativa L.), one of the most important crops cultivated in both tropical and temperate regions, has a high sensitivity to cold stress. Chilling stress limits the N uptake and nitrogen metabolism in rice. To identify the genes and pathways involved in cold tolerance, and specifically associations with the nitrogen metabolism pathway, we have compared the gene and protein expression changes between a cold-tolerant cultivar, Dongnong428 (DN), and a cold-sensitive cultivar, Songjing10(SJ). Using absolute quantification (iTRAQ) with high-throughput mRNA sequencing (RNA-seq) techniques, we identified 5,549 genes and 450 proteins in DN and 6,145 genes and 790 proteins in SJ, that were differentially expressed during low- water temperature (T w ) treatment. There were 354 transcription factor (TF) genes (212down, 142 up), 366 TF genes (220 down, 146 up), including 47 gene families, differentially expressed in the DN under control (CKDN) vs. DN under low-T w (D15DN) and CKSJ vs. D15SJ, respectively. These results indicated that TF genes play a major role in post-translational regulations. Genes related to rice cold-related biosynthesis pathways, particularly the MAPK signaling pathway, zeatin biosynthesis, and plant hormone signal transduction pathways, were significantly differentially expressed in both rice cultivars. Differentially expressed proteins (DEPs) related to rice cold-related biosynthesis pathways and particularly glutathione metabolism were significantly differentially expressed in both rice cultivars. Transcriptome and proteome analysis of the nitrogen metabolism pathways showed that major genes and proteins were down-regulated that participated in γ-aminobutyric acid (GABA) and glutamine synthesis. CPS2; AIG2LD; uncharacterized isomerase basic secretory protease GO terms oxidoreductase activity; terpenoid biosynthetic process; lipid metabolic process; phytoalexin metabolic process; oxidation-reduction process. GO terms of magnesium ion binding; lyase activity; terpenoid biosynthetic process; lipid transferase gene involved in L-ornithine synthesis were up-regulated in the D15SJ-vs-CKSJ group. In the D15DN-vs-CKDN group,2 glutamate synthase genes involved in glutamate synthesis were up-regulated. Compared to the D15SJ, 1 ornithine carbamoyl transferase gene involved in L-ornithine synthesis, showed significant up -regulation in D15DN. regulation of transcription factors. The present study confirmed the known cold stress- associated genes and identified a number of putative


reduced.
We selected the cold-tolerant cultivar Dongnong428 (DN) and cold-sensitive cultivar Songjing 10 cultivar (SJ) to utilize in the following experiments (Fig. 1a). The yield of CKDN was 880 g/m 2 , that of CKSJ was 810 g/m 2 . The 15 d of low-T w treatment decreased the yields in comparison to the controls, by 41% for D15DN (515 g/m 2 ) and 51% for D15SJ (395 g/m 2 ).
We measured the physiological indicators and nitrogen content of the samples that experienced chilling damage. Compared with the untreated samples, the nitrate reductase (NR) activity of the root was significantly decreased in D15DN (37.8 %) and D15SJ (51.5 %) , while the glutamate synthase (GOGAT), glutamine synthetase (GS), glutamate dehydrogenase (GDH), aspartate aminotransferase (GOT), and alanine aminotransferase (GPT) activities were significantly increased in D15DN and D15SJ. This indicated that under the low-T w treatment, the NO 3 reducing ability in the roots of the rice was decreased under low-T w treatment, but the NH 4 + assimilation ability was enhanced; which was the main reason for the increased of N concentration in the roots (21.6%, DN; 14.1%, SJ). However, due to the low-T w treatment, the biomass in the roots of the rice was significantly decreased (30.4%, DN; 36.4%,SJ), which leads to decreased N content in the roots (15.4%,DN; 27.4%,SJ) (Fig. 1b).
Amino acid content analysis showed that the glutamate acid (Glu) content in D15DN was significantly reduced and the Pro content was slightly increased, compared to the CKDN.
In D15SJ, the content of Glu increased and the content of Pro decreased slightly, compared to the CKSJ. Compared with D15SJ, the Glu content in D15DN was 64.7% less, while the Pro content was increased by 42.9% (Fig. 1c).

Genome-wide expression changes of mRNA and protein expression under cold stress
There were significant phenotypic differences in the dry matter accumulation and N metabolism related physiological indicators of the roots between the DN and SJ (Fig.1b). In order to reveal the regulatory mechanisms of the cold stress and explore new regulators, cold responsive transcriptome and proteome changes in the roots of DN and SJ were studied using RNA-seq and iTRAQ. Cold-tolerant cultivar DN and cold-sensitive cultivar SJ were grown normally (Control; CK) and in low-T w conditions (D15  1993). There were 1,184,867 spectra generated from the roots of the two rice varieties. A total of 265,517 of the spectra were matched to known spectra, and 212,755 spectra were matched to unique spectra. There were also 11,666 mapped peptides, 32,750 mapped unique peptides, and 7,281 mapped proteins, with a 1% false discovery rate (FDR). The repeatability of the proteomic analyses according to the CV (coefficient of variation) revealed that 95-100% of the proteins could be found in all the identified proteins, when the CV value was less than 50 % among the replicates.
Of the proteins, 88.8% were detected based on iTRAQ, and corresponding transcripts could be detected in the transcriptome, accounting for 14.8% of the total number of transcripts detected (Fig. 2a). To explore the relationship between the proteins and their corresponding genes, we matched all expressed proteins with their corresponding genes in each treatment and observed a weak correlation, with r-values ranging from 0.14 to 0.31, demonstrating that gene expression cannot fully account for the abundance of protein species, and strong post-translational regulations may exist in the process of protein species production.

Cold response transcriptome changes in the two rice cultivars
Differential gene expression was examined using the R package DEseq2 to quantify and analyze all control and treatment combinations (Love et al, 2014).Over the 15 days of cold exposure, a total of 5,549 genes were differentially expressed in the roots between D15DN and CKDN of which 3,539 genes were down-regulated and 2,010 up-regulated. There were 6,145 genes that were differentially regulated in the roots between D15SJ and CKSJ of which 3,690 were down-regulated and 2,455 were up-regulated (Fig.2b).The number of differentially expressed genes in the two varieties of rice was close after the cold treatments, indicating that there is some similarity between the two varieties in their responses to the cold treatments.
To classify the plant systems influenced during the cold stress, the DEGs were searched against the GO database for enrichment analysis. The GO database enrichment analysis was separated into three ontologies: molecular function (MF), cellular component (CC), and biological process (BP

Cold response proteome changes in contrasting rice cultivars
Applying the cut-off threshold of a 1.2-fold change for increased accumulations and a 0.83-fold change for decreased accumulations, together with the number of unique peptides ≥2, a total of 450 proteins were differentially expressed in the roots between D15DN and CKDN of which 215 were down-regulated and 235 were up-regulated; whereas 790 proteins were differentially regulated in the roots between D15SJ and CKSJ, and 454 were down-regulated and 336 were up-regulated (Fig. 2c).
Gene Ontology enrichment analysis was used to identify the over-represented terms from external encapsulating structure and plasma membrane (Fig. 3b).
The biological process terms associated with the DEPs of CKDN vs. D15DN and CKSJ vs.
D15SJ enriched terms both focused onoxidation-reduction process; single-organism metabolic process; diterpenephytoalexin metabolic process; terpenoid metabolic process; generation of precursor metabolites and energy; carbohydrate metabolic process; response to cold; response to oxidative stress. Specific terms of enrichment for CKDN vs.
D15DN were glycogen biosynthetic process, and for CKSJ VS. D15SJ there was one term, response to cadmium ion (Fig. 3b).
The KEGG pathways enriched using the DEPs in CKDN vs. D15DN and CKSJ vs. D15SJ were diterpenoid biosynthesis, phenylpropanoid biosynthesis, flavonoid biosynthesis, cutin, suberin, and wax biosynthesis, amino sugar and nucleotide sugar metabolism; biosynthesis of amino acids; ascorbate and aldarate metabolism, carbon metabolism, and flavone and flavonol biosynthesis. The specific pathway enrichments for the CKDN vs.
There have been previous report on the genes involved in diterpenoid biosynthesis and the phenylpropanoid biosynthesis pathways, in response to cold stress (Fang et al, 2017).

Combined transcriptome and proteome analyses of the differences in cold response mechanisms of the two rice cultivars
Venn analysis showed that in the transcriptome data, 1015 DEGs were up-regulated in the CKDN \D15DN group, and 904 were down-regulated (Fig. 4a, b).Venn analysis showed that in the proteome data, only 38 DEPs were up-regulated in the CKDN\D15DN group, and 50 were down-regulated ( Fig. 5a, b). There were only 5 up-regulated and 9 down-regulated in both proteome and transcriptome data in the CKDN \D15DN group, indicating a strong post-translational regulation and the complementarity of the transcriptome and proteome analysis. The up-regulated genes contained an aspartyl protease family protein(Os06g0304600); tuliposide A-converting enzyme b6(Os09g0455900); probable polygalacturonase (Os06g0106800); and enolase 1(Os06g0136600); lichenase-2(Os05g0375400). The down-regulated genes contained FAR1; PRXIIE-2; CYP76M7; CXE15; CPS2; AIG2LD; ALDH2C4; uncharacterized isomerase BH0283 (Os01g0266500); and basic secretory protease (Os10g0491000). CYP76M7 participated in the GO terms of oxidoreductase activity; terpenoid biosynthetic process; lipid metabolic process; phytoalexin metabolic process; oxidation-reduction process. CPS2 participated in the GO terms of magnesium ion binding; lyase activity; terpenoid biosynthetic process; lipid metabolic process; and phytoalexin metabolic process. Enolase 1 participated in the GO terms of magnesium ion binding; lyase activity; oxidoreductase activity; oxidationreduction process and CYP76M7, CPS2, and enolase 1 may play important roles in resisting cold stress.

Identification of transcription factors (TFs) responding to cold response
27.4% in the SJ. The amino acid content analysis showed that the glutamate content in the D15DN was significantly reduced and the Pro content was slightly increased, compared to the control. In D15SJ, the content of glutamic acid increased, and the content of Pro decreased slightly, compared to the control. Compared with D15SJ, the glutamate content in D15DN was 64.7% less, while the Pro increased by 42.9%.
Transcriptome analysis of the nitrogen metabolism pathways, identified 9 genes that were down-regulation in the D15SJ-vs-CKSJ group (Fig. 7a, b), including 3 glutamate decarboxylase genes involved in GABA synthesis, 2 glutamine synthetase genes involved in glutamine synthesis, 1 glutamate synthase gene involved in glutamate synthesis from glutamine, 2 delta-1-pyrroline-5-carboxylate synthetase genes involved in L-Glutamyl 5phosphate synthesis, and 1 amino-acid N-acetyltransferase gene involved in in N-Acetyl-Lglutamate synthesis. There were 6 genes that were down-regulation in the D15DN-vs-CKDN group, including 3 glutamate decarboxylase genes involved in GABA synthesis, 2 glutamine synthetase genes involved in glutamine synthesis, and 1 delta-1-pyrroline-5carboxylate synthetase gene involved in L-glutamyl 5-phosphate synthesis. Compared to the D15SJ, no genes showed significant up-or down-regulation in the D15DN.
Proteome analysis of the nitrogen metabolism pathways, identified 5 proteins that were down-regulation or down-regulation in the D15SJ-VS-CKSJ group (Fig. 7

The role of transcription factors in rice responses to cold stress
The key of plant tolerance to cold stress is a complex regulatory network involving TFs and other regulatory genes that control enzymes, regulatory proteins, and metabolites. In this study, genes from five TF families including the C2C2-Dof, C2C2-GATA, WRKY , GRAS The WRKY TFs are involved in the mitogen-activated protein kinase (MAPK) signaling pathway, which is involved in stress-induced defensive responses (Asai et al, 2002 impairs the cold-induced CBF expression and decreases freezing tolerance (Chinnusamy et al, 2003, Ding et al, 2015. In our experiment, ICE1(LOC_Os01g0928000) was up regulated in D15DN more than that in D15SJ.
Many WRKY proteins have been proven to respond to various abiotic stresses, such as cold , Hwang et al, 2005. Several CBF-independent regulators which were cold-induced transcription factors, function in a similar manner to CBFs, and can induce the expression of COR genes under cold stress, including WRKY33 andMYB73 (Park et al, 2015, Ding et al, 2019.In this study, we found that more WRKY and MYB genes were up regulated in CKDN vs. D15DN than in CKSJ vs. D15SJ.

Nitrogen metabolism responses to cold stress
Although the actual roles of Glu in plant osmo tolerance remain controversial, there is supporting evidence for a positive effect on enzyme and membrane integrity, osmotic adjustment, and free radical scavenging. Still, some authors have argued that Glu accumulation under stress is a product of, and not an adaptive response to stress.
Previous study has found that the N uptake of rice was severely affected by the low-T w treatments during reproductive growth (Jia et al, 2019) and that Glu content and glutamate dehydrogenase (GDH) activity were important traits that influenced the grain yield and spikelet sterility, respectively (Jia et al, 2015).In this study, we found that with the low-T w treatment, the Glu content in D15DN was significantly reduced, and the Pro content was slightly increased, compared to the control. In D15SJ, the trends of Glu and Pro content were opposite to those of D15DN. Compared with D15SJ, the Glu content in D15DN was 64.7% less, while the Pro was increased by 42.9%. These results could be attributed to several factors. As indicated by the key enzyme activities of nitrogen metabolism and amino acid content (Fig. 1 b, c), the increase of GOGAT activity and GS activity in D15SJ (69.0% and 59.3%) was greater than that of D15DN (43.8% and 29.6%), while the increase of GDH activity in D15SJ(58.0%) was less than that of D15DN (85.3%), compared with the control. The GDH activity was demonstrated to confer resistance to various stress conditions (Dubois et al, 2003, Robredo et al, 2011 and GDH catalyzes the reversible oxidative deamination of Glu, to supply 2-oxoglutarate and ammonium (Aubert et al, 2001). Secondly, Glu plays a central role in the amino acid metabolism of plants (Seifi et al, 2013), which may orchestrate crucial metabolic functions under abiotic stresses, such as GABA and Pro biosynthesis, which each perform key roles in plant defense mechanisms (Forde and Lea, 2007). The finding that a large amount of the Glu of D15DN was used in the biosynthesis of Pro, which leads to the Glu content in D15DN was less, while the Pro was more than that of D15DN.
Moreover, Compared to D15SJ, the 1 ornithine carbamoyl transferase gene involving in Lornithine synthesis, showed significant up-regulation in D15DN. Indicating that low-T w treatments may promote Glu metabolism, synthesis of ornithine, and Pro in D15DN.

Conclusion
Cold stress during reproductive growth, resulted in the genes and proteins related to the biosynthesis pathways of cold stress being significantly differentially expressed in both rice cultivars investigated (DN and SJ). The DN was more efficient in N metabolism, the MAPK signaling pathway, and the regulation of transcription factors. The present study confirmed the known cold stress-associated genes and identified a number of putative new cold-response genes. The study revealed that the translational regulation under cold stress plays an important role in the cold-tolerant DN. The low-T w treatment affects the N uptake and N metabolism in rice, and promotes Glu metabolism, synthesis of ornithine, and Pro in cold-sensitive SJ.

Plant material and growth conditions
Two japonica rice (Oryza sativa. L. subsp. japonica) cultivars, Dongnong428 (DN) and Songjing10 (SJ), which are currently used for local rice production in China, were utilized in this investigation. DN is a cold resistant cultivar, while SJ is a cold-sensitive cultivar.
A split-plot design was used in this experiment. The whole plot factor had two levels, normal irrigation (CK) and 15 days of low-T w irrigation (D15) during the reproductive growth stage of the rice. The cultivar of the rice (DN and SJ) was the subplot factor, the area of each subplot was 9 m 2 , and there were three replicates of the design.
Approximately 30-to 35-day-old seedlings were transplanted with one seeding per hill at a spacing of 13.3 cm (hill space) by 30.0 cm (row space).The application of the low-T w treatment was started at the beginning of the reproductive growth stage, when the young panicle was approximately 1 cm in length (da Cruz et al, 2006). The low water temperature for the treatment was 17°C with an approximately 20 cm water layer depth, 17°C is close to the temperature that causes cold damage to rice during reproductive growth in Heilongjiang Province, China(Jia et al, 2015). The cold water (at 17°C) was generated automatically using a temperature controlled cold water irrigation system. The low-T w irrigation was applied from 6:00 AM to 8:00 PM every day for the 15-day period.
The experimental approach used for measuring the gross root growth was the previously reported mesh bag method (Steen, 1991). For each root sample, according to the planned to 1600 m/z at a resolution of 70,000; MS/MS fragment scan range was >100 m/z at a resolution of 17,500 in HCD mode; normalized collision energy was 27%; dynamic exclusion time was 15 s; automatic gain control (AGC) for the full MS target and the MS2 target was respectively 3e6 and1e5; the number of MS/MS scans following one MS scan: the 20 most abundant precursor ions were above a threshold ion count of 20000.

Database search and quantification
Protein identification and quantification were performed with Proteome Discoverer Software (Thermo Fisher Scientific), using two included algorithms, Mascot and SEQUEST.
Searches were constructed against the genome database. Peptide and MS/MS tolerance was set to 10 ppm and 0.6 Da, respectively. Enzyme specificity was set to trypsin, with two missed cleavages. Each confident protein identification and quantification required at least one unique peptide. The false discovery rate (FDR) of the identified proteins was ≤0.01. For each of the three biological repeats, spectra were combined into one file and searched. A peptide confidence level of 95%, or an unused confidence score > 1.3, were used as the qualification criteria. The relative quantification of proteins was calculated based on the ratio of peak areas from the MS/MS spectra. Differentially abundant proteins were determined using a Student's t-test. Proteins with a < 0.05 P-value and a fold change ≥1.5 were considered up-regulated and a fold change ≤0.67 was considered to be downregulated.

Proteome bioinformatic analysis
The R package top GO version 2.28.0 was used for the GO enrichment analysis, with the classic algorithm option (each GO term was tested independently            c, Protein-protein interaction analysis of the specific CKSJ vs. D15SJdifferentially expressed genes, preformed using Cytoscape.

Figure 6
Identification of the transcription factors (TFs) responding to cold responses.