Identification of Genes Differentially Expressed in Response to Cold in Pisum sativum Using RNA Sequencing Analyses

Low temperature stress affects growth and development in pea (Pisum sativum L.) and decreases yield. In this study, RNA sequencing time series analyses performed on lines, Champagne frost-tolerant and Térèse frost-sensitive, during a low temperature treatment versus a control condition, led us to identify 4981 differentially expressed genes. Thanks to our experimental design and statistical analyses, we were able to classify these genes into three sets. The first one was composed of 2487 genes that could be related to the constitutive differences between the two lines and were not regulated during cold treatment. The second gathered 1403 genes that could be related to the chilling response. The third set contained 1091 genes, including genes that could be related to freezing tolerance. The identification of differentially expressed genes related to cold, oxidative stress, and dehydration responses, including some transcription factors and kinases, confirmed the soundness of our analyses. In addition, we identified about one hundred genes, whose expression has not yet been linked to cold stress. Overall, our findings showed that both lines have different characteristics for their cold response (chilling response and/or freezing tolerance), as more than 90% of differentially expressed genes were specific to each of them.


Introduction
Cold stress is one of the most important factors that limit plant productivity around the world. Understanding the molecular bases of the cold response is thus essential to breed cold-tolerant varieties. To survive winter frosts, plants need to acquire frost tolerance which depends on the duration and time of the exposition to low temperatures [1] and varies according to species [2,3] and genotypes [4]. Plants can adopt two strategies to overcome frost. The first one consists of escaping the main frost periods, which can be obtained by different adaptive mechanisms in the natural population like developing a rosette form or reducing the hypocotyl length, which has been well documented in Arabidopsis thaliana [5][6][7]. The second one concerns the acquisition of freezing tolerance (FT), through a phenomenon called cold acclimation [8]. Following exposure to low temperatures, plants increase their ability to tolerate cold temperatures. Cold acclimation reveals two mechanisms of tolerance, which include: chilling tolerance and the induction of freezing tolerance (FT). Chilling tolerance represents the ability of a plant to respond to low but positive temperatures inferior to 15 • C, and FT is an induced response where plants acquire an increased freezing tolerance upon a prior low but non-freezing Three samplings were performed in each experiment (LT and N). For the N experiment, the samplings were performed at the 20th (T0), 22nd (T1), and 27th (T2) day of the nursery period. For the LT experiment the samplings were realized at the 20th day (T0) of the nursery period and at the 3rd (T1) and 16th (T2) day of the LT period ( Figure 1). The developmental stage of the plants was regularly checked by the state of leaf unfolding, allowing the sampling of plants at the same developmental stage in both LT and N conditions, with the same final number of degree-days [30] despite differences in number of days until sampling. At the end of the recovery period, we confirmed the phenotype of the two lines, because the Ch plants survived the frost period, unlike the Te ones.
A total of 24 samples were harvested, corresponding to 2 lines × 2 treatments (LT and N) × 3 sampling times (T0, T1, and T2) × 2 biological replicates. For each sample, aerial parts of 3 plants were pooled and were immediately frozen using liquid nitrogen and stored at −80 °C until RNA extraction.

RNA Extraction and High Throughput Sequencing
Total RNA was isolated with the Plant RNAeasy Mini kit (Qiagen) and quantified using a NanoDrop. The paired-end sequencing was performed using an Illumina HiSeq2000 sequencer at the Next Generation Sequencing (NGS) platform of GeT: http://get.genotoul.fr/. About 900 million raw reads generated from the 24 samples were deposited to the National Center for Biotechnology Information (NCBI) in the SRA (Short Read Archive) database (BioProject #PRJNA543764).
Gene ontology (GO) (http://geneontology.org/ Gene Ontology Consortium and http://pantherdb.org/), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.kegg.jp/kegg/tool/map_pathway2.html KEGG Mapper) information was assigned based on the A. thaliana homologous proteins, and GO functional classification was performed (https://www.arabidopsis.org/tools/bulk/go/index.jsp). Given the importance of protein kinases [38][39][40][41][42], transcription factors and transcription regulators (TF) [43][44][45][46] in cold stress responses, a specific Three samplings were performed in each experiment (LT and N). For the N experiment, the samplings were performed at the 20th (T0), 22nd (T1), and 27th (T2) day of the nursery period. For the LT experiment the samplings were realized at the 20th day (T0) of the nursery period and at the 3rd (T1) and 16th (T2) day of the LT period ( Figure 1). The developmental stage of the plants was regularly checked by the state of leaf unfolding, allowing the sampling of plants at the same developmental stage in both LT and N conditions, with the same final number of degree-days [30] despite differences in number of days until sampling. At the end of the recovery period, we confirmed the phenotype of the two lines, because the Ch plants survived the frost period, unlike the Te ones.
A total of 24 samples were harvested, corresponding to 2 lines × 2 treatments (LT and N) × 3 sampling times (T0, T1, and T2) × 2 biological replicates. For each sample, aerial parts of 3 plants were pooled and were immediately frozen using liquid nitrogen and stored at −80 • C until RNA extraction.

RNA Extraction and High Throughput Sequencing
Total RNA was isolated with the Plant RNAeasy Mini kit (Qiagen) and quantified using a NanoDrop. The paired-end sequencing was performed using an Illumina HiSeq2000 sequencer at the Next Generation Sequencing (NGS) platform of GeT: http://get.genotoul.fr/. About 900 million raw reads generated from the 24 samples were deposited to the National Center for Biotechnology Information (NCBI) in the SRA (Short Read Archive) database (BioProject #PRJNA543764).

Differential Expression Analysis and Statistical Tools
For the differential expression analysis, Illumina-cleaned reads from the 24 samples were pseudo-aligned on the de novo transcriptome assembly using Kallisto v0.43.1 [50]. Count data were analyzed using a multifactorial design (line (Ch and Te), treatment (LT and N), and time (T0, T1, and T2)) with the R package DESeq2 [51]. The multiple steps of statistical and clustering analyses led to three sets of distinct DEGs (Differentially Expressed Genes) ( Figure S1). Firstly, only unigenes fulfilling the criterion TRUE (when the number of reads is sufficient to perform the statistical test) in DESeq2 were conserved for the next analyses. Then, unigenes represented by less than or equal to 48 normalized counts among the 24 samples were discarded. The p-values were calculated for nine combinations (ChNT0:TeNT0, ChNT0:ChNT1, ChNT0:ChNT2, ChLT0:ChLT1, ChLT0:ChLT2, TeNT0:TeNT1, TeNT0:TeNT2, TeLT0:TeLT1, and TeLT0:TeLT2). Only unigenes with an adjusted p-value ≤ 0.02 in at least one of the nine combinations were conserved. A sorting was then realized according to the significance of the functional annotation (E-value ≤ 9 × 10 −4 ). The resulting set of unigenes was submitted to three successive analyses of variance (ANOVA) using Multi Experiment Viewer [52] statistical tools, in order to identify unigenes related to the differences between the two lines or related to the cold responses. First, a two factor ANOVA was performed allowing us to retain only DGEs which varied according to the line and/or the treatment. Then a one-way ANOVA was realized, providing a set of unigenes related to the differences between the two lines (Line Response). The remaining genes were submitted to an additional one-way analysis allowing us to decipher genes responding identically to the LT in the two lines (Temperature Common Response: TCR) from those which presented different cold responses in Ch and Te (Temperature Line Specific Response: TLSR).
Afterwards, gene expression patterns were built and classified with the MeV clustering tools. After log2 transformation of normalized count and mean-centered reduced fit, hierarchical clustering (HCL) were performed using Pearson's correlation and average linkage clustering method. GO term enrichment analyses of the different sets of DEG were performed using AmiGO2 (http:/amigo1. geneontology.org/cgi-bin/amigo/term_enrichment).

Statistical Differentially Expressed Genes (DEGs) Corroboration
In order to support the expression of DEGs identified in this study, we used quantitative polymerase chain reaction (qPCR) data from an earlier published experiment [28] which presented only a few differences in the environmental conditions and the sampling times. The T0, T1, and T2 samples of the RNA-seq study were compared to the T0, T6, and T10 of the qPCR analysis, respectively ( Figure 1). Blastn searches (E-value ≤ 1 × 10 −25 , coverage ≥ 200 pb, ≥ 95% of identity) were performed to link the unigenes representing DEGs in the RNA-seq study to the ESTs used to define primers for the qPCR analysis.

Samples, Sequencing, and Assembly Assessment
Twenty-four RNA libraries were built from RNA samples extracted from the two contrasted pea lines, Ch (frost-tolerant) and Te (frost-sensitive), under LT (a low but non-freezing temperature regime is applied before submitting plants to frost) and N (frost is applied without any prior LT treatment). The Illumina sequencing of the 24 samples led us to obtain 886,477,626 paired-end reads. After removing low-quality sequences, a total of 793,583,651 clean reads were assembled into 150,342 contigs with a minimum length of 201 bp (Table S1). Following an analysis of the quality of the transcriptome assembly, 122,194 high-quality contigs, for a total of 118,787,279 bp and representing 73,225 unigenes were conserved for the further analyses. The largest contig was 11,608 bp long and the N50 value of the assembly was 1904 (Table S1). Blastx searches against A. thaliana protein sequences allowed to annotate 34% (24,854/73,225) of the unigenes (E-value ≤ 9 × 10 −4 ). The distribution of the biological processes of the annotated unigenes resembled those reported for A. thaliana genes, suggesting that the construction of the pea RNA-seq libraries did not induce an enrichment of sequences related to a particular class of function ( Figure S2).

Differential Expression Analysis and Clustering
The threshold applied on DESeq2 results (adjusted p-value ≤ 0.02) retained 11,076 unigenes, presenting a modulation of expression within the 24 different samples. Among them, 9676 unigenes were annotated using A. thaliana proteins. A two-way ANOVA performed on the 9676 annotated unigenes permitted us to detect 4981 genes, in which expression was significantly modulated according to the line and/or to the low temperature treatment ( Figure S1). Among them, 33 were compared with data previously obtained by qPCR. We observed that RNA-seq and qPCR data are correlated (R = 0.71) ( Figure 2). the N50 value of the assembly was 1904 (Table S1). Blastx searches against A. thaliana protein sequences allowed to annotate 34% (24,854/73,225) of the unigenes (E-value ≤ 9 × 10 -4 ). The distribution of the biological processes of the annotated unigenes resembled those reported for A. thaliana genes, suggesting that the construction of the pea RNA-seq libraries did not induce an enrichment of sequences related to a particular class of function ( Figure S2).

Differential Expression Analysis and Clustering
The threshold applied on DESeq2 results (adjusted p-value ≤ 0.02) retained 11,076 unigenes, presenting a modulation of expression within the 24 different samples. Among them, 9676 unigenes were annotated using A. thaliana proteins. A two-way ANOVA performed on the 9676 annotated unigenes permitted us to detect 4,981 genes, in which expression was significantly modulated according to the line and/or to the low temperature treatment ( Figure S1). Among them, 33 were compared with data previously obtained by qPCR. We observed that RNA-seq and qPCR data are correlated (R = 0.71) (Figure 2).

Figure 2.
Scatter plot between the two sets of transcripts issued from the present RNA-seq and qPCR data from Reference [28]. The x-axis corresponds to log2 of qPCR ratios and the y-axis shows log2 of RNA-seq ratios for 33 transcripts. The linear relationship between the 2 variables and their correlation coefficient, R = ∑xy / √ ∑x 2 ∑y 2 are reported in the top of the graph.
The examination of the 33 expression patterns reflects the correlation between the data from the two methods of gene expression analysis and supports the robustness of the whole transcriptomic analysis ( Figure S3). Then, the 4981 DEGs were submitted to a one-way ANOVA based on the line factor, which led to identify 2487 genes differentially expressed between Ch and Te but not modulated during the cold treatment (Line Response, LR). This set represents a part of the constitutive differences between the two lines. The hierarchical clustering (HCL) analysis divided this set into two blocks. The first block contained 906 genes, which were more expressed in the Ch than in the Te samples (LR (a), Figure 3A (a)) and inversely, the second one was composed of 1581 genes which were more expressed in the Te than in the Ch samples (LR b, Figure 3A (b)). Then a second one-way ANOVA based on the cold treatment factor was performed with the remaining 2494 (4981 − 2487 = 2494) DEGs ( Figure S1) and revealed 1403 genes related to the "Temperature Common Scatter plot between the two sets of transcripts issued from the present RNA-seq and qPCR data from Reference [28]. The x-axis corresponds to log2 of qPCR ratios and the y-axis shows log2 of RNA-seq ratios for 33 transcripts. The linear relationship between the 2 variables and their correlation coefficient, R = xy/ √ x 2 y 2 are reported in the top of the graph.
The examination of the 33 expression patterns reflects the correlation between the data from the two methods of gene expression analysis and supports the robustness of the whole transcriptomic analysis ( Figure S3). Then, the 4981 DEGs were submitted to a one-way ANOVA based on the line factor, which led to identify 2487 genes differentially expressed between Ch and Te but not modulated during the cold treatment (Line Response, LR). This set represents a part of the constitutive differences between the two lines. The hierarchical clustering (HCL) analysis divided this set into two blocks. The first block contained 906 genes, which were more expressed in the Ch than in the Te samples (LR (a), Figure 3A (a)) and inversely, the second one was composed of 1581 genes which were more expressed in the Te than in the Ch samples (LR b, Figure 3A (b)). Then a second one-way ANOVA based on the cold treatment factor was performed with the remaining 2494 (4981 − 2487 = 2494) DEGs ( Figure S1) and revealed 1403 genes related to the "Temperature Common Response" (TCR) of the two lines. The HCL analysis carried out on these 1403 significant genes could separate 520 ( Figure 3B (a)) and 883 ( Figure 3B (b)) genes up-expressed in the N condition and the LT condition, respectively. The remaining 1091 genes (2494 − 1403 = 1091), which responded differently in Ch and Te under the LT treatment were considered to be associated to the "Temperature Line Specific Response" (TLSR). The HCL performed on this last set ( Figure 3C) showed four distinct expression patterns. The first subset (TLSR (a)) gathered 253 genes which were more expressed in Ch than in Te and down-regulated during the LT conditions in both lines ( Figure 3C (a)). The second subset (TLSR (b)) was composed of 228 genes which were also more expressed in Ch than in Te but up-regulated during the LT conditions ( Figure 3C (b)). The third one (TLSR (c)) contained 131 genes that were less expressed in Ch than in Te and down-regulated during the LT conditions ( Figure 3C (c)) and the last of the four subsets (TLSR (d)) was composed of 479 genes which were likewise less expressed in Ch than in Te but up-regulated during the LT treatment ( Figure 3C (d)). Response" (TCR) of the two lines. The HCL analysis carried out on these 1403 significant genes could separate 520 ( Figure 3B (a)) and 883 ( Figure 3B (b)) genes up-expressed in the N condition and the LT condition, respectively. The remaining 1091 genes (2494 − 1403 = 1091), which responded differently in Ch and Te under the LT treatment were considered to be associated to the "Temperature Line Specific Response" (TLSR). The HCL performed on this last set ( Figure 3C) showed four distinct expression patterns. The first subset (TLSR (a)) gathered 253 genes which were more expressed in Ch than in Te and down-regulated during the LT conditions in both lines ( Figure 3C (a)). The second subset (TLSR (b)) was composed of 228 genes which were also more expressed in Ch than in Te but up-regulated during the LT conditions ( Figure 3C (b)). The third one (TLSR (c)) contained 131 genes that were less expressed in Ch than in Te and down-regulated during the LT conditions ( Figure 3C (c)) and the last of the four subsets (TLSR (d)) was composed of 479 genes which were likewise less expressed in Ch than in Te but up-regulated during the LT treatment ( Figure 3C (d)).

Functional Annotation, Gene Ontology (GO) Term Enrichment, and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathways
For the three sets of DEGs described above (LR, TCR, and TLSR) and distributed into eight subsets, the GO term enrichment analysis is detailed in Figure S4, the mapping on KEGG pathways in Table S2, the list of genes annotated as coding putative kinases in Table S3, and those annotated as coding transcription factors in Figure 4. In addition, a summary of all these results is provided in Table 1 and the list of genes known to be related to cold stress in Table 2. This comprehensive analysis of the function of DEGs allowed us to present an overview of the differences between the two lines and to decipher their respective cold responses.
We observed that the five most significant enriched GO terms related to biological processes were the same for all of the eight subsets. These terms were "metabolic process", "cellular process", "localization", "cellular component organization or biogenesis", and "response to stimulus". This result seems to show that there is no difference in biological processes involved in the different subsets, at least not at this level of functional assignment. On the other hand, the enriched GO terms related to cellular localization were different within the subsets. In two cases, LR (a) and TCR (a), the five most significant terms related to localization concerned chloroplast and plastid components. This was also the case for a part of the terms enriched in the TLSR (a) subset. For all other subsets, the major enriched terms were related to cell, cytoplasm, and/or intercellular part.
The mapping of all DEGs on the KEGG pathways revealed relatively important differences in the number of mapped genes according to the subset (Table S2, Table 1). Indeed, 36% of the DEGs of the TLSR (d) subset were mapped on the pathways versus 80% of those of the TLSR (a) subset. This analysis, unlike the enrichment of GO terms, revealed a diversity of pathway representations between the subsets. For example, the most represented pathways in the LR (b) and TCR (b) subsets were related to the ribosome metabolism, while those of the LR (a) and TLSR (b) subsets were related to the amino acid metabolism.
The analysis, using iTak, allowed us to annotate 193 putative kinases among the DEGs of the three sets, including 96, 54, and 43 in the LR, TCR, and TLSR sets, respectively (Table S3, Table 1). Kinases were the most represented in the TLSR (b) and (d) subsets (4.8% for each) and the less represented in the TLSR (c) subset (1.5%). A total of 397 putative TFs was also annotated, with 188, 114, and 95 in the LR, TCR, and TLSR subsets, respectively ( Figure 4, Table 1). TFs were the most represented in the TLSR (d) subset (9.6%) and the less represented in the TCR (a) subset (4.8%). The most represented family differs according to the subset, bHLH and AP2-EREBP being the most represented in three subsets each (LR (b), TCR (a) and (b) for bHLH and LR (a), TLSR (b) and (c) for AP2-EREBP). We noted that some of the kinases and TFs are known to be related to cold stress, including Myb, CBF, and WRKY (Table 2 and see below).
The functional classification using GO knowledgebase allowed us to classify 145 genes in GO terms related to the cold response, with 128, 11, and 6 in "response to cold" GO:0009409, "cold acclimation" GO:0009631, and "cellular response to cold" GO:0070417, respectively ( Table 1, Table 2). Even if it was expected to identify genes related to these GO terms within the TCR (44 genes: 3.1%) and TLSR (36 genes: 3.3%) sets, it was more surprising to identify 65 genes (2.6%) in the LR set. This could mean that several genes that are known to be induced during cold stress in other species could be constitutively expressed in Ch or Te. The proportion of genes related to these GO terms in the different subsets was between 2.3% (LR (b)) and 6.9% (TLSR (c)).  Table 2. GO matching of the cold response (i.e., "response to cold", "cold acclimation", "cellular response to cold") in eight subsets of differentially expressed genes. The Ch (Champagne) up and Te (Térèse) up columns correspond to "LR", followed by two columns of "TCR" down and up and the last four columns referring to the "TLSR" a, b, c, and d, as described in the legend of Figure 3.

Discussion
The genetic studies previously carried out by our team have shown that the Champagne and Térèse lines differ in their response to cold. Following a period of acclimation, the Ch line becomes tolerant to frost while Te remains sensitive [29,126]. This present study aims to provide molecular elements to explain the cold responses in pea. The statistical treatments of the RNA-seq data of 73,225 unigenes led us to distinguish 4981 DEGs, divided into three main sets according to their expression patterns. The analysis of these genes, their modulation of expression, and their affiliation to metabolic pathways has allowed us to enrich our knowledge on the behavior of these two lines regarding their response to cold.

Differences in Gene Expression Between The Two Lines Independently of the LT Treatment: Predispositions to Face Cold Stress?
The class of genes considered here corresponds to genes which are differentially expressed in the two lines at T0, i.e., genes which are more expressed in Ch than in Te samples, or inversely. The opposition of gene expression patterns was observed throughout the time course study with an almost constant expression level in each line, independently of the LT treatment ( Figure 3A). We can therefore consider these genes as being part of the constitutive differences between the two lines. In both Ch and Te up-expressed subsets ( Figure 3A (a) and 3A (b), respectively), the KEGG pathway annotation revealed that a majority of these genes are related to metabolism (84.4% and 68.9% for the Ch and Te up-expressed subsets, respectively) suggesting important intrinsic differences between the two lines regarding cellular metabolism. In particular, 73 and 49 genes were related to amino-acid and energy metabolisms among genes up-expressed in Ch and Te, respectively. In addition, 102 genes in the Ch up-expressed subset and 39 in the Te up-expressed subset were related to RNA and protein metabolism. The GO-enrichment analyses were consistent with these results, highlighting the assignment of genes related to metabolism and to RNA and RNA surveillance pathways in the Ch and Te up-expressed subsets, respectively. Interestingly, even if the majority of the genes considered here did not present any modulation of expression during the time course study, 50/906 (5.5%) and 108/1581 (6.8%) of genes from the Ch and Te up-expressed subsets respectively, are known to be involved in responses to diverse stimuli, notably to responses to stress. In particular, within the Ch up-expressed subset, we identified genes coding a cold-regulated 413 inner membrane protein 1

Discussion
The genetic studies previously carried out by our team have shown that the Champagne and Térèse lines differ in their response to cold. Following a period of acclimation, the Ch line becomes tolerant to frost while Te remains sensitive [29,126]. This present study aims to provide molecular elements to explain the cold responses in pea. The statistical treatments of the RNA-seq data of 73,225 unigenes led us to distinguish 4981 DEGs, divided into three main sets according to their expression patterns. The analysis of these genes, their modulation of expression, and their affiliation to metabolic pathways has allowed us to enrich our knowledge on the behavior of these two lines regarding their response to cold.

Differences in Gene Expression Between The Two Lines Independently of the LT Treatment: Predispositions to Face Cold Stress?
The class of genes considered here corresponds to genes which are differentially expressed in the two lines at T0, i.e., genes which are more expressed in Ch than in Te samples, or inversely. The opposition of gene expression patterns was observed throughout the time course study with an almost constant expression level in each line, independently of the LT treatment ( Figure 3A). We can therefore consider these genes as being part of the constitutive differences between the two lines. In both Ch and Te up-expressed subsets ( Figure 3A (a) and (b), respectively), the KEGG pathway annotation revealed that a majority of these genes are related to metabolism (84.4% and 68.9% for the Ch and Te up-expressed subsets, respectively) suggesting important intrinsic differences between the two lines regarding cellular metabolism. In particular, 73 and 49 genes were related to amino-acid and energy metabolisms among genes up-expressed in Ch and Te, respectively. In addition, 102 genes in the Ch up-expressed subset and 39 in the Te up-expressed subset were related to RNA and protein metabolism. The GO-enrichment analyses were consistent with these results, highlighting the assignment of genes related to metabolism and to RNA and RNA surveillance pathways in the Ch and Te up-expressed subsets, respectively. Interestingly, even if the majority of the genes considered here did not present any modulation of expression during the time course study, 50/906 (5.5%) and 108/1581 (6.8%) of genes from the Ch and Te up-expressed subsets respectively, are known to be involved in responses to diverse stimuli, notably to responses to stress. In particular, within the Ch up-expressed subset, we identified genes coding a cold-regulated 413 inner membrane protein 1 (AT1G29395), and a 3-oxoacyl-synthase II (AT1G74960), which are known to accumulate in cold/freezing conditions in order to preserve chloroplast membranes integrity [110,111]. We also noted in the Ch up-expressed subset the presence of two genes, the first being a mediator of the RNA polymerase II transcriptional subunit 32 (AT1G11760) which regulates the expression of the second, a dehydration element B1A (CBF3, AT4G25480). Both are well known for their involvement in cold acclimation in A. thaliana [112]. More surprisingly, several genes are also known to be involved in the cold responses within the Te up-subset. For example, we found in this subset three genes coding AGAMOUS-like proteins (AT3G61120, AT4G22950, AT2G45660), known to lead to early flowering in A. thaliana and induced by an extended cold treatment [127]. We also identified genes coding a calmodulin-binding transcription activator protein (AT5G64220), which is involved in the rapid induction of CBF factors [53], and an adenine nucleotide alpha hydrolase-like super family (AT1G09740), which enhances freezing tolerance in Arabidopsis after a short period of cold-acclimation [54]. Moreover, genes coding a raffinose synthase (AT5G20250, [45]), a xyloglucan endotransglucosylas/hydrolase (AT4G03210 [55]), and an enolase (in pea [128] and LOS2, AT2G36530 [56]) were scored in the Te up-expressed subset, all of them being involved in cold stress in A. thaliana. The fact that some genes related to cold stress are up-expressed in Te could also be expected, since this line could possess deficient alleles that would not confer to Te an effective defense against cold.
Concerning the kinases, we observed within the Ch up-expressed subset one gene encoding a MPK3 (mitogen-activated protein kinase 3, AT3G45640) shown to be up-expressed in response to cold stress in A. thaliana [57]. Within the Te up-expressed subset, we identified another MAP3K (MAPK/ERK kinase kinase 1, AT4G08500) that plays a major role in cold stress signaling in A. thaliana [113]. More interestingly, 50 (5.52%) and 138 (8.73%) TFs were identified in Ch and Te up-expressed subsets, respectively. Among them, 50 distinct TF families were represented, including 19 TF families which gathered preferentially expressed genes in Ch and Te (i.e., common TF families) and 7 and 24 families which were represented by genes preferentially expressed in Ch or Te, respectively (i.e., line-specific TF families). Among them, several TFs were previously signaled as differentially expressed in cold stress in legume species [16,18,19,22] or are known to be involved in cold acclimation or freezing tolerance in A. thaliana, e.g., DREB1 (CBF3), CAMTA, NAC, and WRKY [25,129]. The other TFs, including ARID, BSD, mTERF, RWP-RK, S1Fa-like, SOH1, SRS, SWI/SNF-BAF60b, TAZ, and TUB are reported for the first time in pea and may be also involved in the intrinsic phenotypic differences between the two lines.
A recent study, using genetic structure and linkage disequilibrium in a large collection of pea germplasm, also highlighted the line-dependent differences [130]. These differences may be due to the fact that Ch belongs to the winter forage lines cohort and Te is issued from the spring lines set. Overall, our analysis reflects the importance of intrinsic differences between lines at the gene expression level [28,128,131]. Among this "Line Response" set, many genes coding proteins related to cold responses were identified. These results suggest that Ch and Te have constitutive defenses against cold stress, which have been described in the majority of cases as induced in other species and particularly in A. thaliana.

How Pea Faces Cold Stress
The identification of genes showing a significant increase/decrease of expression in Ch and Te during the LT and no fluctuation in the N experiment suggests that both lines have the capacity to undertake molecular modifications in response to cold stress, that can be considered as the chilling response shared by both genotypes. These genes were gathered into the TCR set. In addition, genes showing variation of expression only in Ch during the LT treatment could be attributed to the FT capacity of this frost tolerant line and were clustered within the TLSR set. Elsewhere, the expression pattern variations are more diverse for this set of TLSR, both at the beginning of the LT and also over time, as well as in N condition between these two lines. Furthermore, it is important to note that for both TCR and TLSR sets, we have scored more up-expressed genes in Te than in Ch.

Chilling Response
The TCR set was separated into two subsets, gathering almost the same number of genes which were down-or up-regulated during the LT treatment. Concerning down-expressed genes, it should be noted that the enrichment of GO terms related to chloroplast could provide evidence that genes involved in the photosynthetic system are hugely affected by chilling ( Figure S4 (F)). These observations are in agreement with the fact that cold stress leads to a disruption and/or dysfunction of photosynthesis and causes damages to thylakoid membrane and chloroplastic envelopes [132]. The up-expressed subset, for its part gathered, in particular, several RNA-binding proteins, which are involved in RNA and RNA surveillance pathways, operating in plant responses to abiotic stress [58,133]. Concerning the GO annotations, we observed that 30 (3.40%) and 14 (2.7%) genes of the up-and down-expressed subsets respectively, were associated to cold ( Table 2). For example, genes representing a STARCH EXCESS 1 (AT1G10760 [114]), a Glycine-rich RNA-binding protein7 (AT2G21660 [59]), two low temperature and salt-responsive protein LTI6A and LTI6B (AT3G05880 and AT3G05890 [60]), a Late Embryogenesis Abundant 4-5 (AT5G06760 [61]), and a Galactinol synthase 2 (AT1G56600 of the raffinose pathway [62]) were found in the up-expressed subset and are known to present an increase of expression in response to cold condition in A. thaliana. At last, a transcriptional adapter ADA2b (AT4G16420), which may repress freezing tolerance and does not require the expression of CBF or COR genes in A. thaliana [115] was scored in the up-subset. Moreover, we found also in the up-expressed subset genes coding a VERNALIZATION INDEPENDENCE 4 protein (AT5G61150, homologous to LEO1), involved in vernalization response in many A. thaliana ecotypes [134] and a sensitive to freezing 2 protein (AT3G06510, SFR2, a constitutively expressed b-glucosidase), which is conserved in all land plants [63] and involved in the response to freezing by protecting chloroplast membrane from damages.
Overall, concerning kinases and TFs annotations, we observed that most of the annotated kinases in up-and down-expressed subsets are known to be differentially expressed in cold stress [16,18,19,22]. In particular, we identified a gene coding calcium/calmodulin-regulated receptor-like kinase 1 (AT5G54590) in the up-expressed subset that is involved in freezing tolerance in A. thaliana [113], and 5 casein kinases I and 3 SnRK which were involved in stress response [38,135]. Elsewhere, more TFs were scored in the up-than in the down-expressed subset, with seven families being represented in both subsets, and eight and 31 families being specific to down-or up-expressed subsets, respectively. In particular, we identified specifically into the up-expressed subset a MADS protein (AT4G24540, AGAMOUS-LIKE 24 protein), known to be implicated as a transcription activator mediating floral transition in response to vernalization [7]. Here, we report for the first time the involvement of genes coding ARID, DBP, EIL3, mTERF, RWP-RK, SWI, and TAZ proteins in the chilling response.

Champagne Specific Responses to Cold and Acquisition of Freezing Tolerance
As we have seen previously, the specific response to cold concerns a greater variety of genes than the common response. We observed, whatever the modulation of gene expression levels (increase or decrease), many genes presented different level of expression in Ch and Te at the starting point of the time course study (T0). These observations again highlight the importance of the constitutive differences, even in the specific responses of the two lines. Indeed, the second subset (TLSR (b), Figure 3C (b)) gathered up-regulated genes in Ch and in Te during LT, but which are more expressed in Ch, as a result of a higher expression level at T0 for this line. Similarly, genes within the third subset (TLSR (c), Figure 3C (c)) are down-regulated during LT and are less expressed in Ch than in Te during LT as a result of a lower level of expression in Ch at T0. Hence, genes from clusters TLSR (b) and TLSR (c) could be related to the specific responses of Ch taking part in the freezing tolerance capacity of this line. In a same manner, first and forth subsets ( Figure 3C (a) and (d)) gathered genes that are related to the specific responses of Te. Considering that these responses (qualitatively and/or quantitatively) are ineffective to bring freezing tolerance to Te, they won't be discussed below. We noted that only very few genes presented a modulation of expression between the two sampling times (3 and 16 days) during LT, suggesting that at day three, most of the defenses against cold stress are implemented.
Within the TLSR (b) subset, six genes were assigned to GO terms related to cold including a gene coding a pre-mRNA-processing factor 31 homolog (AT1G60170, PRP31), which is involved in the regulation of expression of cold-responsive genes (CORs) in A. thaliana [120]. The PRP31 protein possesses a HAT domain, which is also found in 4 tetratricopeptide repeat protein (AT3G53560, AT2G37400, AT3G46790 and AT3G23020, TPR and/or PRP1) present in this subset [64]. The five other genes code a RNA-binding family protein (AT1G70200, presenting RRM/RBD/RNP motifs), which is implicated in cold tolerance by 23S ribosomal RNA processing in A. thaliana [64], two chloroplastic extern membrane protein 16-1 (AT2G28900) induced in low temperature in land plants [65], a Fibrillin-1a (AT4G04020) which is involved in response to freezing [66], and a Glycine-rich RNA-binding protein 2 (AT4G13850), which confers freezing tolerance after a cold acclimation period [116]. The KEGG annotation revealed that 28.3% of the genes were mapped on other pathways than those related to "Metabolism" and were rather related to RNA and RNA-related pathways. Thus, the RNA metabolism and more particularly RNA-binding with two Pentatricopeptide repeat-containing (AT1G11290 and AT1G20300, PPR motif), one RNA-binding (AT3G08620, KH domain) family, two RNA-binding (AT1G70200 and AT3G20890, PRM/RBD/RNP motifs), one Helicase (AT3G08620, KH domain), and one Glycine-rich containing domain (AT4G13850, RRM domain) proteins in this subset seem to be an important component of the specific cold response in Ch pea line, as well as in A. thaliana [133,136] and Oryza sativa [137].
The post-translational regulation using kinase proteins is considered as a key feature in plant response against cold stress [138]. Among the 11 up-expressed genes coding kinases in the subset TLSR b, two code calcium-dependent protein kinase 6 (AT4G14580 [57] and AT2G17290 [139]) and one a receptor-like protein kinase FERONIA (AT3G51550 [140]) that are involved in stomatal closure control in relation to cold conditions. Moreover, a gene coding a leucine-rich repeat receptor-like serine/threonine-protein kinase implicated in jasmonic acid and ethylene-dependent systemic resistance (AT3G14840 [141]) is scored. In the earlier studies [28,142], the evidence that jasmonate metabolism could play a role in freezing tolerance was provided. TF constitute 9.21% of the genes of this subset. Excepted SWI/SNF-BAF60b, which is signaled in this study for the first time in cold responses, all other are already known to be up or down regulated under cold treatment in legumes species [16,19,22].
The KEGG annotation within the TLSR (c) subset revealed that most of the genes were associated to "Metabolism". A total of 9 genes were assigned to GO terms related to cold, including genes coding an acyl-CoA-binding domain-containing protein 6 (AT1G31812) that binds phosphatidylcholine in phospholipid metabolism [143], a chlorophyll a/b-binding protein 3-1, chloroplastic (AT1G61520) belonging to the light-harvesting complex in photosystem I [67], two NADPH-dependent aldo-keto reductase, chloroplastic (AT2G37770) that detoxifies a range of toxic aldehydes and ketones produced during stress [68], a bidirectional sugar transporter SWEET17 (AT4G15920) involved in fructose transport [121], and a 3-hydroxyisobutyryl-CoA hydrolase 1 (AT5G65940) that plays a role in peroxisomal metabolism in cold stress signaling and plant tolerance to cold stress, by the degradation of valine [69]. We noticed also the presence of two genes coding a PLAT domain-containing protein 1 (AT4G39730), which functions as positive regulator of abiotic stress tolerance [70], a serine/thereonine protein kinase 2 (AT3G08720), which is a downstream effector of the target of rapamycin signaling pathway (TOR) that presents an increase of protein activity via a phosphorylation induced under cold treatment [71]. All of the TF of TLSR (c) subset were already recorded in cold and freezing tolerance in legumes species, validating once again their participation as key factors of the cold acclimation process in Ch.

Conclusions
In most published studies dealing with cold acclimation, gene expression analyses have been realized at the beginning (1 to 3 h) [12,14] or within the 24 h [19] after the acclimation period. Since we look at gene expression after three and sixteen days of acclimation, we expected to identify novel differentially expressed genes. Hence the expressed genes in this work should be involved in subsequent stages of the cold response and more downstream of metabolic chains.
From a total of 4489 differentially expressed genes, we observed the importance of the constitutive differences in gene expression between the two lines. In particular, we identified more preferentially expressed genes related to RNA metabolism in Te, and to protein metabolism in Ch.
Elsewhere, we observed that most of the genes we identified as involved in the freezing tolerance presented similar modulations (activation/repression) in Ch and Te, but with different levels of expression at the beginning of the time course study. This reveals again the importance of initial differences, in the specific responses of the two lines. According to the examination of clusters gathering genes related to chilling and FT, we were surprised to observe less genes modulated by the cold stress in Ch compared to Te. This probably means that Ch presents constitutive and/or induced mechanisms that are more efficient to get over cold than those implemented in Te.
Furthermore, we have identified many TFs, which are linked here for the first time to cold responses. Overall, genes whose expression is for the first time correlated with cold response could open new horizons in the use of genetic diversity of low temperature responses in pea.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2223-7747/8/8/288/s1, Figure S1: Workflow for statistical analyses, Figure S2: Functional classification by Gene Ontology, Figure S3: Comparison of the qPCR and the RNA-seq data, Figure S4: GO terms enrichment of Biological Processes and Cellular Component, Table S1: Statistical overview of sequencing and transcriptome assembly data, Table S2: KEGG pathways repartition and classification, Table S3: Kinases repartition in each subset of DEGs.