Global Expression Profiling of Low Temperature Induced Genes in the Chilling Tolerant Japonica Rice Jumli Marshi

Low temperature is a key factor that limits growth and productivity of many important agronomical crops worldwide. Rice (Oryza sativa L.) is negatively affected already at temperatures below +10°C and is therefore denoted as chilling sensitive. However, chilling tolerant rice cultivars exist and can be commercially cultivated at altitudes up to 3,050 meters with temperatures reaching as low as +4°C. In this work, the global transcriptional response to cold stress (+4°C) was studied in the Nepalese highland variety Jumli Marshi (spp. japonica) and 4,636 genes were identified as significantly differentially expressed within 24 hours of cold stress. Comparison with previously published microarray data from one chilling tolerant and two sensitive rice cultivars identified 182 genes differentially expressed (DE) upon cold stress in all four rice cultivars and 511 genes DE only in the chilling tolerant rice. Promoter analysis of the 182 genes suggests a complex cross-talk between ABRE and CBF regulons. Promoter analysis of the 511 genes identified over-represented ABRE motifs but not DRE motifs, suggesting a role for ABA signaling in cold tolerance. Moreover, 2,101 genes were DE in Jumli Marshi alone. By chromosomal localization analysis, 473 of these cold responsive genes were located within 13 different QTLs previously identified as cold associated.


Introduction
The inherent ability of plants to endure low temperatures affects both their geographical distribution and overall productivity. Many cereal crops from temperate regions, e.g., wheat, barley, rye and oat, have evolved efficient protection mechanisms to tolerate freezing [1]. This tolerance is acquired through a process known as cold acclimation, which occurs at low, but above-zero temperatures [2]. On the other hand, the highly important staple food, rice, like many other tropical plants, is chilling sensitive and does not survive freezing temperatures even after cold acclimation. However, chilling susceptibility varies substantially among cultivars of the same species [3], indicating that even chilling sensitive crops, at least to a certain degree, can acclimate to cooler temperatures. For example, japonica rice cultivars show a greater tolerance to chilling than indica cultivars, although variations within japonica exist [4].
Rice is mainly grown in warm climates (.25uC), but exposure to low temperature (LT) is common for rice cultivated in temperate zones or at high elevations in several regions of Europe, South Asia and Southeast Asia. It is, therefore, of outstanding scientific and economic interest to understand chilling tolerance, identify key regulatory components of chilling acclimation, resolve mechanistic differences between tolerant and sensitive rice and apply this knowledge in the development of new, more tolerant rice cultivars.
It is known that cold acclimation leads to physiological and metabolic changes in cell and tissue structures as a result of an extensive reprogramming in gene expression [2,5,6,7]. A large number of genes that are differentially expressed during cold acclimation have been identified and characterized in important cold hardy cash crops like wheat (Triticum aestivum) [8], barley (Hordeum vulgare) [9,10] and oat (Avena sativa) [7,11]. However, the plant species most extensively studied during cold stress and acclimation is a non-crop plant, Arabidopsis thaliana, where global transcriptional profiling experiments have identified several cold responsive genes [12,13,14,15,16,17]. The complexity of genetic re-programming upon cold stress has also been demonstrated by various bioinformatics approaches [15,18,19,20,21].
Previously, Rabbani et al. [22] used a rice cDNA microarray of 1,718 ESTs and identified 36 cold responsive genes in two-week old seedlings of chilling tolerant rice Nipponbare (spp. japonica) exposed to +4uC for 24 hours. Cheng et al. [23] used a rice cDNA microarray of 5,855 unique ESTs and identified 121 cold responsive genes in 10 days old seedlings of the chilling tolerant rice CT6748-8-CA-17 (spp. japonica) treated at +10uC for up to 24 hours. In another study, Oda et al. [24] used the 44K Agilent oligonucleotide microarray to compare two japonica cultivars Sasanishiki and Hitomebore exposed to low temperature stress (19uC) at the reproductive stage. Microarray analysis of anthers from the two cultivars led to the identification of 356 differentially expressed genes in either or both cultivars. Yun et al. [25] analyzed genes induced by chilling stress (+10uC) in Nipponbare using microarrays representing 40,000 genes and identified 8,668 differentially expressed genes. Mittal et al. [26] performed microarray analysis of the indica rice Pusa Basmati that was cold stressed at +5uC and identified 924 differentially expressed genes. Zhang et al. [27] performed comparative microarray analysis of a chilling tolerant rice cultivar (LTH; japonica) and a chilling sensitive rice cultivar (IR29; indica) and showed that although the early response to low temperatures was similar in the two cultivars, the genes that were expressed at the later time points belonged to substantially different functional categories. These studies clearly show that chilling tolerance varies among rice cultivars and that many rice genes respond to low temperature stress.
In order to gain new insights into cold stress response in rice, the main aim of this study was to conduct a global cold (+4uC) responsive gene expression profiling of the Nepalese highland rice cultivar Jumli Marshi (JM) (spp. japonica). This rice is grown in the district Jumla located at an altitude of up to 3,050 m in Nepal. The average maximum and minimum temperature in the region is +21uC and +4uC, respectively [28]. JM has been the preferred variety grown in this district for several decades, thus making Jumla the highest and coldest region in the world for commercial cultivation of rice. The fact that JM can be grown at this altitude while maintaining productivity corresponding to 60% of the average Nepal rice productivity per hectare [28] implies that JM is indeed chilling tolerant.
Obviously, JM has developed ways to protect itself from cold stress. Identifying genes involved in the underlying molecular mechanisms could reveal new insights into how cold tolerance is achieved. This will ultimately enable development of new cultivars with improved cold tolerance. In this work, we performed transcriptome analysis of JM under cold stress and identified 4,636 differentially expressed genes. To our understanding, this is the first report on the global analysis of cold responsive genes in Nepalese highland rice during early cold stress (,24 h, +46C).

Plant Material and Growth Conditions
Regular growth condition: Seeds of two rice cultivars, Oryza sativa, ssp. japonica, cv. Jumli Marshi (JM) and ssp. indica, cv. IR64 (IR64), were first soaked in water for 16 hours at room temperature and thereafter grown on standard soil in 14 hours photoperiod, with a day/night air temperature of 25uC/20uC and 250 mmol m 2 s 21 light.
Cold condition. At mid-day, three weeks old plants were transferred to growth chambers (Percival) in the same photoperiod, but with an air temperature of 4uC and a light intensity of 100 mmol m 2 s 21 . Pooled leaf tissue from five individual Jumli Marshi plants were harvested at 0, 0.5, 2, 4, 8, and 24 hours, frozen in liquid nitrogen and stored at 280uC until further analysis.

Chlorophyll Fluorescence Measurements
Chlorophyll fluorescence was measured with the portable chlorophyll fluorometer PAM 2000 (Heinz Walz GmbH, Germany) and the photosystem II efficiencies F v /F m = (F m -F 0 )/ F m were estimated as per the manufacturer's instructions. Plants were dark acclimated for one hour before taking measurements from up to 20 individual seedlings.

Data Preparation and Analysis
Total RNA was extracted from JM leaf tissue with TRIZOL reagents (Invitrogen) according to the manufacturer's protocol and purified by RNeasy MinElute Cleanup Kit (Qiagen). The RNA quality and concentration was measured using Agilent 2100 BioAnalyzer and Nanodrop ND-1000. Biotinylated target cRNA was prepared from 4 mg of total RNA following the manufacturer's specifications (Affymetrix). The samples were then hybridized to Affymetrix GeneChipH Rice Genome Arrays, which contain probes to query ,51K transcripts representing both japonica and indica cultivars. The chips were thereafter washed and stained in a GeneChipH Fluidics Station 450. Scanning was carried out with GeneChipH Scanner 3000 and image analysis was performed using GeneChipH Operating Software. Two biological replicates were analysed per time point. The CEL files were submitted to ArrayExpress with the accession number E-MEXP-3718.
Data was processed using Bioconductor [32] in R v2.14. Raw CEL files were background corrected with the GCRMA method and quantile normalized using the Bioconductor package affyPLM v1.30 [33]. Probe-set present/absent calls were calculated with the mas5calls method in the AffyBatch package [34]. Probe sets that were marked present in at least one of the samples were considered for further analysis. Probe sets with IQR greater than the median IQR of all probe sets were selected with the genefilter v1.36 package. To identify differentially expressed probe sets, a linear contrast matrix was built between the control and the cold-treated samples using the Limma v3.10 package [35]. Differentially expressed genes (DEGs) were identified with the empirical Bayes method in the Limma package using default parameter settings. A Benjamini-Hochberg corrected p-value ,0.05 was set as the significance threshold. Thereafter, only probe sets with a unique RAP OS ID were retained for further analysis.
MapMan annotations for Oryza sativa (spp. japonica) genes (v1.0) were downloaded (www.mapman.gabipd.org), and annotations for the DEGs were extracted using a custom Perl script. Hierarchical clustering was performed with Cluster 3.0, using the average linkage method and Pearson correlation as similarity measure. Dendrograms were generated with Java TreeView 1.1.6. Clustering was also performed with Short Time-series Expression Miner (STEM, version 1.3.8) [36] using default parameters.

JM has Higher Cold Tolerance than IR64
To confirm the higher cold tolerance in JM compared to lowland rice, three weeks old seedlings grown at regular conditions were transferred to cold conditions and exposed to +4uC for three days (see Materials and Methods). Plants were then allowed to recover for two weeks in regular growth conditions, and their viability was visually estimated. The results confirmed that JM is chilling tolerant while IR64 is not since all IR64 plants wilted and died while the JM plants recovered from the cold stress (Figure 1a, 1b).
To quantify the stress levels in the two cultivars, chlorophyll fluorescence signals were measured. The ratio between the variable (F v ) and maximum (F m ) fluorescence signals reflects the efficiency of photosystem II. Plants were first dark-acclimated for 1 hour, and the fluorescence was measured (PAM 2000 fluorometer; Heinz Walz GmbH, Germany) just before transferring the plants to +4uC (0 days), and then after 1, 2 and 3 days in the cold. Readings from 50 individual seedlings were collected, and the experiment was repeated twice. Prior to cold exposure, F v /F m ratios were similar (,0.8) in both JM and IR64 plants, indicating that all plants were physiologically healthy. After one day in cold, the F v /F m ratio in IR64 fell to 0.14, while, in JM, the decrease was relatively small (0.64). The differences between the averages were significant at p,0.0001 (Student's t-test). After two days in cold, the F v /F m ratio reached near zero in IR64, while, in JM, it reached 0.55 ( Figure 1c). These results show that JM tolerates cold stress significantly better than IR64.

Identification of Differentially Expressed Genes under Cold Stress in JM
Since JM is chilling tolerant, it was chosen for global gene expression profiling under cold stress. Three weeks old plants grown under regular conditions were moved to cold conditions (+4uC) at mid-day and pooled leaf tissues from five individual plants were harvested at six different time points: 0, 0.5, 2, 4, 8 and 24 hours.
Microarray data analysis revealed that 4,636 genes were significantly differentially expressed (Table S1). The complex transcriptional response pattern to early cold stress in JM can be visualized by hierarchical clustering (Figure 2). Two major clusters were formed, one with 1,490 up-regulated genes, and the other with 3,146 down-regulated genes. There were 183 genes (148 up, 35 down) DE by at least two folds (log 2 ) within two hours of cold exposure. Most common annotations in this set were genes with unknown function (83 genes), regulation of transcription (36 genes), protein degradation or post-translational modification (16 genes) and signaling (8 genes) ( Figure S1). Out of the 36 genes with transcription factor (TF) activity, the most common TF families were AP2/EREBP (12 genes), WRKY (5 genes) and bHLH (4 genes). The most highly induced genes within 2 hours were OsDREB1B (10 log 2 folds) and OsDREB1A (8.9 log 2 folds). Quantitative real-time RT-PCR measurements of six genes ( Figure 3 and Figure S2) showed good correlation with the microarray data.

DREB1 Genes are Similarly Expressed in JM and IR64
DREB1s encode transcription factors that belong to the AP2 family of DNA binding proteins. Several experiments in many different plant species have shown that DREB1A, DREB1B and DREB1C are induced within 30 minutes of exposure of plants to  cold. For example, over-expression of DREB1s in Arabidopsis resulted in increased freezing tolerance of non-acclimated transgenic plants compared to the non-transformed control plants [37]. Thus, the mRNA expression levels of the three DREB1 genes were measured by quantitative RT-PCR in JM and IR64. The results revealed that all three genes were highly induced (.5 folds) within 2 hours in both JM and IR64, although, the expression levels were higher in IR64 than in JM ( Figure 3). Previously Zhang et al [38] identified 22 genes as part of the OsDREB1c regulon. In JM, 18 of these 22 genes were DE upon cold stress indicating the critical role of OsDREB1c regulon under cold tolerance in JM ( Figure S3).

Comparative Analysis of Microarray Data from Four Rice Cultivars Reveals Genes that are Differentially Expressed only in Chilling Tolerant Rice
Comparative analysis with previously published rice microarray data was done to identify genes that are differentially expressed (DE) only in chilling tolerant rice. Zhang et al. [27] performed microarray data analyses of the chilling tolerant rice cultivar LTH (spp. japonica) and the chilling sensitive rice cultivar IR29 (spp. indica). Mittal et al. [26] performed microarray data analysis of the rice cultivar Pusa Basmati (PB1, ssp. indica). PB1 is mostly cultivated in the North-Western part of India [39] in the warm season May-November and thus is considered chilling sensitive.
In this work, comparative analysis was done between JM and the three previously published datasets and the genes that were DE by at least 2 folds in at least one of the four rice cultivars were studied further ( Figure 4). The analysis identified 182 genes that were DE by at least two folds in all four rice cultivars. These genes were termed Common Cold Induced (CCI) ( Figure 4, Table 1). The CCI gene list includes several well-known cold induced genes induced by 2 folds within 2 hours after induction including OsDREB1C, OsDREB1G, OsWRKY71, OsNAC3 and also genes induced later than 2 hours including OsWRKY1 and OsNAC4. The comparative analysis also identified 511 genes that were only expressed in the two chilling tolerant cultivars JM and LTH. These 511 genes were termed as Cold Induced in Tolerant cultivars (CIT). Expression levels of three genes from the CIT group were validated by quantitative RT-PCR ( Figure S2). As the CIT genes are DE only in the chilling tolerant cultivars, they may play a significant role in chilling tolerance in rice. A third set with 2,102 genes was only DE in JM. These were denoted as JM Only (JMO). As JM is cultivated at altitudes of up to 3,050 m, further understanding of the JMO genes may help to identify novel mechanisms for survival of plants in these environmental conditions. Therefore, the JMO genes were further analyzed in this study. A complete list of genes in different groups is in Table  S2.

Identification of over-represented cis-elements in CCI
and CIT gene sets. Clustering was performed with STEM [36] to identify the CCI and CIT genes that were significantly coexpressed. The clustering of CCI genes resulted in two statistically significant clusters with 127 and 14 genes each, whereas the CIT genes were grouped into six clusters ( Figure 5). Promoter analyses of genes in the eight clusters by the Osiris tool [40] lead to the identification of over-represented cis-elements/motifs previously known to be involved in cold regulation ( Table 2). For example, CCI cluster 1 contains the well-known cold regulated motifs ABRE and DRE/CRT, and eleven other previously known motifs. This suggests a cross-talk between various transcription factor families during cold stress. Among the six CIT clusters, known over-represented motifs were only identified in cluster 3. In the remaining five clusters, no previously known motifs could be identified, suggesting that yet unknown regulons are involved in regulating many of the CIT genes.
To identify new putative motifs, gene promoters from CIT clusters 1, 2, 4, 5 and 6 were analyzed using the Element algorithm [41]. The 1 kb upstream region was analyzed for each promoter, and the top five hits were retained (Table 3). This analysis showed that while a sub-set of the over-represented words were entirely or partially similar to previously known cis-elements, the remaining words did not belong to any previously known cis-elements and thus are newly identified putative cis-elements (Table S3).
Comparative analysis of functional annotations. MapMan annotations common in the CCI and CIT gene sets were identified. Examples of such annotations are regulation of transcription, posttranslational modification, protein degradation, development and calcium signaling. Interestingly, the CIT gene set also included additional annotations such as signaling in receptor kinases, protein synthesis and RNA binding ( Figure 6). MapMan annotations were also obtained for all 4,636 DEGs (Table S4). The top three annotations were protein (645), RNA (574) and miscellaneous (227) (Figure 7). Analysis of annotations that were highly abundant in JM and LTH compared to chilling sensitive PB1 and IR29 showed that amino acid metabolism, RNA transcription, transport of sugars,   photosynthesis light reaction and secondary metabolisms were the top five over-represented annotations in the cold tolerant cultivars (Figure 8).

Several TF Families are DE in JM
Annotations were also downloaded from the Plant Transcription Factor database (PlantTFDB) [42] to assess the most prominent cold responsive transcription factor (TF) families. In JM, 277 (5.9%) of the DEGs belonged to 37 represented TF families. TFs from well-known cold stress associated families like AP2/ERF (34 genes), NAC (24 genes), MYB-like (21 genes), bHLH (19 genes) and WRKY (14 genes) were among the overrepresented families (Figure 9, Table S1). Hierarchical clustering of the AP2/ERF TFs shows that most of these genes are either highly up-or down-regulated within 24 hours of cold stress and thus form two major clusters based on their expression ( Figure 10). OsDREB1s are among the most highly induced genes within 30 minutes of exposure to cold.

Significant Changes Occur in Hormone, Signaling and Carbohydrate Metabolism Associated Genes
A number of signaling, as well as hormone regulation pathways, have previously been shown to be prominent in cold stress responses [43,44]. Genes involved in different hormone pathways were identified from the MapMan annotations. Among pathways coupled to cold stress were auxin (indole-3-acetic acid) (36 genes), abscisic acid (ABA) (15 genes), cytokinin (13 genes) and ethylene (13 genes) (Figure 11a). Genes related to calcium signaling (40 genes) and the MAP kinase cascade (8 genes) were also differentially expressed (Figure 11b), as well as genes involved in sugar synthesis and degradation pathways, including sucrose degradation (11 genes), trehalose synthesis (8 genes) and starch degradation (8 genes) (Figure 11c).   HO 2 ? _ , OH . , ROOH, ROO . , and RO . [45]. ROS are highly reactive and induce deleterious effects on membranes that lead to ion leakage and also cause damage to proteins, DNA and lipids. Ultimately these effects result in cell death [45]. It has been shown previously that in the rice cv. Nipponbare, H 2 O 2 levels increase within 1.5 hours of chilling stress (+10uC) [25]. To counteract the increased ROS levels, the cell produces scavengers such as peroxidases, catalases, ascorbate, glutathione, superoxide dismutase, glutaredoxins and thioredoxins that detoxify ROS [45]. In JM, 36 genes associated with ROS scavenging were cold induced. Among peroxidases, nine genes were differentially expressed, two up-regulated after eight hours and the remaining seven genes down-regulated after two to eight hours. Two glutathione peroxidases were also differentially expressed, one gene up-and one down-regulated after 24 hours, while glutathione S transferase was induced with two genes up-regulated at 8 hours and three down-regulated after four hours of cold stress. In addition, twenty genes belonging to glutaredoxins and thioredoxins were differentially expressed, with four genes up-regulated within the first two hours of cold stress and the remaining up-regulated during the 24 hour assay period. Out of the remaining 16 genes, one was upregulated after 24 hours and the others were down-regulated ( Figure 12, Table S5). Of the 36 genes, eleven genes were in the JMO set (underlined in Figure 12), one (Os01g0667900) was in the CCI set and four (Os03g0762300, Os11g0284900, Os12g0188700, Os08g0378900) were in the CIT set. Overall, this suggests differences in the ROS scavenging mechanisms in different rice cultivars.

Cellular Components Play Distinct Roles under Cold Stress
Organelles and sub-cellular components are highly affected by cold stress [46]. For example, in Arabidopsis, 184 nuclear proteins [47], 43 chloroplast proteins [48] and 38 plasma membrane proteins were identified as differentially expressed under cold stress [49]. Mitochondria, which are the prime sites for ROS production in abiotic stress, regulate ROS levels through their energy dissipating systems [50]. To evaluate if sub-cellular components play an active role during cold stress in JM, a functional annotation analysis was performed for DEGs encoding proteins predicted to be either organelle localized or associated to a cellular structure (hereon referred to as 'orgLoc genes') ( Figure 13, Table  S1).
The analysis showed that orgLoc genes from various cellular components were differentially expressed at a significant level during cold stress. The percentage of induced genes for each subcellular component or organelle was: cell wall (18%), plasma membrane (17%), nucleus (16%), cytoplasm (13%), endoplasmic reticulum (13%), mitochondria (12%), chloroplast (11%), vacuole (10%), Golgi apparatus (9%) and cytoskeleton (9%) (Figure 13). Moreover, several genes encoding plasma membrane receptors were also cold induced. A functional analysis of the orgLoc genes suggested distinct roles for different sub-cellular components and organelles during cold stress. For example, cold induced alterations in transport of metabolite process were found in the plasma membrane, mitochondria, chloroplast and vacuole, and signaling process in the plasma membrane, mitochondria and cytoplasm. Differential regulation in auxin signal transduction process could be seen in the plasma membrane, while, in the nucleus, auxin, ABA, cytokinin and ethylene signal transduction process were observed. Metabolism process of several compounds was also seen in the mitochondria, endoplasmic reticulum, cytoplasm, chloroplast and cell wall while differences in redox activity process were seen mainly in mitochondria and chloroplasts. Together this strongly suggests that the buildup of tolerance to low temperature stress involves many different components and that these components act synergistically.

Identification of Cold Responsive Genes in QTLs Associated with Cold Stress
Previous studies have identified 22 different QTLs associated with cold stress in rice (hereon referred to as 'cold QTLs'). To assess whether any of these QTLs overlap with stress response genes in JM, all 22 QTLs were downloaded from the Gramene database [51] and putative genes within the QTLs were identified and checked for homology to the JM DEGs.
This showed that 473 of the DEGs were present in 13 of the cold QTLs ( Figure 14, Table S6) while no DEGs were found in the remaining 9 QTLs. MapMan annotations of the genes present in the 13 cold QTLs suggest distinct functional roles of these QTLs. The most common annotations were protein (62 genes), RNA (54 genes), signaling (27 genes), miscellaneous (22 genes) and transport (20 genes). All 13 QTLs also contained genes with yet unknown molecular functions. Of the 473 genes, 232 genes were found in the JMO set, while 40 and 12 genes were found in the CIT and CCI sets respectively (Table S7). The fact that several different QTLs contained DEGs with different functions suggests that cold tolerance is a multigenic trait. To increase cold tolerance by breeding, several different processes would, therefore, need to be successfully transferred to the same genome. As the cold QTLs were identified through genetic studies, it is highly likely that several of the cold induced genes located within the identified QTLs are critical for cold stress tolerance also in JM.

Discussion
Low temperatures negatively affect both productivity and yield of rice grown in colder regions of the tropics. However, over the years, farmers and breeders have selected cultivars that are more Table 3. Over-represented words in the 1 kb upstream promoter regions were identified using the Element algorithm [41]. chilling tolerant, resulting in rice cultivars that can even cope with cold stress to some extent. Since we were interested in the molecular mechanisms underlying this adaptation, we selected a chilling tolerant rice variety Jumli Marshi (JM) as a model and studied its transcriptional responses to cold stress. JM is a locally adapted endogenous variety that is routinely grown at up to 3,050 m in the Jumla district of Nepal. As a consequence of the high altitude, this variety is commonly exposed to temperatures down to or below 4uC during the early developmental stages in the growing season. JM is commonly used as a parental donor in various rice breeding programs worldwide for the development of new chilling tolerant rice cultivars. Chilling tolerant Nepalese rice cultivars Manjushree-2 and Khumal-8 were developed with JM as one of the parents [52]. JM was also successfully used as a parental donor for the development of upland chilling tolerant rice FOFIFA-172 in Madagascar [53]. JM has mostly been crossed with indica cultivars for the development of new chilling tolerant cultivars. Rice cultivar Khumal-8 is such an example where the parents are JM (spp. japonica) and IR36 (spp. indica). A study on the pedigrees of 28 Nepalese rice cultivars showed that 65.7% of the ancestors were of indica type and only 2.8% were japonica [54] again indicating that the parental types have mostly been indica Thus, in general a more detailed molecular comparison of japonica and indica subtypes is necessary for increased understanding of physiological, molecular and evolutionary aspects of chilling tolerance in rice. A few studies where chilling tolerant japonica rice cultivars have been compared to chilling sensitive indica ones exist [25,27], but more comparisons with different cultivars are needed. Due to their superior agronomical characteristics, indica cultivars are most commonly used in Nepalese rice breeding programs, while, japonica cultivars are valued mainly for their chilling tolerance characteristics. Thus, comparative analysis of chilling sensitive indica and chilling tolerant japonica rice will facilitate in understanding of the chilling tolerance mechanisms and identify markers for breeding. In this work, a comparative transcriptome analysis was done between chilling tolerant upland Nepalese rice JM with a chilling tolerant japonica (LTH), and the two chilling sensitive indica rice cultivars PB1 and IR29. Adding chilling sensitive japonica rice to this analysis would have been of interest, but unfortunately there is no publicly available microarray data generated from chilling sensitive japonica rice.
Many of the cold induced genes derived in this study have previously been associated to cold stress in other systems [14,16,17,23,26,27,55,56]. For example, DREB1s act as up-stream regulators during cold stress and are induced several folds within 30 minutes of cold treatment in several plant species [14,16,26,27,55,57,58,59,60,61,62]. Thus, many of the early cold  induced genes seen in cold hardy plants are present in rice and the overall regulation of the upstream genetic regulators is also conserved. This indicates that some of the differences between chilling and freezing tolerant crops is not in the primary activation of cold acclimation pathways but may lie in the second wave of activation of downstream protection mechanisms or the existence of CBF-independent regulons in the chilling tolerant rice cultivars. However, the OsDREB1C regulon targets [38] are cold induced in both JM and IR64 ( Figure S3b), which further strengthens the hypothesis that CBF-independent regulons play a critical role in chilling tolerance in rice.
In a previous study, Yun et al. [25] identified an early response putative regulatory network in japonica rice upon chilling stress at 10uC and suggested a strong integration of defense and growth related responses. Several well-known cold stress responsive TFs were differentially expressed including AP2, bZIP, MYB, WRKY, bHLH and NAC families. Functional annotations such as transport, signaling, defense, transcription and apoptosis were over-represented. In this work, although JM was cold stressed at 4uC, functional annotations common under chilling stress (10uC) were also found to be common under cold stress, indicating an overlap in regulatory mechanisms in chilling and cold stress. In a transcriptome profile comparison between LTH and IR29, Zhang et al. [27] identified 1,256 genes to be up-regulated only in LTH upon cold exposure. They also identified several functional annotations specific to either LTH or to IR29 or common to both. The results from this work further highlight genes, functional annotations and regulatory motifs that are specifically induced in JM.
Chlorophyll measurements suggest that IR64 plants suffer greatly within 48 hours of cold stress indicating that genes induced within the first 24 hours are critical for plant survival. This study identified 4,636 genes that were significantly differentially expressed within 24 hours upon cold stress in JM. Functional annotations of these genes showed that many of them are involved in gene regulation, post-translational modification, signaling, transport, development and metabolism. Comparative analysis with one other chilling tolerant rice and two chilling sensitive rice resulted in identification of 182 genes that were cold induced in all four rice cultivars (CCI genes). Also, 511 genes were induced in the two chilling tolerant cultivars (CIT genes). Promoter analysis of the CCI and CIT clusters gave further insight into similarities and differences among different cultivars. Motifs binding to CBF/DRE or ABRE TFs were over-represented in the promoters of CCI genes, thus suggesting that well-known cold induction pathways are conserved in rice cultivars differing in their chilling tolerance. Interestingly, ABRE elements were over-represented in the promoters of CIT genes, but CBF/DRE elements were not. This  indicates that there are differences in the gene regulatory mechanisms in the chilling tolerant and sensitive lines. These differences could either be activation of parallel pathways that are independent of CBF regulation, or activation of pathways that are downstream of the CBF regulatory network. A further difference in the response of chilling tolerant and sensitive rice cultivars are the functional annotations of the DE genes. The top three annotations that were abundant in JM compared to the chilling sensitive cultivars were amino acid metabolism, RNA transcription and sugar transport.
Gene expression profiling under cold stress in Arabidopsis showed that the AP2/ERF family of TF genes was the most represented [63]. In japonica rice, it was shown that molecular responses triggered by oxidative signals are critical for survival under low temperature stress [25]. Under chilling stress, the chloroplast electron transport chain is over-reduced, resulting in increased production of ROS [64]. Less effective mitochondrial electron transport during cold stress can also lead to ROS formation [65], which can be detoxicated by different scavenging antioxidant enzymes produced by the plants. ROS scavengers, such as catalase   [66], ascorbate peroxidase [67], glutathione reductase [68] and glutathione peroxidase [69] have been shown to play a role in eliminating free oxygen radicals during chilling stress. In rice, putative ROS mediated regulatory modules under chilling stress (+10uC) have been identified, and were suggested to be independent of DREB and ABA regulons [23]. In this work, genes involved in production of peroxidases, glutathione peroxidases, glutathione S-transferases, glutaredoxins and thioredoxins were significantly differentially expressed (Figure 12), indicating that ROS scavenging mechanisms are active in JM during cold stress (+4uC).
This study also shows that various cellular components play distinct roles during cold stress. Several orgLoc genes from nucleus, plasma membrane, endoplasmic reticulum and other components were differentially expressed in JM. The analysis reflects the physiological processes and downstream molecular mechanisms that possibly occur in the cellular components under cold stress. The plasma membrane acts as a primary barrier between the cytoplasm and the extracellular space and plays a critical role in cell integrity, signaling and transport [70]. In this work, plasma membrane orgLoc genes were involved in transport and the auxin signal transduction process. The cell wall provides cell integrity and support, and cell wall orgLoc genes were involved in the cell wall modification and organization process, reflecting changes occurring in the cell wall upon cold stress. The mitochondria and chloroplast are involved in metabolism of several compounds as well as energy generation from the electron transport chain. Again, several orgLoc genes with various functions differentially expressed in these cellular components indicate critical roles for these under cold stress. An overview of the processes taking place in each component during cold stress is shown in Figure 13.
Although efforts have been made to introduce increased chilling tolerance in rice by traditional breeding, the progress has been limited because chilling tolerance is a multigenic trait and involves several layers of protection to shield the plants from the stress. However, a number of QTLs that are statistically significantly associated with an increase in chilling tolerance have been identified through years of breeding. A list of such 'cold QTLs' from rice is summarized in Table S6. Since these QTLs are coupled to chilling tolerance, it is of interest to identify cold induced genes located within them. In a previous study, 445 cold  Table S1. doi:10.1371/journal.pone.0081729.g013 induced genes from rice were located to 21 QTLs [27]. In this work, 473 of the cold induced JM genes were found to be located within 13 of the 22 previously identified cold QTLs ( Figure 14, Table S6). Of these 473 genes, 232 genes were found in the JMO set, while 40 and 12 genes were found in the CIT and CCI sets respectively (Table S7). A more substantial increase in chilling tolerance is most likely only possible by a combined contribution from multiple QTLs. Thus, this and other similar studies on QTLs will facilitate the development of new markers and will further our understanding of chilling tolerance mechanisms in rice.

Conclusion
Jumli Marshi is a highly adapted upland chilling tolerant Nepalese rice variety and is used extensively in breeding for new chilling tolerant rice cultivars both in Nepal and other countries. In this study, a global expression profiling was done to identify cold induced genes in JM. Comparative analysis resulted in identification of 2,101 genes DE only in chilling tolerant JM rice. From the comparison of JM and LTH, it can be concluded that there are significant differences in the way the two chilling tolerant japonica rice cultivars respond to cold, which further emphasizes the importance of including locally adapted landraces in scientific and breeding programs. In addition, in JM, 473 cold induced genes localized in previously described cold associated QTLs were identified and among these, 232 genes were DE only in JM. Further validation of these genes may lead to the development of new molecular markers for breeding of rice cultivars with increased chilling tolerance.