Transcriptome Response Mediated by Cold Stress in Lotus japonicus

Members of the Lotus genus are important as agricultural forage sources under marginal environmental conditions given their high nutritional value and tolerance of various abiotic stresses. However, their dry matter production is drastically reduced in cooler seasons, while their response to such conditions is not well studied. This paper analyzes cold acclimation of the genus by studying Lotus japonicus over a stress period of 24 h. High-throughput RNA sequencing was used to identify and classify 1077 differentially expressed genes, of which 713 were up-regulated and 364 were down-regulated. Up-regulated genes were principally related to lipid, cell wall, phenylpropanoid, sugar, and proline regulation, while down-regulated genes affected the photosynthetic process and chloroplast development. Together, a total of 41 cold-inducible transcription factors were identified, including members of the AP2/ERF, NAC, MYB, and WRKY families; two of them were described as putative novel transcription factors. Finally, DREB1/CBFs were described with respect to their cold stress expression profiles. This is the first transcriptome profiling of the model legume L. japonicus under cold stress. Data obtained may be useful in identifying candidate genes for breeding modified species of forage legumes that more readily acclimate to low temperatures.


INTRODUCTION
Legumes (Fabaceae) are the primary source of plant protein for human consumption and in livestock feed, as well as key components of natural and agricultural ecosystems. Members of the Lotus genus, which includes about 120 species around the world, are particularly well known due to their elevated adaptability to marginal environmental conditions. This feature has made species of the Lotus genus popular as alternative forage in South America and Australia, as well as for dune revegetation and the reclamation of soil that has been burned or contaminated by heavy metals (Escaray et al., 2012). Within the genus, L. japonicus has become a model legume (Handberg and Stougaard, 1992), having been extensively used in abiotic stress Sainz et al., 2010;Babuin et al., 2014) and root nodulation studies (López et al., 2008;Li et al., 2014); this utility mainly results from particular genomic features that make it especially useful for recently developed functional genomics techniques (Sato et al., 2008;Fukai et al., 2012;Urbañski et al., 2012).
Even though forage Lotus species are considered cool season plants (Blumenthal and McGraw, 1999), dry matter production is in fact minimal during colder seasons at high latitudes (Bullard and Crawford, 1995;Halling et al., 2004), causing substantial economic losses (Wheeler et al., 2000;Thakur et al., 2010). It is therefore important for breeders to understand the mechanisms through which these plants tolerate cold and freezing. Previous research has focused on the cold stress response of Arabidopsis and has included both genetic and biochemical approaches (Mckown et al., 1996;Hannah et al., 2005;Kaplan et al., 2007). However, few studies have addressed legumes (Lucau-Danila et al., 2012;Dinari et al., 2013).
Plants have various ways in which they respond to and tolerate cold stress. These include changes in the composition, structure, and function of the plasma membrane, the synthesis of cryoprotectant molecules (soluble sugars and low-molecularweight nitrogenous compounds, such as proline), and an increase in the scavenging activity of reactive oxygen species (ROS; Steponkus, 1984;Apel and Hirt, 2004;Wang et al., 2013). Cold temperatures also induce distinct secondary metabolic pathways. In particular, flavonoids have been proposed as a primary target of study for understanding this phenomenon; this large family of metabolites, synthesized via the phenylpropanoid pathway, has been related to photoprotection, cold hardiness, drought resistance, and antioxidative ability (Christie et al., 1994;Chalker-Scott, 1999).
Most cold stress responses are due to changes in gene expression, for which many transcription factors (TFs) have already been identified (Chinnusamy et al., 2007). The most well documented pathways involve a class of DREB/CBF TFs that are upstream regulated by several proteins, such as ICE1, MYB15, and ZAT12 (Shinozaki et al., 2003;Maruyama et al., 2009); the genes that are regulated through these pathways are collectively identified as cor genes. A large number of these genes encode proteins with known enzyme activity, specifically molecular chaperones and LEA proteins (Thomashow, 1999;Maruyama et al., 2009).
In this study, we describe the first transcriptome profiling by combining RNA-Seq with DREB1/CBF and cor gene expression analysis of the model legume L. japonicus after subjecting it to cold stress. To the best of our knowledge, this is the first report of L. japonicus response under low temperature stress. The results of this study can therefore guide the development of new agricultural forage practices and serve as a template for future studies that expand this understanding to other species within the genus.

Plant Materials and Growth Conditions
Seeds from the L. japonicus ecotype Gifu B-129 were scarified by stirring in pure sulfuric acid for 3 min, washed 10 times with sterile distilled water, and then sowed in Petri dishes containing 0.8% agar solution. Plates were incubated for 7 days in a growth chamber under a 16/8 h photoperiod at 24 • C/21 • C ± 2 • C (day/night) and 55/65 ± 5% relative humidity. Light, at a photosynthetic flux density of 250 µmol m −2 s −1 , was provided by F 40W Grolux fluorescent tubes. Seedlings were transferred to sterilized sand-perlite (2:1), placed in the previously used growth chamber under the same conditions, and irrigated with half-strenght Hoagland solution (Hoagland and Arnon, 1950).
Plants with 4-6 fully developed leaves, corresponding to roughly 3 weeks of development, were used in all experiments.

Treatments and Experimental Design
The experimental design was completely randomized, with three biological repetitions per treatment. For both the control and stress tests, the 3 week-old seedlings were placed in a Percival E-30B (Percival Scientific, Perry, IA, USA) growth chamber under a 16/8 h photoperiod (day/night). Lighting conditions were identical to those used in Section Plant Materials and Growth Conditions Plants under cold stress treatment were kept at 9 • C during daylight simulation and 5 • C at night, whereas the controls were kept at conditions identical to Section Plant Materials and Growth Conditions.
For RNA-Seq, shoots were harvested after 1 complete 24 h cycle, immediately frozen in liquid N 2 , and stored at −80 • C until their RNA was completely extracted. For DREB1/CBF and cor gene expression analysis, shoots were harvested after 0, 1, 3, 8, and 24 h of cold stress treatment.

RNA Extraction and Expression Analysis
Total RNA was extracted using a Plant Spectrum Total RNA Kit (Sigma) according to the manufacturer's instructions. RNA was checked for quality and quantified using agarose gel electrophoresis and spectrophotometric analysis. Total RNA samples (5-10 µg/sample) were shipped on dry ice to the Instituto de Agrobiotecnologia de Rosario (Rosario, Argentina) for RNA-Seq and analysis. mRNA was purified using oligo-Dt and cDNA was synthesized. Samples were prepared for RNA sequencing as described in the Illumina TruSeq R RNA Sample Preparation Guide (July 2012). High-performance, paired-end (2 × 100 bp) sequencing was performed on an Illumina Hiseq 1500. Low-quality RNA-Seq reads (Qscore < Q30) were analyzed using FastQC (Version 0.11.2) and discarded (Andrews, 2010). Then, a total of 174520000 reads were aligned against the L. japonicus genome (Gene Model, Release 2.5; Sato et al., 2008) using TopHat (Version 2.0.12; Trapnell et al., 2012).
To identify differentially expressed genes, a pair-wise comparison between normalized gene expression values (FPKM) for both conditions was performed using a t-test at the 99.99% confidence level using Cufflinks (Version 2.2.1; Trapnell et al., 2012). Illumina reads generated from all 6 samples are available at the NCBI BioProject browser, accession number PRJNA288510, BioSample accessions SAMN03801565, SAMN03801566, SAMN03801567, SAMN03801568, SAMN03801569, and SAMN03801570.
The overall RNA extraction procedure was identical for relative quantification of DREB1/CBF and cor genes. The absence of DNA from the RNA samples was tested by null polymerase chain reaction (PCR) amplification of the universal rDNA primer pair ITS1/ITS4 . Then, cDNA was synthesized from 3 µg of total RNA using M-MLV Reverse Transcriptase (Promega) and 100 pmol of the oligo-Dt primers, as per the supplier's instructions.

RNA-Seq Transcript Analysis
An MA-plot was constructed using R (Version 3.1.3) for Windows, which contained both the CummeRbund and ggplot2 packages (Team, 2014). All transcripts were used to make the plot.
Gene Ontology (GO) enrichment analysis was carried out to reveal the biological processes differentially expressed under the stress conditions, using Fisher's exact test against the L. japonicus entry in the legumeIP database (Al-Shahrour et al., 2004;Li et al., 2012). Blast2GO software (https://www.blast2go.com/) was used and applied only for those transcripts that showed significant differential expression, defined, with respect to fold change (FC), as those for which log 2 FC ≥ |2| (Conesa et al., 2005).
Identification of the Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways modified under low temperature treatment was used to corroborate the GO enrichment analysis (Kanehisa et al., 2010). KEGG annotation was done with the KEGG Automatic Annotation Server (http:// www.genome.jp/kaas-bin/kaas_main), again only for those transcripts for which log 2 FC ≥ |2| (Moriya et al., 2007). Classification was also performed using the LegumeIP database .
Up-regulated TFs were identified and classified through PlantTFDB (http://planttfdb.cbi.pku.edu.cn/) and the L. japonicus entry in the Legume IP database Jin et al., 2014 Kanehisa et al., 2010;Li et al., 2012).Transcripts not functionally annotated by one or more of the above methods were assigned functional annotations based on a BLASTX search vs. the Genbank NR protein database for E-values less than 10 −3 ; the lowest E-value sequences were chosen as representative.
To validate the results obtained from RNA-Seq, the expression data of 10 randomly chosen genes were analyzed by quantitative real time PCR (qRT-PCR; Supplementary Table 1)

qRT-PCR
To validate the RNA-Seq analyses and to quantify the relative expression of DREB1/CBF and cor genes, primers were designed with the help of Primer3 software (Untergasser et al., 2012). The primer pairs are given in Supplementary Table 1. An aliquot of 5 µL of 1:8 diluted cDNA was used in the qRT-PCR reactions, made using 15 µL of the FastStart Universal SYBR-Green Master Mix (Rox, Roche) and 2.5 pmol of each primer. Three biological replicates, every one accompanied by two technical replicates, were performed for each sample and gene. Cycling parameters were used that consisted of an initial step of 95 • C for 10 min and a two-step cycle of 95 • C for 30 s and 60 • C for 1 min, repeated 40 times. This was followed by the dissociation protocol. Amplifications were performed on an Mx3005P qPCR apparatus (Stratagene, CA, USA). The average threshold cycle was determined for each transcript. Gene quantification was determined based on the relative expression of the target gene vs. the reference gene EF-1α . For comparative purposes, average threshold cycle values of the control samples were used as reference.
The statistical analysis of relative gene expression was performed using the InfoStat/L program (Di Rienzo et al., 2008).

Sequencing and Read Assembly
Almost 68% of the transcripts resulting from read alignment mapped with the reference genome, allowing for the identification of 40993 genes (Figure 1). In total, 9781, or approximately 24%, were differentially expressed between treatments (p < 0.01, represented as red points in Figure 1). Of those genes, 1077, or about 11%, showed more than a two-fold change in their expression ratio. Within this group, 713 genes were up-regulated under cold stress treatment, while 364 were down-regulated; 133 and 75 of these genes, respectively, have no annotation in the L. japonicus genome. These 1077 genes were used for functional classification and annotation. It is worth noting that the percentage of mapped reads obtained was expected due to unavailability of complete genome information and the presence of highly repetitive sequences in the L. japonicus genome (Sato and Andersen, 2014). down-regulated genes, as defined by those for which log 2 FC≥ 2 and log 2 FC ≥ −2, respectively. GO enrichment values (represented as percentages) were calculated by generating the ratios for the number of observed over expected GO term repeats.

GO Enrichment Analysis
GO terms were successfully identified for 87.2% of the transcripts. The ten GO terms most highly represented in the upand down-regulated gene list, with false discovery rate p-values smaller than or equal to 0.01, were identified (Figures 2A,B, respectively).
The main up-regulated GO term was phospholipid transport, followed by features related to plasma membrane composition, specifically with respect to the anchored component of the membrane and response to endoplasmic reticulum stress. In addition, six of the ten up-regulated GO terms were related to the cell wall, both specifically to the cell wall and to the glucuronoxylan, xylan, hemicellulose, cell wall polysaccharide, and cell wall macromolecule metabolic processes. The gluthatione transferase activity GO term, related generally to glutathione metabolism, was also highly up-regulated. Meanwhile, the light harvesting complex, light harvesting, and photosynthesis GO terms were significantly down-regulated.

KEGG Metabolic Pathway Identification
A total of 579 (66.6%) of the annotated genes matched a KEGG pathway, with the highest percentages pertaining to the metabolism, organismal system, and human diseases categories for both the up-and down-regulated genes ( Figures 3A,B, respectively). After metabolism, organismal system was the most represented category in the KEGG analysis (Figure 3). Many of the up-regulated genes in this category were classified as  Table  2). Unexpectedly, the least represented category in both cases was that of environmental information processing. However, some TFs were identified in the environmental adaptation pathway (Supplementary Table 2); almost all of them were up-regulated.
The different genes were then further classified within these categories (Figure 4); it is worth noting that within the organismal system classification, environmental adaptation was most heavily represented, with 44 genes significantly affected. Overall, the KEGG results support conclusions drawn from the GO data. Glutathione S-transferase genes were up-regulated (Table 1), as were most of those for starch and sucrose metabolism. In addition, the transcription of some of the genes involved in arginine and proline metabolism were upregulated by more than a factor of two. Meanwhile, lightharvesting complex proteins were down-regulated, as were the chloroplast and thylakoid membrane terms, including those of the plastoglobule, plastid thylakoid, and chloroplast thylakoid. Finally, cutin, suberine, and wax biosyntheses in the lipid metabolism category were down-regulated.

TF Identification and Classification
The expression levels of 41 plant TFs were up-regulated by low temperature. These TFs were sorted into 14 different families ( Figure 5), with seven of these classified as MYB TFs; one of these (transcript XLOC_019119) was one of the two genes that had not been previously annotated (Supplementary Table  3). Additionally, five of the up-regulated TFs belonged to the NAC family, while six were part of the AP2/ERF family. One calmodulin binding transcription activator (CAMTA) TF (chr4.CM0307.390.r2.d) was found up-regulated, along with a number of calmodulin genes (Supplementary Tables 2, 3).

Differentially Expressed Genes with no Previous Annotation
Overall, 50 (37.6%) and 24 (32.0%) of the up-and downregulated genes, respectively, did not hit any sequence in the nonredundant protein GenBank database. Meanwhile, 27 (20.3%) and 20 (26.7%) of the up-and down-regulated genes, respectively, hit for an E-value greater than 10 −3 . The rest of the transcripts matched a sequence with an E-value lower than 10 −3 ; however, from those, 22 (39.3%) and 10 (24.4%) of the up-and downregulated genes, respectively, hit to an unknown protein. In all, 34 (60.7%) and 21 (75.6%) of the remaining up-and down-regulated transcripts, respectively, matched with a characterized sequence in the GenBank database (Supplementary Tables 4, 5).
Aside from the two previously discussed, non-annotated transcripts that were identified as TFs (Supplementary Table  3), other up-regulated sequences gave transcripts that hit with ribonucleases, retrotransposon proteins, serine/threonineprotein kinases, a sucrose-phosphate synthase 4, fatty acid protein elongation, and a cellulose synthase-like protein (Supplementary Table 4). Meanwhile, the down-regulated genes contained transcripts that matched with the ribonuclease H family, a serine/threonine protein kinase, a mitochondrial arginine transporter, and a nuclear pore complex protein (Supplementary Table 5).

DREB1/CBF Gene Expression Profiling
After 24 h of cold stress, L. japonicus showed no expression differences in this family of regulatory genes, though their response may be evident at shorter times. For this reason, relative expression profiles were studied.

DISCUSSION
Because metabolism was the most expressed category in the KEGG classification, we further analyzed these pathways. The organismal system category, the second most expressed, was also analyzed.

Lipid Metabolism
The regulation results help to account for observations related to lipid metabolism in previous studies. With respect to upregulation, earlier work had demonstrated molecular changes in the lipid membrane of plants undergoing cold and freezing stress (Iba, 2002;Welti et al., 2002) due to altered lipid biosynthesis and biomembrane rearrangement (Smolenska and Kuiper, 1977), as well as specific fatty acid changes (Willemot et al., 1977). Meanwhile, with respect to down-regulation, plants have been shown to decrease their cuticular wax content when subjected to cold stress (Shepherd and Wynne Griffiths, 2006), which is consistent not just with these results but with the negative correlation observed between frost tolerance and high epicuticular wax load in Salix sp. (Hietala et al., 1997).
Changes in the composition of wax due to stress must also be analyzed, as these features play a major role in stressor response; for example, two transgenic Arabidopsis plants with significantly enhanced drought tolerance were shown to differ in both wax composition and differential freezing response (Zhang et al., 2005(Zhang et al., , 2007.

Cell Wall Metabolism
Genes related to the cell wall are regulated by many different abiotic stresses, particularly water deficit stress and wounding; during severe low temperature conditions, extracellular ice formation could result in cell dehydration, collapse, or lysis due to ice crystal extension (Gall et al., 2015). This explains the strong up-regulation of genes involved in cell wall modification during low temperature stress, both in the literature (Domon et al., 2013;Dong et al., 2013) and in this study.

Xenobiotic Biodegradation
Many authors have highlighted the important role of diverse gluthatione-dependent enzymes in detoxification of xenobiotics and ROS, as well as in the adaptation to different abiotic stresses, particularly cold stress (Dixon et al., 1998;Anderson and Davis, 2004;Halusková et al., 2009). During cooling, H 2 O 2 and other ROS may accumulate, causing oxidative stress and cellular damage. In these conditions, reduced glutathione could first be oxidized by ROS, then regenerated by the activity of NADPHdependent gluthatione reductase [EC:1.6.4.2]. Unsurprisingly, enhanced glutathione content and glutathione reductase activity have been correlated to cold stress in different plant species (Kocsy et al., 2001). In fact, overexpressing gluthatione Stransferase [EC:2.5.1.18] in transgenic rice enhances germination and growth at low temperature (Takesawa et al., 2002). Therefore, the up-regulation in glutathione transerase activity in this study comes as no surprise.

Photosynthesis and CO 2 Fixation
The down-regulation of terms related to photosynthesis and CO 2 fixation are consistent with previous studies. Low temperature conditions are known to cause a reduction in maximum quantum yields for CO 2 uptake, losses in the photochemical efficiency of photosystem II, and, with prolonged exposure to excessive light, a decreased rate of light saturated photosynthesis (Kratsch and Wise, 2000;Renaut et al., 2005). This phenomenon, photoinhibition, is associated with photodamage in plant cells and is usually used as a marker for cold tolerance (Demmig-Adams and Adams Iii, 1992;Huner et al., 1993;Long et al., 1994).
Some authors describe photosynthesis as a sensor that is capable of regulating the imbalance between energy uptake and CO 2 fixation (Huner et al., 1993). From this perspective, photoinhibition should be the consequence of a plant's capacity to adjust photosynthetically to the prevailing environmental conditions, rather than the result of damage or injury (Huner et al., 1993;Ensminger et al., 2006). However, photodamage still should not be discounted.
Proteomics analysis of Thellungiella halophila at low temperature demonstrated that half of the identified coldresponsive proteins were related to chloroplast physiology and function (Gao et al., 2009), suggesting at least partial regulation of cold stress tolerance through chloroplast metabolism. A comparable analysis of rice gave similar results (Hashimoto and Komatsu, 2007), further corroborating our data.

Starch Metabolism
Several studies have demonstrated a strong correlation between sugar concentration and cold stress tolerance (Tarkowski and Van den Ende, 2015). In fact, some key enzymes of this metabolism, specifically sucrose synthase [EC:2.4.1.13] and sucrose phosphate synthase [EC:2.4.1.14], can be regulated by low temperatures (Guy et al., 1992;Sasaki et al., 2001); however, this only applies to the first of these enzymes in this study ( Table 1).
The role of starch and sucrose metabolism in this context has been reinforced in multiple mutant and transgenic plant studies. A mutant plant at the eskimo1 (esk1) locus, which exhibits higher levels of soluble sugars, was frost tolerant (Xin, 1998). On the other hand, an Arabidopsis mutant (sfr4) that was deficient in tolerating freezing conditions after cold acclimation did not present cold-induced elevation of sucrose and glucose levels (Mckown et al., 1996). Meanwhile, Gilmour et al. (2000) studied the overexpression of the CBF3 TF in Arabidopsis, which confers cold tolerance, and found that total sugar levels were higher than those in control plants.
It could be argued that sucrose is a storage carbohydrate that can be easily catabolized when needed, making it useful either under stress or after leaving stress conditions (Pollock and Lloyd, 1987;Guy et al., 1992). Furthermore, the increase in these metabolite concentrations could be a way of maintaining carbon flux by compensating for CO 2 assimilation and avoiding photoinhibition (Savitch et al., 2000). On the other hand, this upregulation could also be a way to cope with osmotic impairment generated by low temperature conditions (Anchordoguy et al., 1987;Palonen and Junttila, 1999) In such cases, plants can respond through an osmotic adjustment (Xiong and Zhu, 2002;Beck et al., 2007), in a process for which sugars could be a key component (Guy et al., 1992).

Arginine and Proline Metabolism
Proline is a major organic osmolyte that accumulates in a variety of plant species in response to environmental stresses, such as drought, salinity, and extreme temperatures (Hare and Cress, 1997;Ashraf and Foolad, 2007). In plants, L-proline is synthesized from L-glutamic acid via 1 -pyrroline-5-carboxylate (P5C) through P5C synthetase (P5CS) and P5C reductase. Some authors have described a transcriptional regulation of the P5CS gene, the rate-limiting factor in proline biosynthesis, under osmotic stress (Yoshiba et al., 1995(Yoshiba et al., , 1997. Yoshiba et al. (1995) did not find any induction of P5CS during cold stress of Arabidopsis thaliana. However, this contrasts with the up-regulation of this gene observed in this study, combined with an increase in proline content in different L. japonicus accessions after 7 days of cold stress (data not shown). Differences in P5CS gene transcriptional regulation are to be expected, due to differences in plant species, genotype, timing, and intensity of stress imposition.
Although the actual roles of proline in plant osmotolerance remain controversial, there is supporting evidence for a positive effect on enzyme and membrane integrity, osmotic adjustment, and free radical scavenging (Kishor et al., 2005). Still, some authors have argued that proline accumulation under stress is a product of, and not an adaptive response to stress (De Lacerda et al., 2003;Maiale et al., 2004). This uncertainty calls for further study of this feature, at least in the case of L. japonicus.

Biosynthesis of Secondary Metabolites
Most phenolic compounds in plants are derived from the phenylpropanoid pathway and have many different physiological roles. Increases in phenolic compound content under abiotic stress, particularly with respect to flavonoids, have been extensively described (Christie et al., 1994;Dixon and Paiva, 1995), making our results consistent with previous observations. Among the enzymes involved in the phenylpropanoid biosynthesis pathway, phenylalanine ammonia-lyase [EC:4.3.1.24] is one of the most relevant (Rivero et al., 2001). This enzyme catalyzes the transformation of L-phenylalanine into trans-cinnamic acid, and has been shown to increase in activity in response to thermal stress (Leyva et al., 1995). However, we did not find any transcriptional regulation of this particular enzyme. Nevertheless, post-transcriptional regulation is possible.
Meanwhile, phenols are oxidized by peroxidases, which increase in activity in response to different types of stress (Jansen et al., 2001;Michalak, 2006). Most flavonoids outperform wellknown antioxidants, such as ascorbate, but different flavonoid oxidation processes could involve ROS scavenging as well (Hernández et al., 2009). Thus, peroxidase activity may play an essential role in these phenomena (Yamasaki et al., 1997).
Chalcone synthase [EC:2.3.1.74] and chalcone isomerase [EC:5.5.1.6] are codified by the CHS and CHI genes, respectively, and are key components of the flavonoid biosynthetic process; as such, their transcription is enhanced under cold stress (Wu et al., 2014), as was observed in this case.

BR Biosynthesis and Steroid Hormone Biosynthesis
BRs are steroidal plant hormones implicated in the promotion of plant growth and development, as well as in the stress response (Sasse, 1997). In fact, evidence suggests that BRs have a protecting effect at low temperatures, in that they activate coldstress related genes COR47 and COR78 (Müssig et al., 2002;Janeczko et al., 2007;Bajguz and Hayat, 2009). However, recent studies corroborate this one in suggesting that BR biosynthesis is down-regulated at low temperatures in certain plants (Gray and Heath, 2005;Hannah et al., 2005;Kaplan et al., 2007).
This response, seemingly the opposite of what would be beneficial, was reversed by Kim et al. (2010), who studied Arabidopsis with a BR-insensitive 1 (bri1) mutation. These plants suffered from defective BR signaling, which seemingly attributed to improve tolerance to cold stress.
The role of BR in the abiotic stress response is still not clear, but our results support a possible physiological relevance of this down-regulation.

Organismal System Category
Heat-shock proteins are molecular chaperones that cope with stress-induced denaturation of other proteins (Feder and Hofmann, 1999;Sun et al., 2002;Wang et al., 2004). Indeed, Sabehat et al. (1998) showed that an induced expression of two heat-shock proteins in Lycopersicon esculentum L. correlated with protection against some of the symptoms of chilling injuries.
Meanwhile, interleukin-1 receptor-associated kinase 4 [EC:2.7.11.1] (K04733) is a central part of the immune and inflammatory response in mammals (Janssens and Beyaert, 2003). The signaling domains of this kinase share significant homology with some of the plant proteins that mediate activation of mitogen-activated protein kinase (MAPK) pathways (Neill and Greene, 1998). These plant proteins have been related to pathogen attack response (Dardick and Ronald, 2006), as well as abiotic stress tolerance (Lehti-Shiu et al., 2009;Rodriguez et al., 2010). Particularly, Zou et al. (2006) showed an increased expression of a receptor-associated kinase of maize (ZmPti1) at low temperatures. In our study, the up-regulation of some of the aforementioned genes may reflect the induction of many mitogen-activated protein kinase signaling pathways as a consequence of cold stress imposition. However, the specific gene functions remain to be further studied.

Up-Regulated TFs
Genes that are members of the AP2/ERF family are plant-specific TFs that share a conserved DNA-binding domain and participate in the abiotic stress response via specific binding to DRE/CRT, a cis-acting element in the promoter's target genes . AP2/ERF can be divided into four different subfamilies, namely AP2, RAV, ERF, and DREB . Some of the most described AP2/ERF TFs in the cold stress response network belong to the DREB1/CBF (A-1) subgroup (Yamaguchi-Shinozaki and Shinozaki, 2005;Nakashima and Yamaguchi-Shinozaki, 2006). Among these, CBF3, CBF1, and CBF2 are highly similar in amino acid sequence and are quickly and transiently induced by cold stress. Their expression is ABA independent and their products activate multiple stressinducible target genes . Interestingly, none of the AP2/ERF TFs that are up-regulated were DREB1/CBFs (Supplementary Table 3). Different authors have described characteristic expression profiles of these TFs within the first 8 h of stress (Novillo et al., 2004;Vogel et al., 2005), indicating that their expression should be studied at shorter intervals. Expression levels of some of these measured genes are analyzed in this way in the following section.
On the other hand, many studies demonstrate participation of the MYB TFs in different cold stress response pathways (Chinnusamy et al., 2007). Particularly, some are involved in the regulation of CBF TFs through the ABA independent signaling pathway (Chinnusamy et al., 2007). In fact, MYB15 negatively regulates the expression of the CBFs, as has been shown in a mutant and transgenic Arabidopsis with impaired MYB15 expression (Agarwal et al., 2006). Some ABA-inducible MYB proteins may participate cooperatively in the ABA-dependent expression of different cold regulated genes as well (Shinozaki and Yamaguchishinozaki, 2000;Zhu et al., 2005;Dai et al., 2007). Because these MYB proteins are synthesized after endogenous levels of ABA accumulate, they do not play a signaling role until the later stages of the stress response. As a consequence, MYB TFs may participate in both the short and long term response to cold stress, and in an ABA-independent and/or dependent pathway.
Most of the other up-regulated TF families, such as WRKY, NAC, C2H2, and HSF, participate in different abiotic stresses responses (Chinnusamy et al., 2007;Chen et al., 2012;Scharf et al., 2012). Particularly, the NAC family may contribute to ABAdependent gene expression under various stresses, including cold (Tran et al., 2004;Yamaguchi-Shinozaki and Shinozaki, 2005;Nakashima et al., 2007Nakashima et al., , 2012. Meanwhile, the calcium-calmodulin-CAMTA complex participates in the transcription regulation of different target genes. In this sense, cellular calcium(II) levels act as a secondary messenger that modulates diverse physiological processes that are important for stress adaptation. In particular, a rapid calcium influx is required for proper cold acclimation (Doherty et al., 2009).
In addition, a possible function of the calcium-calmodulin-CAMTA complex was suggested in the DREB1/CBF signaling pathway (Thomashow, 2010). In fact, Doherty et al. (2009) proved that CBF1, CBF2, and ZAT12 were direct targets of CAMTA3 (Doherty et al., 2009). The evidence that CAMTA TFs and calcium signaling are upstream of the DREB1/CBF regulon make them relevant in the first few hours of stress imposition. However, this does not mean that calcium signaling does not play a role at longer timeframes. Thus, an induction in chr4.CM0307.390.r2.d gene expression may imply that calcium participates in the cold stress response in a DREB1/CBF independent pathway or downstream from their target genes, a hypothesis that is reinforced by our calmodulin gene expression data.
The DREB1/CBF response The DREB1/CBF regulatory genes are an AP2/ERF TF subgroup that binds to the DRE/CRT regulatory element present in the promoters of target genes (Yamaguchi- Shinozaki and Shinozaki, 2005;Mizoi et al., 2012). These genes have been well characterized as key components of the ABA-independent cold stress response, mainly in Arabidopsis (Thomashow, 2010;Medina et al., 2011). In fact, a distinctive CBF pattern expression has been described for within the first 8 h (Novillo et al., 2004;Vogel et al., 2005). Remarkably, until now, no CBF expression studies have been performed for the Lotus genus.
The CBF1 and CBF3 results of this study are in agreement with those obtained by Novillo et al. (2004) and Cook et al. (2004), which showed maximum expression levels within the first 2 h of cold stress in Arabidopsis plants. Different authors have also described a similar expression pattern for CBF2, reaching a maximum only several hours later than CBF1 and CBF3 (Cook et al., 2004;Novillo et al., 2004Novillo et al., , 2007. The observed differences in our study may be explained by the fact that most pertinent studies were completed using Arabidopsis, with only a few using legume species like Medicago spp. (Pennycooke et al., 2008;Zhang et al., 2011). In fact, there is no conclusive evidence for this expression pattern in L. japonicus. Interestingly, CBF2 has been described as a negative regulator of CBF1 and CBF3 (Novillo et al., 2004;Medina et al., 2011). As a consequence, the initial reduction observed in its expression level in L. japonicus might actually support the increase in CBF1 and CBF3.
Like CBF2, ZAT12 has been described as a negative regulator of DREB1/CBFs, though it is generally induced in this timeframe in Arabidopsis (Vogel et al., 2005;Zhai et al., 2010). However, there is no evidence of this TF behavior in legumes.
On the other hand, one TF that positively regulates DREB1/CBF TF expression is the MYC-like bHLH protein ICE1 (Yamaguchi- . The relevance of this TF in cold stress response has been demonstrated by ice1 Arabidopsis mutants in which DREB1/CBFs expression was deregulated, leading to a reduced expression of many of their downstream cold-responsive genes (Chinnusamy et al., 2003;Lee et al., 2005). However, it was shown that a cold induced posttranslational modification of ICE1, rather than an increase in expression, allowed for the positive regulation of the DREB1/CBF regulon (Sharma et al., 2005). Further studies should be carried out to clarify the role of ICE1 in the cold stress response of L. japonicus.
Finally, the expression patterns of RD29A and COR47 are consistent with the fact that these are known cor genes that drive the cold stress response in plants several hours after stress conditions begin (Thomashow et al., 2001).

CONCLUSION
Our study provides novel insights into the molecular mechanisms of L. japonicus cold stress response.Overall, 1077 differentially expressed genes were obtained using RNA-Seq data, consisting of 713 that were up-regulated and 364 that were down-regulated. These genes were classified by their functional annotation through BLASTx, GO, and KEGG analysis.
Cold stress in L. japonicus modulates gene expression of lipid and cell wall metabolism, together with transcription regulation of genes related to flavonoid and phenylpropanoid biosynthesis. Genes related to xenobiotic degradation, heatshock proteins, and starch and proline metabolism were up-regulated. In contrast, photosynthesis and chloroplast related genes were down-regulated, demonstrating a dependency of cold acclimation on the photosynthetic process.
Different types of TFs were up-regulated as well, with the MYB, AP2/ERF, and NAC families being the most numerous. Interestingly, two putative novel TFs were identified that participate in cold response. Finally, the expression profiles of some DREB1/CBFs and cor genes were described for the first 8 h of stress. To our knowledge, this is the first report of DREB1/CBF expression under low temperature conditions in this genus.
Our results constitute the first transcriptome profiling of the model legume L. japonicus under cold stress. Data obtained allowed for the identification of possible target metabolisms, information that could be used to improve stress tolerance in L. japonicus, not to mention other members of the genus, in the future.

AUTHOR CONTRIBUTIONS
PC, SM, OR, and FE conceived and designed the experiments. PC performed the experiments. PC and FE analyzed the data. OR contributed reagents, materials, and analytical tools. PC, SM, OR, and FE wrote the paper. All authors read and approved the final manuscript.

AUTHOR INFORMATION
PC is a CONICET fellow. SM, OR, and FE are career research members of CONICET.

ACKNOWLEDGMENTS
This work was supported by the following grants: PICT of Agencia Nacional de Promoción Científica y Tecnológica (ANPCYT, Argentina), Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina), and Universidad de San Martín (UNSAM). 3PA. We would also like to thank Daniel Davis for English language and general editing of the manuscript.