Comparative Transcriptome Analysis of Shoots and Roots of TNG67 and TCN1 Rice Seedlings under Cold Stress and Following Subsequent Recovery: Insights into Metabolic Pathways, Phytohormones, and Transcription Factors

Cold stress affects rice growth, quality and yield. The investigation of genome-wide gene expression is important for understanding cold stress tolerance in rice. We performed comparative transcriptome analysis of the shoots and roots of 2 rice seedlings (TNG67, cold-tolerant; and TCN1, cold-sensitive) in response to low temperatures and restoration of normal temperatures following cold exposure. TNG67 tolerated cold stress via rapid alterations in gene expression and the re-establishment of homeostasis, whereas the opposite was observed in TCN1, especially after subsequent recovery. Gene ontology and pathway analyses revealed that cold stress substantially regulated the expression of genes involved in protein metabolism, modification, translation, stress responses, and cell death. TNG67 takes advantage of energy-saving and recycling resources to more efficiently synthesize metabolites compared with TCN1 during adjustment to cold stress. During recovery, expression of OsRR4 type-A response regulators was upregulated in TNG67 shoots, whereas that of genes involved in oxidative stress, chemical stimuli and carbohydrate metabolic processes was downregulated in TCN1. Expression of genes related to protein metabolism, modification, folding and defense responses was upregulated in TNG67 but not in TCN1 roots. In addition, abscisic acid (ABA)-, polyamine-, auxin- and jasmonic acid (JA)-related genes were preferentially regulated in TNG67 shoots and roots and were closely associated with cold stress tolerance. The TFs AP2/ERF were predominantly expressed in the shoots and roots of both TNG67 and TCN1. The TNG67-preferred TFs which express in shoot or root, such as OsIAA23, SNAC2, OsWRKY1v2, 24, 53, 71, HMGB, OsbHLH and OsMyb, may be good candidates for cold stress tolerance-related genes in rice. Our findings highlight important alterations in the expression of cold-tolerant genes, metabolic pathways, and hormone-related and TF-encoding genes in TNG67 rice during cold stress and recovery. The cross-talk of hormones may play an essential role in the ability of rice plants to cope with cold stress.


Introduction
Rice is the most important staple food in the world, especially in Asia.Two subspecies of rice, Oryza sativa ssp.japonica (temperate rice) and ssp.indica (tropical rice), are widely grown in different areas.Rice seedlings frequently experience cold injury, which affects their growth and yield.In general, indica rice tends to be more sensitive to low temperatures.Thus, to further improve rice quality and production and to overcome the limiting factor of cold, a thorough understanding of cold stress tolerance mechanisms in rice is needed, especially the differential means of cold stress perception and responses to this type of stress in the indica (e.g., TCN1) and japonica (e.g., TNG67) rice varieties.
To adapt to environmental stresses, energy conservation and metabolic homeostasis are pivotal for all organisms.Under cold stress, various physiological and biochemical responses are altered in plants, such as the inhibition of photosynthesis, respiration and protein translation, accumulation of reactive oxygen species (ROS), alterations in metabolite profiles and osmolyte adjustment.Therefore, energy deprivation is likely a consequence of stress damage, which is often associated with reduced photosynthesis or respiration, ultimately resulting in growth arrest and cell death.Under abiotic stress, plants can reprogram or reconfigure their primary metabolism to redistribute energy resources for survival [1].Alterations in primary metabolism involving sugars and sugar alcohols, amino acids and tricarboxylic acid (TCA) cycle intermediates are general trends in abiotic stress responses [2].
Many growth and developmental processes in plants are affected by the balance and coordination of different plant hormones.Fluctuations in stress-responsive phytohormone levels are central to integrating stress signaling and regulating stress responses [3].However, the involvement of plant hormones in abiotic stress, especially cold stress tolerance, in rice remains poorly understood.To reveal the activities of plant hormones in rice seedlings at low temperatures, cold damage and levels of plant hormones, including abscisic acid (ABA), ethylene (ET) and polyamine, were evaluated in seedlings of 2 rice cultivars with contrasting responses to cold, TNG67 and TCN1 [4].The TNG67 seedlings were remarkably cold-tolerant compared with those of TCN1, as reflected by electrolyte leakage, the tetrazolium chloride reduction assay results and the survival ratio.After incubation at 5°C for 3 hr, the stomata of TNG67 immediately closed, but those of TCN1 did not, indicating the presence of wilt symptoms in TCN1 [5].In the cold-tolerant cultivar TNG67, ABA levels rapidly increased in roots and shoots in response to cold stress, but this did not occur in the cold-sensitive cultivar TCN1.Interestingly, exogenous addition of ABA enhanced cold stress resistance in TCN1 rice seedlings.The levels of both 1-amino-cyclopropane-1-carboxylic acid (ACC) and ET were decreased in TNG67 and TCN1 in response to cold treatment, and TCN1 was not able to restore the ET level after removal from the cold [4].The level of another plant hormone, putrescine (Put), and the activity of arginine decarboxylase (ADC) were enhanced in both the shoots and roots of TNG67.The levels of spermidine/spermine and activity of S-adenosylmethionine decarboxylase (SAMDC) were increased in TNG67 shoots.In TCN1, although the level of Put and activity of ADC were increased in shoots, the magnitudes of these increases were smaller than those observed in TNG67.Further, the Put level in roots was decreased after cold treatment [6].transcriptome studies have been performed in Arabidopsis thaliana and rice, a significant differences in organ-specific transcript profiling has been observed [20].Both of these types of plants clearly use distinct combinations of genes for abiotic stress responses and regulation.
In the present study, to gain insights into the key components that determine rice cold stress tolerance, we performed transcriptome analysis with a custom designed oligonucleotide array, a rice OneArray microarray platform (Phalanx Biotech Group Inc. [21]), using shoot and root tissues from TNG67 and TCN1 rice plants.Here, we discuss new findings from data analyses of DEGs, metabolic processes, changes in phytohormone-related gene expression, and important TFs.

Plant materials
Seeds of the rice cultivars TCN1 and TNG67 were obtained from the Taichung and Kaohsiung experimental research stations in Taiwan.Seeds were immersed in water for germination and then hydroponically cultivated in Kimura B solution at pH 5.0.Seedlings were grown in a phytotron at 30°C day/25°C night and 70~80% relative humidity.The plants were transferred to a growth chamber at 4°C for a cold treatment after growth to the three-leaf stage.After a 3-or 24-hr cold treatment or recovery to normal growth conditions for 24 hr, rice roots and shoots were harvested and immediately frozen in liquid nitrogen for RNA extraction.

Measurement of electrolytes and Fv/Fm
For electrolyte leakage assay, shoots and roots of five rice seedlings were washed and immersed in deionized water and then kept in the dark for 24 hr at 25°C.The amount of electrolytes that leached into the solution was measured with a conductivity meter (Extech EC500 Waterproof ExStik pH/Conductivity Meter).The change in conductivity was calculated as the ratio of the percentages of electrolyte leakage before and after autoclaving.Conductivity after autoclaving is assumed to represent complete (100%) electrolyte leakage.The maximum quantum efficiency of PSII photochemistry (Fv/Fm) was recorded as an indicator of PSII capacity.A JUNIOR-PAM chlorophyll fluorometer (Heinz Walz) was used to determine Fv/Fm at a position located 2 cm from the leaf tip of the third leaf of each rice seedling with or without cold treatment at various time points.

Microarray analysis
Transcriptome analysis was conducted using a 22-k Rice Whole Genome OneArray v1 oligo microarray platform (Phalanx Biotech Group, Inc., Taiwan), which was designed with high specificity and sensitivity for both japonica and indica rice.This array was constructed using MSU v6.1 database cDNA sequences, BGI indica (93-11) coding sequences and BGI japonica (Syngenta) coding sequences.Probes were designed and selected to cover 90% of the wellannotated genes found in both rice subspecies.The chip contains 22,003 probes that are randomly distributed across the array, each comprising a 60-mer oligonucleotide.In addition to 824 quality control probes, 21,179 probes covering 20,806 and 13,683 genes in japonica and indica were present on the chip, which has been well designed and validated [11].After combining protein-coding and non-protein-coding loci with rice FLcDNA, the exact number of genes shared by both cultivars was determined to be 13119.The gene array list with detailed descriptions can be found at http://www.phalanx.com.tw/products/RiOA_Probe.php.
Total RNA extractions from shoots and roots were performed using TRIzol reagent (Invitrogen) according to manufacturer's instructions.RNA quality was determined based on 260/ 280 ratio and a 260/230 ratio of > 2.0.One microgram of total RNA was reverse transcribed and amplified using an Amino Allyl aRNA Amplification Kit (Phalanx Biotech Group, Taiwan) and labeled with Cy5 dyes (Amersham Pharmacia, Piscataway, NJ, USA).Fluorescent targets were hybridized to the rice whole-genome OneArray with Phalanx hybridization buffer using a Phalanx Hybridization System.After 16 hr of hybridization at 50°C, non-specific binding targets were washed away via three different washing steps (wash I, 42°C for 5 min; Wash II, 42°C for 5 min and 25°C for 5 min; and Wash III, rinse 20 times).The arrays were scanned using a GenePix 4200A 01 Autoloader (635 nm, power of 100, and PMT of 520; and 532 nm, power of 10, and PMT of 460), and fluorescence intensity was quantified.The original dataset was deposited in the NCBI GEO database (accession no.GSE57895).
The signal intensity of each spot was loaded into a Rosetta Resolver System (Rosetta Biosoftware) for data analysis.The error model of this system allows for the removal of both systematic and random errors from data.We filtered out scanned spots with a flag value of < 0. Spots that passed the criteria were normalized by 50% media scaling normalization method.DEGs were considered those with changes in expression of greater than two-fold after normalization to the control and with significant results for the t test (p < 0.05) based on 6 replicates (3 biological repeats x 2 technical repeats) for each treatment compared with the control.These DEGs were verified by quantitative PCR (qPCR) and correlation analysis (R = 0.857; p < 0.05) (S1 Fig) .The heuristic cluster method was used for clustering analysis.

Quantitative RT-PCR
Integrated RNA was reverse transcribed into cDNA using SuperScript III Reverse Transcriptase (Invitrogen) for qPCR analysis.The synthesized cDNA was mixed with SYBR Green PCR Master Mix (Applied Biosystems), and real-time PCR was conducted with an Applied Biosystems ABI 7500 sequence detection system.The corresponding gene-specific primers used in gene expression analysis are listed in S1 Table.

Hierarchical cluster analysis
Hierarchical cluster analysis of the DEGs in the shoots and roots was performed with Gene Cluster 3.0 (http://bonsai.hgc.jp/~mdehoon/software/cluster/)followed by Alok Saldanha's Java TreeView v1.1.6r2for graphical analysis of the results.

MapMan analysis
The averaged signals for the cold treatment (3 biological replicates x 2 technical replicates) at different time points were expressed relative to those of the control samples, converted to a log 2 scale and displayed using MapMan v3.5.1 (http://gabi.rzpd.de/projects/MapMan/)[22].

Network analysis of co-expressed genes
Microarray network analysis was performed using Rice Co-expression Network (http://www.ricearray.org/coexpression/network.shtml).The network type selected was "Abiotic stress", the correlation coefficient cutoff was set at 0.8, and the depth to search co-expressed genes was set at 2. After online analysis, an SIF file was created and downloaded.We used Cytoscape v3.1.1 (http://www.cytoscape.org/) to visualize and edit the SIF file and generate graphs.

Results
Physiological and phenotype analyses showed that TNG67 is a coldtolerant rice variety TNG67 and TCN1 show contrasting cold tolerance at 4°C [4].To assess the effect of cold on photosynthetic capacity, we measured Fv/Fm values in leaves of TNG67 and TCN1 rice plants.After recovery from the cold for 24 hr, the Fv/Fm values were low in the TCN1 rice leaves but were similar to those measured under normal growth conditions in the TNG67 leaves (Fig 1A).To dissect the differences in the dynamic cold response between TNG67 and TCN1, we measured electrical conductivity and Fv/Fm parameters in both types of rice seedlings after various durations of cold treatment.The membrane leakage of TCN1 rice seedlings immediately increased following cold treatment for 3 hr and plateaued at 12 hr.However, the electrical conductivity of TNG67 rice seedlings remained relatively unchanged (Fig 1B).In response to cold treatment, the TCN1 but not the TNG67 leaves were pale, wilted and chlorotic (Fig 1C).Therefore, the TCN1 rice was more sensitive to low temperatures during the early phase (3 hr) and lost viability after 5 days of recovery.

Re-establishment of homeostasis during recovery is important, and shoots or roots may play a distinct role in regulating cold stress tolerance in rice
To compare the global gene expression changes in the rice cultivars under cold stress, we performed comparative transcriptome analysis with a customer-created oligonucleotide microarray and RNA from roots or shoots of TNG67 and TCN1 rice seedlings at early (3 hr; C3) and late (24 hr; C24) cold treatment and recovery time points (return to normal growth for 24 hr; Re24).DEG analysis was conducted with 3 individual biological and technical replicates, and the gene expression results revealed greater than two-fold changes in the expression of a number of genes after normalization to the control and significant t test results (p < 0.05) (S1 and S2 Tables).Comparisons of gene expression between the shoots of the two cultivars revealed a total of 834 DEGs in TNG67 and 579 DEGs in TCN1 at C3.Among them, 383 DEGs were shared by both cultivars.Following 24 hr of cold stress (C24), 2032 DEGs were identified in TNG67 in addition to 2168 DEGs in TCN1, and the two varieties shared 507 DEGs.After recovery, 1334 DEGs were present in TCN1, but only 290 DEGs were detected in TNG67, and only 148 DEGs were shared by both varieties (Fig 2A).In contrast, 1117 and 1269 DEGs were identified in the roots of TNG67 and TCN1, respectively, at C3 (705 DEGs in common), and 2652 and 3280 were detected at C24, respectively (1870 DEGs in common).After recovery, we discovered 1842 DEGs in TNG67 and 3096 in TCN1, 1029 of which were shared by both cultivars (Fig 2B).In response to cold treatment, the number of DEGs identified in roots was greater than those in shoots in both TNG67 and TCN1 (Fig 2B).This result indicated that root transcript levels were more easily regulated compared with those in shoots.Thus, roots may be more responsive to cold stress than shoots in both cultivars.Additionally, we found markedly increased numbers of DEGs in the shoots (by approximately five-fold) and roots (by approximately two-fold) of TCN1 during recovery.The expression patterns of the DEGs in the shoots and roots suggest that these 2 different tissues may play distinct roles in regulating rice cold stress tolerance.
Based on the similarity of the gene expression patterns between the roots and shoots of the two cultivars, we analyzed the microarray data using a hierarchical clustering method.After the 3-and 24-hr cold treatments and subsequent recovery of TNG67 and TCN1 shoots and roots, the DEGs could be separated into 3 distinct clusters (Fig 3).Cluster I contained those that were repressed under cold stress but induced after recovery.Cluster II comprised shoot but not root DEGs that were induced in response to cold stress.Cluster III included DEGs that were induced by cold in both shoots and roots.Cold treatment and recovery distinctly affected gene expression at various time points.
Pathways that mediate protein metabolism, modification, translation, stress response, and cell death are important for cold stress tolerance in rice To group DEGs with similar functions and identify functional categories with different DEGs in rice in response to cold stress, we performed gene ontology (GO) analysis of the microarray dataset.The numbers of DEGs in almost every GO Slim category of biological process were higher in TNG67 than in TCN1 in response to cold stress.
In shoots (Table 1), some of the GO Slim categories with upregulated genes were significantly enriched in TNG67 but not in TCN1 during the early stage (3 hr) of cold stress, including protein metabolic process (GO:0019538), protein modification process (GO:0006464), response to stress (GO:0006950), and cell death (GO:0016265).GO:0006091 (generation of precursor metabolites and energy) was an exception in TCN1, in which it was significantly enriched compared with TNG67.During the late stage of cold stress (C24), the numbers of upregulated genes in each GO Slim category were higher in TNG67 than in TCN1.However, the numbers of downregulated genes in each GO category were lower in TNG67 than in TCN1.After recovery, the DEGs in many GO Slim categories remained significantly enriched in TCN1.During recovery, only GO:0000160 [two-component signal transduction system (phosphorelay)] was significantly enriched among the upregulated genes in TNG67 and may be important for maintaining homeostasis during recovery (S4 Table ).GO:0006979 (response to oxidative stress), GO:0042221 (response to chemical stimulus) and GO:0005975 Table 1.Gene ontology (GO) annotation of differentially expressed genes (DEGs) assigned to the biological process (BP) category after cold treatment and recovery in shoots.The bars show the numbers of genes in each GO Slim category.GO enrichment analysis was conducted using agriGO (http://bioinfo.cau.edu.cn/agriGO/), and GO terms with a p < 0.05 were considered to have significantly enriched expression in the cluster.

Shoot TNG67 TCN1
Up-regulated (carbohydrate metabolic process) were significantly enriched among the downregulated genes in TCN1.Malfunctions in these processes during recovery may lead to death following cold stress.
In roots (Table 2), GO:0019538 (protein metabolic process), GO:0006464 (protein modification process), and GO:0006412 (translation) were significantly enriched at both the DNA and protein levels among the upregulated genes in TNG67 but were not or were less enriched in TCN1 after 3 hr of cold stress.After removal from the cold for 24 hr, the DEGs were significantly enriched in the GO:0006457 (protein folding), GO:0008643 (carbohydrate transport) and GO:0006952 (defense response) categories in TNG67 but not in TCN1 (S5 Table ).

Expression of genes involved in TCA cycle, glycolysis, proteasome, phenylpropanoid and lignin biosynthesis pathways was significantly induced in TNG67 rice in response to cold stress
To further investigate alterations in metabolic pathways that occur during cold stress and recovery, we performed MapMan analysis to visualize the DEG data from microarray analysis.MapMan has been applied with Rice Net tools for systematic analysis of the early heat stress transcriptome in rice seedlings [23].For shoots, the metabolic overview map (Fig 4) revealed a greater enhancement of gene expression associated with light reactions and mitochondrial electron transport (S6 Table) in TCN1 than in TNG67.In contrast, the expression of genes related to the TCA cycle and glycolysis pathways was elevated in TNG67 compared with TCN1.The detailed genes ID/description identified aforementioned are provided in EXCEL files as supplementary S6  Table 2. GO annotation of DEGs assigned to the biological process (BP) category after cold treatment and recovery in roots.The bars show the numbers of genes in each GO Slim category.GO enrichment analysis was conducted using agriGO (http://bioinfo.cau.edu.cn/agriGO/), and GO terms with a p < 0.05 were considered to have significantly enriched expression in the cluster.

Root TNG67 TCN1
Up-regulated  compared with TCN1 following 3 hr of cold treatment.In particular, E2 ligase gene expression was predominant in TNG67 following 24 hr of treatment (S6 Table ).
In roots, the secondary metabolism map generated by MapMan analysis showed a tendency toward the increased activities of the phenylpropanoid and lignin biosynthesis (S6 Table ) pathways in TNG67.Expression of the genes in these two pathways was greater in TNG67 than in TCN1 after 3 hr of cold stress.Additional genes were upregulated after 24 compared with those upregulated after 3 hr of cold stress in both cultivars; however, more genes were repressed in TCN1.In contrast, more genes were repressed in TNG67 after recovery compared with those repressed after 24 hr of cold stress (Fig 5).
ABA, polyamines, auxin and jasmonic acid (JA)-related genes were preferentially expressed in TNG67 shoots and roots and were closely associated with cold stress tolerance Plant hormone cross-talk occurs under various abiotic stresses and substantially regulates plant growth and development in response to stress [3,24].However, fluctuations in the levels of stress-responsive plant hormones occur during different stages, including biosynthesis, degradation, transportation, storage, perception and signal transduction.To obtain a better understanding of the roles of different plant hormones in rice in response to cold stress, we examined the expression profiles of various plant hormone-related genes in root and/or shoot tissues of rice cultivars using our microarray dataset.Expression profiles of ABA-related genes ABA is a major stress hormone involved in plant responses to abiotic stresses.Previously, the ABA concentration has been found to rapidly increase in the shoots and roots of TNG67 but not in those of TCN1 in response to cold stress [5].However, the molecular mechanism underlying the involvement of ABA in the cold stress response has not been well characterized.
Based on our microarray analysis (Fig 6), the expression of LOC_Os09g38320 (PSY3) and LOC_Os01g12710 (SDR) was induced by cold temperatures to a greater extent in the shoots of TNG67 compared with the levels observed in those of TCN1.PSY3 expression induces ABA biosynthesis under abiotic stresses [25].These genes may contribute to the rapid biosynthesis and response of ABA in TNG67 shoots under cold stress.In TNG67, LOC_Os01g03750 (ABA4) and LOC_Os12g42280 (NCED2) transcripts were highly induced at low temperatures in both shoots and roots compared with TCN1.In addition, the expression levels of LOC_Os06g40170 (PLDα6), LOC_Os06g40180 (PLDα7), LOC_Os06g40190 (PLDα8), LOC_Os10g38060 (PLDα13), and LOC_Os12g39630 (SAPK9), which are involved in the ABA signal transduction pathway, were higher in the shoots of TNG67 than in those of TCN1.Two ABA-glucosyltransferase (ABA-GTase) genes, LOC_Os07g32010 and LOC_Os07g32060, which convert the storage form of ABA into its free form, were downregulated in TNG67 shoots.Furthermore, the expression levels of LOC_Os02g47470 (ABA8OX1) and LOC_Os09g28390 (ABA8OX3), which encode different isoforms of ABA 8'-hydroxylase genes responsible for the deactivation of ABA, were higher in shoots compared with roots.Therefore, the transcriptional regulation of ABA-related gene expression is one factor that controls cold stress tolerance in TNG67.

Expression profiles of ET-related genes
Lee et al. [4] have shown that ACC and ET levels are decreased in TNG67 and TCN1 in response to cold treatment, and the ET level is increased after recovery in TNG67, but not in TCN1.Our microarray data (Fig 6 ) showed that most ET biosynthesis and signaling genes were upregulated in response to cold.ET biosynthesis-related genes, including LOC_Os05g04510 (SAMS1), LOC_Os04g48850 (ACS2) and LOC_Os05g05680 (ACO5) and the signaling genes LOC_Os02g57530 (ETR2; 1) and LOC_Os07g06130 (EIN2; 1), were significantly upregulated in both cultivars and tissues.Some ET biosynthesis genes, such as LOC_Os09g27820 (ACO1), LOC_Os05g05680 (ACO5) and LOC_Os05g05670 (ACO6), were upregulated in TNG67 roots at low temperatures and displayed prolonged expression during recovery.Therefore, the cold stress-mediated ET response may not be regulated at the transcriptional level in rice.

Expression profiles of polyamine-related genes
Polyamines are abiotic stress modulators that interact with ABA and ET to affect plant abiotic stress tolerance [26].Put is specifically involved in freezing tolerance and cold acclimation in Arabidopsis by regulating ABA levels in response to cold stress [27].Thus, we assessed the expression profiles of polyamine-related genes.In response to cold stress, the levels of Put and ADC activity were increased in both the shoots and roots of TNG67.Furthermore, spermidine/spermine and SAMDC levels were increased in the shoots of TNG67 and were increased to a lesser extent in the shoots of TCN1 [6].The expression of ADC genes, including LOC_Os06g04070 (ADC) and LOC_Os08g33620|LOC_Os06g04070 (ADC), was highly induced in TNG67 compared with TCN1 in both shoots and roots.The expression pattern of SAMDC genes, including LOC_Os02g39790 and LOC_Os04g42090, was similar to that of ADC genes.In addition, the expression of polyamine transporter genes in shoots, including LOC_Os01g41420, LOC_Os04g47780 and LOC_Os12g08130, was higher in TNG67 than in TCN1, particularly after 24 hr of cold stress (Fig 6).Therefore, the ADC-mediated pathway of Put biosynthesis appears to be an important mechanism of cold tolerance in TNG67.

Expression profiles of gibberellin-related genes
Reductions in the endogenous gibberellic acid (GA) level [28] and DELLA activity [29] enhance the survival of Arabidopsis during salt stress.However, whether GA is linked to cold stress tolerance in rice is unknown.We examined the expression of GA-related genes and found the suppressed expression of GA biosynthesis-related genes, including LOC_Os06g02019 (KAO), LOC_Os03g63970 (GA20ox1), and LOC_Os01g66100 (GA20ox2), under cold conditions; in particular, LOC_Os01g66100 (GA20ox2) expression was more markedly downregulated in TNG67 than in TCN1 in both shoots and roots.In contrast, the decreased expression of genes in shoots, including LOC_Os01g55240 (GA2ox3) and LOC_Os04g44150 (GA2ox6), was upregulated in TNG67 compared with TCN1 (Fig 6).Thus, the reduction in GA concentration via transcriptional gene regulation may enhance cold stress tolerance in TNG67.

Expression profiles of cytokinin (CK)-related genes
Recent data have suggested that the cross-talk between ABA and CK, and particularly the ABA:CK ratio in xylem sap, are important for stress signaling [30].In addition, the CK twocomponent system is involved in the stress response to low temperatures in Arabidopsis [31].Thus, we evaluated the CK-related gene expression profile in our array dataset.In the CK biosynthetic pathway, we identified only one gene, LOC_Os01g40630 (CK riboside 5'-monophosphate phosphoribohydrolase, LOG), that was induced in TNG67 shoots after 24 hr of cold stress.With regard to CK degradation, we discovered 2 catabolic enzymes, the CKX (CK dehydrogenase)-related genes LOC_Os10g34230 (CKX3) and LOC_Os01g71310 (CKX4), that were expressed at higher levels in the shoots of TNG67 compared with those of TCN1 in response to cold stress.However, LOC_Os10g34230 (CKX3) was preferentially induced in TCN1 roots.The expression of the majority of the genes associated with the CK signal transduction pathway was preferentially repressed in TNG67 compared with TCN1, including LOC_Os04g36070 (RR1) and LOC_Os07g26720 (RR7) in shoots and LOC_Os03g50860 (HK4), LOC_Os02g58350 (RR3), LOC_Os01g72330 (RR4) and LOC_Os07g26720 (RR7) in roots (Fig 6) . CK may participate in antagonistic cross-talk with ABA to negatively regulate cold stress tolerance in rice.During the recovery stage, the CK synthesis gene LOC_Os05g47840 (IPT7) was upregulated in the roots of TNG67 but not in those of TCN1; however, this gene was downregulated in response to cold stress in both cultivars.In shoots, the gene expression profile of CK signal transduction, including the OsRR type-A genes LOC_Os04g36070 (RR1), LOC_Os02g35180 (RR2), LOC_Os01g72330 (RR4), LOC_Os04g57720 (RR6) and LOC_Os11g04720 (RR10), showed a similar pattern.

Expression profiles of auxin-related genes
Auxin is tightly linked to cold stress-induced changes in plant growth and development [32].Endogenous levels of indoleacetic acid (IAA) and JA are differentially modulated in rice under abiotic stress [33].We analyzed the expression of genes related to auxin biosynthesis and found that LOC_Os04g38950 (ASB1) and LOC_Os04g03980 (YUCCA7) were induced under cold treatment in both tissues of the rice cultivars.LOC_Os07g25540 (YUCCA6) was induced in both tissues of TCN1 but only in the roots of TNG67.In contrast, LOC_Os10g04860 (AAO1) was upregulated in TNG67 shoots and roots but was increased only in TCN1 roots.Increased auxin expression under low temperatures may help to ameliorate cold injury in rice.Furthermore, the expression of the auxin signaling-related genes LOC_Os01g08320 (IAA1), LOC_Os01g53880 (IAA6), LOC_Os05g48590 (IAA19), and LOC_Os06g47150 (ARF18) was induced in both cultivars and tissues.The expression of most of these genes was higher in TNG67 than in TCN1.Other genes, such as LOC_Os05g05800 (TIR1; 1), were induced in TNG67 shoots and roots [LOC_Os04g59430 (ARF13) and LOC_Os01g18360 (IAA4)] (Fig 6).Thus, auxin expression may be higher in the shoots and roots of TNG67 than in those of TCN1 under cold stress.Interestingly, GUS staining was increased in the roots of DR5-GUS transgenic rice seedlings treated with cold stress (S3 Fig).

Expression profiles of jasmonate-related genes
JA enhances the freezing tolerance of Arabidopsis by positively regulating the ICE-CBF/DREB1 pathway [34].We dissected the JA-related gene expression profile in response to cold stress in rice.The levels of JA biosynthesis genes, including LOC_Os08g04800 (DAD1; 2), LOC_Os12g37260 (LOX2; 2) and LOC_Os08g35740 (OPR7), were higher in TNG67 than in TCN1 in shoots, and that of LOC_Os03g12500 (AOS1) was increased in roots.LOC_Os12g37260 (LOX2; 2) was the only gene that was specifically expressed in TNG67 shoots.These findings may indicate a pivotal role of JA synthesis and the contribution of increased JA content in triggering cross-talk between JA and cold stress.Our results are consistent with a previous report [34] demonstrating higher JA signaling-related expression [LOC_Os03g08320 (JAZ1), LOC_Os10g25290 (JAZ2), LOC_Os07g42370 (JAZ3), LOC_Os03g08310 (JAZ6), LOC_Os04g32480 (JAZ9), and LOC_Os10g25230 (ZIM18)] in TNG67 than in TCN1 (Fig 6).The addition of 50 μM IBU (JA biosynthesis inhibitor) to hydroponically cultured TNG67 and TCN1 rice seedlings in response to cold stress for 24 hr induced the hypersensitivity to cold stress of both cultivars (S4 Fig) .NAC, WRKY and ERF/AP2 TFs are good candidates for involvement in cold stress tolerance The expression of plant TFs is typically quickly altered in response to environmental stimuli, leading to the establishment of a pleiotropic phenotype to enhance stress tolerance.To elucidate the transcriptional regulation of cold stress-responsive genes, we identified DEGs encoding TFs under cold stress.For Venn diagram analysis, we determined the gene expression patterns of DEGs encoding TFs according to the Database of Rice Transcription Factors (http://drtf.cbi.pku.edu.cn/)(Fig 7).Following 3 hr of cold treatment, TNG67 and TCN1 shared 71 up-and 27 downregulated TF-encoding DEGs in shoots and 61 up-and 63 downregulated TF-encoding DEGs in roots.While 51 upregulated and 26 downregulated TFs were uniquely expressed in the shoots of TNG67, only 10 each of upregulated and downregulated TFs were uniquely expressed in those of TCN1.However, after 24 hr of cold stress, 30 upregulated and 30 downregulated TFs were uniquely expressed in TNG67 roots, and 32 upregulated and 73 downregulated TFs were uniquely expressed in TCN1 roots.Sixty-one up-and 63 down-regulated DEGs were shared by both cultivars.After 24 hr of cold treatment, greater numbers of TF-encoding DEGs (260 and 270 in the corresponding shoots and roots) were observed in these 2 cultivars.TNG67 and TCN1 shared fewer TF-encoding DEGs during the recovery stage.Twenty-three and 123 common DEGs were shared in the shoots and roots, respectively, of TNG67 and TCN1 (Fig 7).The detailed genes list of aforementioned TFs can be found in supplementary tables (S7 and S8 Tables).Nevertheless, TCN1 possessed more TFencoding DEGs during the recovery stage.TNG67 only showed 6 up-and 12 downregulated TFs in shoots compared with TCN1.TNG67 roots displayed the specific expression of 46 upand 24 downregulated TFs.Only one NF-YB TF (LOC_Os03g29970) was upregulated in both the shoots and roots of TNG67.In addition, we stringently selected TF-encoding DEGs that were specifically or preferentially expressed in TNG67 (induced or repressed by greater than two-fold compared with TCN1) in response to the cold treatment and subsequent recovery (Fig 7).We identified 45 and 30 TFs in shoots and roots (S9 and S10 Tables) during the early stages of cold treatment (3 hr).The major categories were NAC-type TFs in shoots and WRKY-type TFs in roots.AP2/ERF-type TFs were uniquely expressed in shoots and roots and specific TFs were not shared between the two tissue types.In response to cold stress, AP2/ERFtype TFs represented the greatest number of DEGs in the shoots and roots of both cultivars.During recovery, bHLH-type TFs represented the major DEGs in the roots of both cultivars, although the expression of these genes was downregulated (Fig 8).The MYB and MYB-related TFs were the major upregulated TFs following cold recovery in the roots of TNG67.Several TFs were differentially expressed only during recovery in TNG67.In shoots, a GRAS TF (LOC_Os05g42130) was upregulated in TNG67 but downregulated in TCN1.bHLH (LOC_Os01g38610) and IAA20 (LOC_Os06g07040) were downregulated under cold stress but upregulated after recovery.In contrast, two AP2/ERF (LOC_Os01g21120 and LOC_Os02g52670) TFs and OsWRKY15 (LOC_Os01g46800) were upregulated in response to cold stress but downregulated after recovery.In roots, a TNG67-specific MYB TF (LOC_Os02g45080) was upregulated after recovery.In addition, AP2/ERF (LOC_Os01g64790) and HMG (LOC_Os02g44930) TFs were upregulated in response to cold stress but downregulated after recovery.

Co-expression gene regulatory network
Because the key DEGs detected, including 45 and 30 TFs in TNG67 shoots and roots, may represent crucial genes involved in cold stress tolerance in rice, we performed gene co-expression network analyses.We built a co-expression network with these TFs based on Rice Oligonucleotide Array Database (http://www.ricearray.org/coexpression/coexpression.shtml).Among the query genes, co-expressed genes related to abiotic stress with a correlation coefficient (r) cutoff of 0.8 were selected.A total of 45 TFs in shoots were inputted as query genes, generating 5 modules representing some TFs that were significantly linked to abiotic stress.An HMGbox DNA-binding protein (HMGB protein) (LOC_Os08g01100) and a CO-like protein that is orthologous to barley HvCMF11 (LOC_Os10g41100) were found to be highly correlated with abiotic stress-related genes, and the PHD ( TFs that are known to be induced by cold, drought and submergence, including CO-like protein, OsIAA23, bHLH, SNAC2 and other NAC genes, were preferentially expressed in the shoots of TNG 67 but not in those of TCN1.In contrast, TFs such as SNAC2, MYB-containing DNA binding protein, OsWRKY1v2, OsWRKY53, OsWRKY71 and OSWRKY24, were identified in the roots of TNG67 but not in those of TCN1.A literature search revealed that SNAC2 is induced by drought, salinity, cold, wounding, and ABA treatment.Most importantly, the overexpression of SNAC2 in Zhonghua 11 2 (japonica) enhances its survival after severe cold stress (4-8°C for 5 days) [11].These findings undoubtedly explain why the expression of this TF (SNAC2) was highly induced specifically in the shoots and roots of cold-tolerant TNG67 but not in those of TCN1.Moreover, OsIAA23-mediated auxin signaling has been shown to be crucial for the maintenance of the quiescent center in root tips (Jun et al., 2011).Additionally and interestingly, OsWRKY24, 51, 71, and 72 are induced by ABA and are involved in the ABA/GA balance in aleurone cells [35].
In the present study, we detected OsWRKY24 and 71 preferentially induced genes in TNG67 roots in response to cold stress, providing additional information regarding the role of ABA in rice cold stress tolerance, as has been shown in a previous physiological study.The other two WRKY TFs (OsWRKY 1v2 and OsWRKY53) were shown to be involved in salicylic acid and JA-mediated disease defense.OsWRKY 1v2 has ever been shown to be expressed at the early phase of chilling stress in japonica rice, while WRKY53 is expressed at the late stage (phase 2) which has been associated with disease resistance [36].The expression patterns of these WRKY TFs may represent cross-talk occurring between abiotic and biotic stress responses.These TFs in the co-expression network may represent dominant roles in the regulation of the cold stress tolerance of TNG67.

Discussion
Rice, particularly indica subspecies, is susceptible to chilling stress.Low temperatures can affect rice development during the germination, vegetative growth and reproductive stages.The screening of cold-tolerant rice varieties tends to be unsuccessful because of the poor correlation of cold resistance with different developmental periods [37].In addition, cold stress tolerance is a quantitative trait that is determined by various quantitative trait loci.To understand the cold stress tolerance mechanisms used by rice to cope with cold sensitivity and to develop coldtolerant rice cultivars, we conducted a comparative transcriptome study to identify changes in gene expression in response to cold and subsequent recovery between the roots and shoots of japonica (TNG67) and indica (TCN1) rice varieties.We aimed to discover some tolerancerelated genes, important TFs and pathways that may benefit the breeding of potential rice varieties that are capable of withstanding cold.

Differential expression in response to cold and during late recovery phases in TNG67 and TCN1
Comparative transcriptome analysis revealed more DEGs during the early cold response (3 hr, 4°C) in TNG67 than in TCN1 in both shoots and roots, with fewer DEGs in TNG67 than in TCN1 during the subsequent recovery phase (Fig 2).These results are consistent with those reported by Zhang et al. [14], who have shown that despite an intrinsic difference in constitutive cold tolerance-related gene expression in cold-resistant rice, the rapid and efficient reversion of gene expression during recovery remains an important factor in determining rice cold stress tolerance.Our results also indicate that TNG67 can reconfigure its genome for rapid gene expression and the re-establishment of homeostasis via metabolite adjustment when relieved of stress.Our gene clustering analysis clustered the DEGs into 3 distinct categories for the cold response and recovery phases (Fig 3).Genes that were specifically expressed in TNG67 during the early cold stress and recovery stages may represent cold tolerance-related genes that warrant further functional study.
Efficient energy, nutrient resource conservation and recycling, water use and cell wall remodeling are essential for rice cold stress tolerance Plants are sessile and unable to escape from environmental stresses.Thus, to acclimate, they must develop sensitive mechanisms for perceiving fluctuations in environmental cues and maintain homeostasis and the flexibility to adapt to environmental changes.Many stress conditions affect cellular ATP generation, protein synthesis and gene expression.The ability to ensure nutrient remobilization during energy deficiency by triggering the correct downstream transcriptional responses and metabolic adjustments is pivotal for plant survival.TCN1 rice was found to be sensitive to low temperatures, as reflected by a reduced photosynthetic capacity and increased membrane leakage (Fig 1).Although GO:0015979 (photosynthesis) was significantly enriched in both TNG67 and TCN1, the number of genes assigned to this GO term in TCN1 was greater than that in TNG67.MapMan analysis also revealed enhanced gene expression associated with light reactions in TCN1 (Fig 4).In addition, the increased expression of mitochondrial electron transport-related genes was observed in TCN1 compared with TNG67.Under cold stress, insufficient activities of the light reaction in photosynthesis and the electron transport chain deprives the cell of ATP energy, which resulted in the accumulation of more ROS in TCN1 than in TNG67.However, MapMan biological pathway analysis showed that TNG67 increased the expression of glycolysis-and TCA cycle-related genes to avoid energy starvation during cold stress.Gene set enrichment analysis (GSEA) of roots showed that carbohydrate transport (GO:0008643) was enhanced under cold stress in TNG67 but not in TCN1.TNG67 may be able to change carbohydrate-coordinated partitioning from roots to shoots so that the shortage in the carbon supply due to diminished photosynthesis can be compensated for under cold stress.These results suggest that TNG67 possesses a more efficient mechanism to conserve energy and overcome the limited photosynthetic capacity under cold stress.
We also identified several over-represented GO terms, including GO:0019538 (protein metabolic process), GO:0006464 (protein modification process), GO:0016265 (death), GO:0006950 (response to stress) and GO:0016265 (cell death), in TNG67 during the early stage of cold stress.GO:0016265 (death) was significantly repressed in TCN1 in response to 24 hr of cold stress.Translational regulation at the protein level and apoptosis or PCD may be important for rice cold stress tolerance, especially the protein modification process (GO:0006464).Previous reports have demonstrated a close association between PCD-related genes and increased cold stress tolerance in japonica rice with introgression of an indica allele [13].Cell death may be considered as not only simply a damaging consequence of cold stress but also part of a regulated process that permits flexibility in the cold stress response.Cell death is also an important physiological process for the removal of senescent and damaged cells during abiotic stress [38].Apoptosis (GO:0006915) may act as a signal for TNG67 under cold stress and trigger another set of genes associated with GO:0006950 (response to stress) to initiate cold tolerance mechanisms.In addition, the regulation of apoptosis is highly associated with the ubiquitin conjugation system [39].We found that in response to cold treatment, the expression of genes related to E2, E3 ligase, and the 20S proteasome were highly induced in TNG67 but not in TCN1 (S2 Fig) .Under stress, protein translation can be arrested, and protein levels can be regulated via a proteasome-dependent degradation process [40].Apoptosis and ubiquitin may play roles in metabolic reprogramming by recycling or reutilizing degraded substrates for the de novo synthesis of proteins to increase cold stress tolerance.
We identified three genes assigned to GO:0006950 (response to stress), LOC_Os04g16450, LOC_Os05g14240 and LOC_Os06g22960, that were predicted to be similar to aquaporin proteins and were induced in response to cold stress in both TNG67 and TCN1, but predominantly in TNG67 roots.Previous studies have demonstrated that under cold stress, the ability to control water absorption in roots is a more important characteristic to overcome stressrelated damage than the regulation of transpiration in the leaf [41].Our MapMan analysis revealed a tendency toward the increased synthesis of phenylpropanoid and lignin in the roots of TNG67 compared with those of TCN1 under cold stress.Cell wall remodeling was found to be a long-term adaption mechanism for plants to cope with multiple biotic or abiotic stresses [42].Whether this is also true for cold stress tolerance in TNG67 remains to be determined in future studies.
ABA, polyamine, auxin and JA hormone pathways act together or alone to regulate cold-responsive TFs and cold tolerance-related gene expression in rice Recent findings have demonstrated the importance of phytohormone cross-talk in the responses and resistance of plants to multiple abiotic stresses [24].In particular, the integration or coordination of hormone-related gene expression under different abiotic stress conditions has become essential to the understanding of plant abiotic stress tolerance mechanisms.Despite discrepancies in reported levels of regulation, from mRNAs to proteins and even metabolites, global transcriptome analysis of the expression, pathways and modules of coexpressed genes still represents a powerful approach for dissecting cross-talk among plant hormones during environmental stress.
In previous studies, Lee et al. [5] have demonstrated an accelerated increase in the concentration of ABA in TNG67 within 4 hr of cold treatment.The expression of ABA biosynthetic genes, including LOC_Os09g38320 (PSY3), LOC_Os01g1271tabolits0 (SDR), LOC_Os12g42280 (NCED2) and LOC_Os01g03750 (ABA4), was higher in the shoots of TNG67 compared with those of TCN1 and may explain why TNG67 had a higher ABA level during cold stress.Among these genes, LOC_Os12g42280 (NCED2) responded more strongly to cold stress in the roots of TNG67 than in those of TCN1 after 3 hr of exposure, and this effect was enhanced at 24 hr.This gene is a key contributor to ABA biosynthesis in TNG67 and TCN1 roots.The roots may represent the major location of the synthesis of ABA, which is transported to shoots during cold stress.When ABA is applied exogenously to plants, an immediate increase in its levels is observed in leaves [43].However, ABA levels do not increase significantly until leaf turgor has dropped to zero, despite the presence of drying soil.Under such conditions, plants that exhibit stomatal closure possess an efficient system for the transport of ABA from root to guard cells through the xylem [44].Furthermore, genes related to ABA signal transduction pathways were preferentially expressed in TNG67 shoots.During cold stress, the expression levels of the 2 ABA-GTase genes LOC_Os07g32010 and LOC_Os07g32060 were significantly decreased in TNG67 shoots; therefore, the increased ABA level may have been a result of the inactivation of ABA-GTase to switch ABA from a storage to a free form.Furthermore, the expression of LOC_Os02g47470 (ABA8OX1) and LOC_Os09g28390 (ABA8OX3), which are responsible for the deactivation of ABA, was induced to a greater extent in shoots compared with roots.In maize, the ABA catabolism rate increases by 11-fold in response to water stress.Pharmacological experiments have shown that this phenomenon may be beneficial for de novo ABA synthesis in roots and critical for rapid early adaptation when there is a requirement for ABA.These studies have revealed an important role of ABA in cold stress tolerance in TNG67.The role of roots in ABA-mediated cold tolerance and the mechanism(s) underlying the long distance root-to-shoot transport of ABA require further clarification.
The gene expression pattern revealed by our microarray analysis is inconsistent with that of a previous study investigating ACC and ET production in response to cold stress [4].A recent study has reported that ET is involved in the transcriptional regulation of CBF regulons and negatively affects freezing stress tolerance in Arabidopsis [45].This observation may aid in explaining our findings.As an alternative explanation, there may have been competition between polyamines and ET for SAM as a common substrate to produce methylthioadenosine.SAM activity shifts toward the polyamine biosynthesis pathway to promote the accumulation of polyamines under cold stress.These genes must be synthesized in advance during cold stress so that the corresponding gene products are quickly translated after removal of plants from the cold.In TNG67, the ACC and ET levels increased during recovery.However, in TCN1, even the ACC level increased during recovery, and it could not be converted into ET.The gradual increase in the ET level in response to cold stress, which reverted after recovery, may represent an important index for evaluating the capacity of cold stress tolerance in rice [4].The microarray data further showed that ET biosynthesis genes, such as LOC_Os09g27820 (ACO1), LOC_Os05g05680 (ACO5) and LOC_Os05g05670 (ACO6), which are responsible for converting ACC to ET, were upregulated in TNG67 roots at low temperatures and displayed prolonged expression during recovery.Ethylene exerts an antagonistic function against ABA, which further supports its role in detecting the disappearance of stress in rice plants.
Under cold stress, ornithine decarboxylase (ODC) genes were repressed, but ADC genes were induced.These results indicated that the Put biosynthesis pathways were shifted in the direction of ADC.In Arabidopsis, polyamine biosynthesis requires the activity of the rate-limiting enzyme ADC to synthesize Put [46], a key molecule that contributes to cold tolerance in Arabidopsis [47].ADC genes, including LOC_Os06g04070 (ADC) and LOC_Os08g33620| LOC_Os06g04070 (ADC), were preferentially induced in TNG67 compared with TCN1 in both shoots and roots.The transcriptional regulation of ADC corresponds to its activity and the Put level, which increased rapidly in TNG67; however, ODC activity remained unchanged after the cold treatment [6].The promoter sequence of OsADC1 was found to contain a CRT/ DRE element, which may mediate early and transient increases in OsADC1 expression in response to cold stress.In addition, the OsADC2 promoter contains 5 ABA-responsive elements that can increase OsADC2 expression after the de novo biosynthesis of ABA [48].At low temperatures, Put interferes with ABA biosynthesis and gene expression and affects the ABA level [27].The free spermidine: spermine levels and the activity of SAMDC are increased in chilled TNG67 [6].We found that SAMDC-encoding genes, including LOC_Os02g39790 and LOC_Os04g42090, were induced in TNG67 compared with TCN1, providing more aminopropyl donors for spermidine synthase (SPDS) or spermine synthase (SPMS), resulting in higher levels of free spermidine: spermine.Cross-talk among polyamines, ABA and ET coordinately establish the state of rice cold stress tolerance.
GA biosynthesis-related genes, including LOC_Os06g02019 (KAO), LOC_Os03g63970 (GA20ox1), and LOC_Os01g66100 (GA20ox2), were suppressed under cold, and LOC_Os01g66100 (GA20ox2) expression was higher in TNG67 than in TCN1 in both shoots and roots.In contrast, deactivated genes in shoots, including LOC_Os01g55240 (GA2ox3) and LOC_Os04g44150 (GA2ox6), were upregulated to a greater extent in TNG67 than in TCN1.Overall, GA seems to decrease the ABA level during cold stress, which is consistent with the antagonistic activity that has been reported between their levels in response to cold temperatures.In bromegrass, enhanced cold tolerance due to exogenous ABA is eliminated by exogenous GA [49].GA exerts its physiological function by degrading its receptor, DELLA.Under salinity and water stress, the endogenous accumulation of ABA and ET stimulates DELLA gene expression and leads to DELLA-dependent growth inhibition [24].This wellillustrated example integrating ABA, ET and GA at the level of DELLA may provide insights regarding rice cold stress tolerance.
Predicting changes in the CK level in response to cold stress is difficult because the expression of the genes responsible for its biosynthesis and degradation were downregulated and upregulated simultaneously.However, the majority of genes related to the CK signal transduction pathway were more strongly repressed in TNG67 than in TCN1; for example, LOC_Os04g36070 (RR1) and LOC_Os07g26720 (RR7) in shoots and LOC_Os03g50860 (HK4), LOC_Os02g58350 (RR3), LOC_Os01g72330 (RR4) and LOC_Os07g26720 (RR7) in roots.AHK3 is a cold-inducible A-type Arabidopsis response regulator (ARR) that is involved in cold signaling.The phenotype of arr7 Arabidopsis mutants includes tolerance to freezing [31].AHK2, AHK3, and AHK4 function as negative regulators of ABA signaling and osmotic stress signaling in both ABA-dependent and-independent pathways [50].Cytokine oxidase (CKX) transgenic Arabidopsis plants with CK deficiency show increased sensitivity to ABA and enhanced tolerances to drought and salt stresses [30].It will be interesting to further assess whether TNG67 achieves cold stress tolerance by repressing CK signaling and increasing sensitivity to ABA.
The expression of the JA biosynthesis genes LOC_Os08g04800 (OsDAD1; 2), LOC_Os12g37260 (OsLOX2; 2), and LOC_Os08g35740 (OsOPR7) was upregulated in the shoots, and that of LOC_Os03g12500 (OsAOS1) was upregulated in the roots of TNG67.The JA level may be increased to a greater extent in TNG67 than in TCN1 in response to cold stress.A previous study has shown that the addition of MeJa to TCN1 alleviates injury due to cold [7].A study of a cold-tolerant introgressed line, K354, has shown that enhanced cold resistance is associated with OsFAD7 upregulation [13].In Arabidopsis, AtFAD7 expression is induced locally by wounding, leading to the rapid and sustainable accumulation of JA [51].Our preliminary study also revealed OsFAD7 induction in TNG67 in response to cold stress.Additionally, the ratio of dienoic and trienoic fatty acids to total fatty acids was higher in TNG67.TNG67 is potentially able to synthesize JA quickly under cold stress.Jasmonates may also induce intracellular alkalinization and guard cell closure [52,53].In Arabidopsis, ABA may participate in cross-talk with MeJA signal transduction to close stomata.Furthermore, MeJA may stimulate the expression of AtNCED3, a crucial enzyme for ABA biosynthesis [54].Theoretically, if TNG67 accumulates JA faster than TCN, it should be able to synthesize more ABA, which may explain why TNG67 was able to assimilate more ABA than TCN1 at low temperatures and was resistant to cold stress.Recently, jasmonate has been shown to play a crucial role as an upstream signal in the ICE-CBF/DREB1 pathway to increase Arabidopsis freezing tolerance [34].Indeed, our microarray data showed that OsCBF1 (LOC_Os08g31580) was more strongly induced in TNG67 than in TCN1 in both shoots and roots.Furthermore, downstream RAP2.6-like (LOC_Os09g28440) was induced in the roots of TNG67 but was repressed in those of TCN1 after 3 hr of cold treatment.In addition, ZAT10-like (LOC_Os12g39400) was repressed in TCN1 roots.Thus, TNG67, but not TCN1, may transcriptionally activate JA biosynthesis-related genes and enhance gene expressions in the ICE-CBF/DREB1 regulon in response to cold stress.
Auxin has been established as a plant hormone that regulates growth and development.The role of auxin in abiotic stress has been gradually elucidated.In Arabidopsis, inflorescence gravitropism is inhibited by cold stress via auxin regulation [55].Cold stress can block the asymmetric redistribution and intracellular cycling of PIN3 and inhibit the intracellular cycling of PIN2.In Arabidopsis roots, cold stress may reduce shootward auxin transport and alter the intracellular auxin gradient [56].Therefore, this gradient may exert effects on plant growth and development to determine cold tolerance [32].According to our microarray data, the expression levels of auxin biosynthesis and signaling-related genes were preferentially induced in TNG67.This finding may indicate the involvement of auxin in the cold tolerance capacity of TNG67.
Low temperatures can alter the gene expression profiles of many hormones, and the complicated interaction or cross-talk among different hormones can act coordinately or alone to regulate the mechanisms underlying cold stress tolerance in rice.Among these phytohormones, ABA, polyamines and JA are the major determinants mediating the capacity for cold tolerance of TNG67.
Potential roles of genotype-dependent and tissue-specific TFs, especially NACs and WRKYs, in rice cold stress tolerance Recently, NAC proteins have been found to be involved in plant responses to pathogens and environmental stimuli.In the present study, NAC TFs were the predominant cold stressinduced genes detected during the early stages of cold stress in rice.For example, OsNAC6/ SNAC2 (LOC_Os01g66120) was induced after 3 hr of cold treatment in the shoots and roots of TNG67 compared with those of TCN1.OsNAC6/SNAC2 expression is induced by cold, drought, salinity stress, JA, wounding and blast disease [11,57].In addition, we identified another OsNAC3 (LOC_Os07g12340) as a drought-induced gene, and transgenic rice overexpressing this gene have been reported with improved dehydration tolerance [58].Thus, OsNAC3 may participate in cross-talk during cold and drought stresses.
The C2H2-type Zn-finger TF OsZOS11-10 (LOC_Os11g47630) was highly induced in TNG67 shoots.This result is consistent with previous reports showing that OsZOS11-10 is highly induced by cold stress similar to OsDREB1 and is negatively regulated by OsDREB1B [59].The lack of an immediate negative effect of OsZOS11-10 on OsDREB1 expression may be explained by competition with other TFs or regulation at the post-transcriptional or translational level.
Three CO-like TFs [LOC_Os06g16370 (the same as Hd1), LOC_Os08g15050, and LOC_Os10g41100) were induced in TNG67 shoots in response to cold stress.CO-like TFs are central regulators in the photoperiod pathway [60].In addition, the expression of MaCOL1, a CONSTANS-like gene in banana fruit, is significantly enhanced by abiotic and biotic stresses, such as chilling and pathogen infection [61].In addition, the peak transcript levels of DREB1/ CBF genes are high and low after the exposure of plants to low temperatures following subjective dawn and evening, respectively.These results indicate that the cold induction of DREB1/ CBF genes is circadian-gated [62].Further, LOC_Os06g16370 and LOC_Os10g41100 were highly expressed under cold stress compared with other TFs and may be important for rice cold stress tolerance.
An R2R3-type MYB gene, OsMYB2 (LOC_Os03g20090), is a transactivator that is involved in salt, cold, and dehydration tolerance in rice.OsMYB2 overexpression in rice enhances the expression of genes encoding proline synthase and transporters [63].The relative gene expression levels of proline synthase and transporters did not differ significantly between TNG67 and TCN1.However, OsMYB2 was uniquely upregulated in TNG67 roots during the early stages of cold stress in our study.Further experiments are needed to confirm the expression of proline-related genes and accumulation of proline in TNG67 and TCN1 roots.
Six genes in TNG67 and TCN1 were expressed in both shoots and roots, including LOC_Os01g53220, a heat shock factor that has been reported to be upregulated in rice cultivars that are resistant to the brown planthopper [64].This heat shock factor may play a role in the interaction between abiotic and biotic stresses in rice because cross-talk between the abiotic and biotic stress responses has been increasingly reported in recent years.
To clarify the relationship between hormones and cold-induced TFs, we further analyzed the expression of these genes using RiceXPro database (http://ricexpro.dna.affrc.go.jp/) (S9 and S10 Tables).After 3 hr of cold treatment, approximately 85% to 60% of TF expression in roots was mainly affected by JA, ABA and auxin (IAA), whereas only 40% to 25% of TF expression in shoots was affected by these hormones (S8 Fig) .Gibberellin and BR had almost no effect.Auxin may partially contribute to these effects by coordinating with ABA and JA in roots to regulate these candidate TFs (S10 Table ).
Genes involved in cytokine signal transduction may be important for reestablishing homeostasis during recovery from cold stress We found that the GO Slim category GO:0000160 [two-component signal transduction system (phosphorelay)] was specifically enriched in the shoots of TNG67 but not in those of TCN1 after recovery from cold exposure.Genes in this category are all OsRR type-A genes, including OsRR1 (LOC_Os04g36070), OsRR2 (LOC_Os02g35180), OsRR4 (LOC_Os01g72330), OsRR6 (LOC_Os04g57720) and OsRR10 (LOC_Os11g04720) (Fig 6).The expression of type-A ARR genes is known to be induced by cytokines [65].Indeed, in our study, the LOC_Os05g47840 (IPT7) gene that encodes the rate-limiting enzyme in cytokine biosynthesis was up-regulated in the roots of TNG67 but not in those of TCN1 during the recovery stage.Moreover, in Arabidopsis, the arr5, arr6, and arr7 mutants displayed elevated cold tolerance, suggesting that type-A ARRs function as negative regulators in response to cold stress [31].In response to gradual relief from cold stress, the cold response of TNG67 may be suppressed in shoots through an increase in IPT7 gene expression during the recovery stage.Moreover, CK can increase the activities of the ET biosynthesis genes ACS4 [66] and ACS5 [67] via a post-transcriptional mechanism.These results further indicate that the increased ET level observed in TNG67 during the recovery stage may have been regulated at the post-transcriptional level.We also found that GO:0015977 (carbon fixation), GO:0006457 (protein folding), GO:0008643 (carbohydrate transport) and GO:0006952 (defense response) were specifically enriched after 24 hr of recovery in the roots of TNG67 but not in those of TCN1.These results indicate that roots may play an important role during the recovery stage in response to cold stress.

Conclusions
Our transcriptome study revealed a preliminary model for the resistance of TNG67 to cold stress ( Fig 10).Following exposure to cold, the rapid induction of ABA and JA in TNG67 accelerated stomatal closure to rapidly prevent water loss.These 2 hormones may interact with one another to alter the expression of a specific set of TFs in shoots (mainly NACs).In roots, auxin may participate in cross-talk with these 2 hormones to adjust the expression of other sets of TFs, such as WRKYs, to promote cold stress tolerance in rice.Polyamine is another key component of cold tolerance in TNG67.In response to cold stress, the level of ET was reduced and then elevated after recovery in these seedlings.Ethylene may facilitate the detection of stress relief because it has an antagonistic function to ABA.CK may also be involved in this process through interacting with ET.Additionally, the activation of apoptosis-related genes may be an important signal to trigger cold stress tolerance.This process can eliminate redundant or harmful cells and recycle the constituents of cells, and it is probably a more energy-efficient means of responding to environmental stimuli.In TCN1 under cold stress, the expression of photosynthesis-related genes increased as a consequence of energy starvation.However, this process causes ROS accumulation and leads to oxidative stress with cold-induced water loss.In contrast, the expression levels of TCA cycle-and glycolysis-related genes were greater in TNG67 than in TCN1.TNG67 may compensate for decreased energy production without excessively increasing photosynthetic reactions to avoid oxidative stress.Our study expands upon the current understanding of cold stress tolerance in rice, and our findings may be used to facilitate the breeding of cold-tolerant rice.

Fig 2 .Fig 3 .
Fig 2. Venn diagram of DEGs in the shoots and roots of TNG67 and TCN1 rice seedlings exposed to cold (3 and 24 hr) and allowed to recover for 1 day (24 hr).The upregulated genes are shown above the blue line with the expression levels written in red, while the downregulated genes are shown below the blue line with the expression levels written in dark grey.DEGs with significant expression changes in the shoots (A) and roots (B) of TNG67 and TCN1 rice seedlings.Genes with a greater than two-fold change in expression compared with the control and a significant t test (p < 0.05) result were considered DEGs.doi:10.1371/journal.pone.0131391.g002

Fig 4 .Fig 5 .
Fig 4. MapMan overview of changes in the expression of genes involved in metabolism in the shoots of rice exposed to cold for (A) 3 hr and (B) 24 hr and (C) following recovery for 24 hr after cold treatment.Genes involved in different metabolic processes are shown in the main panel in dark grey, while putatively associated genes are shown in light grey.Blue indicates that gene expression was induced and red indicates that gene expression was repressed compared with the control.doi:10.1371/journal.pone.0131391.g004

Fig 6 .
Fig 6.Representative gene expression profiles of hormone biosynthesis and deactivation signaling pathways in TNG67 and TCN1 shoots and roots.Only DEGs with significant changes in gene expression are shown in the heat map.Each color spot reflects the expression of a corresponding gene.Red indicates high levels of gene expression, green represents low levels, and gray indicates that no signal was detected on the microarray chip for the assessed gene probe.doi:10.1371/journal.pone.0131391.g006

Fig 7 .Fig 8 .
Fig 7. Venn diagram of upregulated (upper area above the blue line, expression levels written in red) and downregulated transcription factors (TFs) (area below the blue line, expression levels written in dark grey) in (A) the shoots and (B) roots of TNG67 and TCN1 at different time points of cold stress.Genes with a greater than two-fold change in expression compared with the control and a significant t test (p < 0.05) result were considered DEGs.doi:10.1371/journal.pone.0131391.g007 LOC_Os02g03030) TF belonged to a module containing HMGB (LOC_Os08g01100) (S5 Fig).When we excluded these 3 genes in the 2 main modules and then re-analyzed the data, the 2 largest modules contained AUX/IAA (OsIAA23; LOC_Os06g39590), bHLH (LOC_Os01g39330), and 2 NAC TFs (SNAC2; LOC_Os01g66120 and LOC_Os05g34830) in the same module (S6 Fig).In addition, 4 modules were detected after inputting 30 root TFs as query genes.One module contained bHLH (LOC_Os02g47660; orthologous to AtbHLH063), and the other had NAC (SNAC2; LOC_Os01g66120).MYB (LOC_Os07g48870) was present in the smallest module.Interestingly, we found 4 WRKY TFs [OsWRKY 1v2 (LOC_Os01g14440); OsWRKY53 (LOC_Os05g27730); OsWRKY71 (LOC_Os02g08440) and OsWRKY24 (LOC_Os01g61080)] in the same module (S7 Fig).To confirm the specific species-, tissue-and cold-induced gene expression of these TFs, gene-specific primers were designed for the above-mentioned 13 genes, and qRT-PCR was conducted (Fig 9).The results obtained following exposure to low temperature stress for 3 hr were in accordance with those of module analysis, as shown in Fig 9.

Table 1 .
(Continued) Table.The proteasome map (S2 Fig) revealed the significantly increased expression of genes related to E2, E3 ligase, and the 20S proteasome in TNG67