Changes in Lolium perenne transcriptome during cold acclimation in two genotypes adapted to different climatic conditions

Activation of numerous protective mechanisms during cold acclimation is important for the acquisition of freezing tolerance in perennial ryegrass (Lolium perenne L.). To elucidate the molecular mechanisms of cold acclimation in two genotypes (‘Veyo’ and ‘Falster’) of perennial ryegrass from distinct geographical origins, we performed transcriptome profiling during cold acclimation using RNA-Seq. We cold-acclimated plants from both genotypes in controlled conditions for a period of 17 days and isolated Total RNA at various time points for high throughput sequencing using Illumina technology. RNA-seq reads were aligned to genotype specific references to identify transcripts with significant changes in expression during cold acclimation. The genes induced were involved in protective mechanisms such as cell response to abiotic stimulus, signal transduction, redox homeostasis, plasma membrane and cell wall modifications, and carbohydrate metabolism in both genotypes. ‘Falster’ genotype, adapted to cold climates, showed a stronger transcriptional differentiation during cold acclimation, and more differentially expressed transcripts related to stress, signal transduction, response to abiotic stimulus, and metabolic processes compared to ‘Veyo’. ‘Falster’ genotype also showed an induction of more transcripts with sequence homology to fructosyltransferase genes (FTs) and a higher fold induction of fructan in response to low-temperature stress. The circadian rhythm network was perturbed in the ‘Veyo’ genotype adapted to warmer climates. In this study, the differentially expressed genes during cold acclimation, potentially involved in numerous protective mechanisms, were identified in two genotypes of perennial ryegrass from distinct geographical origins. The observation that the circadian rhythm network was perturbed in ‘Veyo’ during cold acclimation may point to a low adaptability of ‘Veyo’ to low temperature stresses. This study also revealed the transcriptional mechanisms underlying carbon allocation towards fructan biosynthesis in perennial ryegrass.


Background
A period of low temperature stress (cold acclimation) can induce protective mechanisms, leading to morphological, physiological, and biochemical changes that are required for the acquisition of freezing tolerance in coldtolerant plants [1]. These protective mechanisms include both alterations in gene expression, and metabolic re-adjustments. Previous studies have demonstrated the induction of genes involved in protective mechanisms such as cell redox homeostasis, signal transduction, cell wall and plasma membrane modifications, and metabolic re-adjustments during cold acclimation in a wide range of plant species [2][3][4]. Some cold regulated genes encode cryoprotectant proteins [5] and protect the cells from dehydration associated with freezing. Water-soluble carbohydrates (WSCs) also play protective roles during abiotic stresses. They can act as cryoprotectants [6], osmoprotectants [7] and also signaling molecules [8]. They might also play roles in neutralization of reactive oxygen species (ROS) [9] and membrane stabilization [10] during abiotic stresses. The likely scenario is that cold acclimation is a complex process with multiple components which are controlled by multiple regulatory mechanisms [11,12].
In temperate grasses, a large proportion of the genome is cold-responsive [27,28] as they are highly adaptive to the cold conditions. Induction of genes encode proteins such as cold-regulated, dehydration-responsive, and ice recrystallization inhibition proteins involved in protective mechanisms, has been shown in perennial ryegrass in response to low-temperature stress [28,29]. However, stress response characteristics vary between different genotypes, especially between genotypes with very different geographic origins. Comparisons of transcriptomic data between such genotypes provide information about the plant adaptations to cold environments. We have recently demonstrated the improved cold stress tolerance and changes in fructan composition in 'Veyo' and 'Falster' genotypes of perennial ryegrass during cold acclimation [15]. 'Falster' is a Danish ecotype that is well adapted to cold climates, and 'Veyo' a Mediterranean variety well adapted to warmer climates [30]. 'Falster' showed a better adaptation during cold acclimation and faster recovery after freezing compared to 'Veyo' [15]. Furthermore, the genotypes differ in that 'Falster' must undergo a period of low temperature (vernalisation) in order to flower. 'Veyo' does not require a period of vernalisation to flower. A recent study has shown that both 'Veyo' and 'Falster' respond differently on a transcriptional level during vernalisation [31]. It would therefore be expected that 'Veyo' and 'Falster' would also respond differently on a transcriptional level during cold acclimation. Here we used these two types of perennial ryegrass to study the transcriptome profiles during cold acclimation using high throughput sequencing technologies and thereby gain a deeper insight into molecular mechanisms of cold acclimation. The specific aims of this study were: i) to identify candidate genes differentially expressed in perennial ryegrass during cold acclimation, ii) to identify molecular pathways differentiated between genotypes adapted to cold and warm climates, and iii) to identify the transcriptional mechanisms underlying carbon allocation towards fructan biosynthesis during cold acclimation.

Differential expression of genes during cold acclimation
High-throughput RNA sequencing, generated~2 Gb reads per sample. A total number of 157,264,629 reads of 50 bp were generated for the genotype 'Veyo' and a total of 151,608,297 reads were generated for 'Falster'. In the 'Veyo' Trinity assembly, 50% of the total assembly was present in contigs of at least 1712 bp. In the 'Falster' Trinity assembly, 50 % of the total assembly was present in contigs of at least 1671 bp. The longest assembled contigs in 'Veyo' and 'Falster' had 15,228 bp and 15,362 bp, respectively. The average contig lengths in 'Veyo' and 'Falster' were 1078.60 bp and 1052.38 bp, respectively.
In total 1, 45,805 transcripts in 'Veyo' and 1, 44,062 transcripts in 'Falster' were identified. In a series of pairwise comparisons between the days, we identified 1874 differentially expressed transcripts in 'Veyo' and 2567 in 'Falster'. There were 263 differentially expressed transcripts common for both 'Veyo' and 'Falster' (Fig. 1). Hierarchical cluster analysis performed using transcripts that were significantly differentially expressed (P ≤ 0.05) between any pairwise comparison (Additional file 1).
A total of 24 clusters of expression profiles were identified using K-means algorithm with distinguishable expression profiles during cold acclimation in each 'Veyo' or 'Falster' (Figs. 2 and 3, Additional files 2 and 3). The sudden drop of temperature from 20°C to 7°C lead to the rapid up-or down-regulation of transcripts in both 'Veyo' (Fig. 2, clusters C 2, C 6, C 10, and C 14) and 'Falster' (Fig. 3, clusters C 7, C 8, C 9, C 12, C 13, C 15, and C 24), indicative of transcripts involved in abiotic stress responses. The best BLAST hit descriptions of differentially expressed transcripts in 'Veyo' and 'Falster' are shown in Additional files 4 and 5, respectively. In general, there was an excellent correlation between the RNA-Seq results and the quantitative RT-PCR results (Additional file 6).
GO enrichment analysis using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool [32] identified major biological, metabolic, and cellular processes significantly enriched (P ≤ 0.05) in 'Veyo' and 'Falster' during cold acclimation (Additional files 7, 8, 9, 10,11 and 12). Terms that were enriched in both 'Veyo' and 'Falster' were mainly related to the processes such as cell redox homeostasis that is important to maintain the redox environment of cells under stress conditions, signal transduction, response to abiotic stimulus, and carbohydrate metabolism. Analysis of terms under "cellular compartment" showed the enrichment of proteins associated with plasma membrane in "Veyo" and "Falster" (Additional files 11 and 12). These proteins included calcineurin B-like protein and elicitor-responsive protein, plasma membrane intrinsic protein aquaporins, cellulose synthase, ras-related protein RIC2, and secretory carrierassociated membrane protein. Terms such as extracellular region, chloroplast, golgi apparatus, and cell-wall were also found to be significantly enriched in both 'Veyo' and 'Falster'. The DAVID analysis also showed the protein domains that were significantly enriched (P<0.05) in the differentially expressed transcripts (Additional files 13 and 14). Enriched SMART domains revealed the presence of functional domains such as serine-threonine protein kinase catalytic (S_TKc), serine-threonine phosphatases, family 2C, catalytic (PP2Cc), heat shock factor (HSF), cyclin, and formin homology (FH) domains in both genotypes.
Divergence of gene expression-profiles between "Veyo" and "Falster" genotypes during cold acclimation Our results show more differentially expressed transcripts in 'Falster' compared to 'Veyo' during cold acclimation. Gene Ontology annotation (GO) of the possible functions of unigenes using Blast2GO showed transcripts related to stress, signal transduction, response to abiotic stimulus, and metabolic processes such as lipid, amino acid, and carbohydrates ( Fig. 4). 'Falster' showed more transcripts in each of these categories compared to 'Veyo'. Analysis of protein domains using DAVID tool revealed the enrichment of EF-hand, calcium binding domain only in 'Falster' (Additional files 14).
In order to make comparisons between the expression profiles of 'Veyo' and 'Falster' , we used common reference L. perenne transcriptome [33] for differential analysis. Principle component analysis (PCA) showed that the transcription profiles were separated according to genotype on PC1, which explained 66 % of the variance (Fig. 5). The second PC separated samples according to treatment and explained 12 % of the variance. There was greater separation on PC2 for samples from the 'Falster' genotype, indicating a stronger transcriptional differentiation during cold acclimation.
There were two KEGG pathways that were significantly enriched in the 'Veyo' differentially expressed gene set; these were 'Phenylpropanoid biosynthesis' (P = 2.80E-02) and 'Circadian rhythm' (P = 2.90E-02). Down regulation of genes encode enzymes involved in phenylpropanoid biosynthesis such as phenylalanine ammonia-lyase, 4coumarate-CoA ligase, and O-methyl transferase, and also the down regulation of genes encode two-component response regulators involved in circadian rhythm [34], in response to the sudden drop of temperature from 20°C to 7°C, was observed in 'Veyo' (Fig. 2, cluster C 10). Only a single KEGG pathway was significantly enriched in the 'Falster' differentially expressed gene set; this was 'Proteasome' (P = 2.00E-03). The genes encode proteasome regulatory subunit proteins up regulated in 'Falster' during cold acclimation (Fig. 3, cluster C 8). Analysis of terms under "biological processes" showed that terms such as rhythmic process (P = 6.30E-03), two-component signal transduction system (P = 2.60E-02), nuclear division (P = 3.30E-02), and mitosis (P = 3.30E-02) were only enriched in 'Veyo' (Additional file 7).

Induction of genes involved in fructan metabolism and carbon allocation towards fructan biosynthesis during cold acclimation
Among the differentially expressed transcripts, 345 transcripts were related to carbohydrate metabolism or carbohydrate binding (Additional files 4 and 5). Transcripts potentially involved in fructan biosynthesis were classified into four groups, based on their deduced amino acid sequence similarity to the previously characterized genes. These four groups are FTs, FEHs, invertases, and other genes potentially involved in carbon flux diversion towards fructan biosynthesis (Additional file 15).
Seven transcripts from 'Falster' that were up regulated (Fig. 3, clusters C 8 and C 13), showed deduced amino acid sequence similarity with the functionally characterized Lp1-SST [20], Lp6-SFT, and Lp6G-FFT [22], of which five showed over 90 % deduced amino acid sequence similarity over an alignment of 500 amino acids (Additional file 15). Three transcripts from 'Veyo' (Fig. 2, cluster C 14), and four transcripts from 'Falster' (Fig. 3, clusters C 3 and C 7) that were down regulated, showed deduced amino acid sequence similarity with FT-like protein (LpFTL) [35]. Seven transcripts from 'Veyo' were classified as invertases, of which two were CWIs. One transcript from 'Veyo' (Fig. 2, cluster C 14) and two transcripts from 'Falster' (Fig. 3, cluster C 7) classified as soluble acid invertases, were down regulated during cold acclimation. Two transcripts from 'Falster' which were down regulated (Fig. 3, cluster C 7), showed over 88 % deduced amino acid sequence similarity with 6-FEH from Timothy (Phleum pratense) [36]. Another transcript from 'Falster' that was up regulated (Fig. 3, cluster C 8) showed 99 % sequence similarity with Lp6-FEH. Fig. 4 The pie charts of the biological processes of unigenes differently expressed during cold acclimation. The pie charts of the biological processes (GO level 4) of unigenes differently expressed in "Veyo" (a) and "Falster" (b) genotypes during cold acclimation. The colours depict distinct categories. Number of transcripts associated to each category is shown in brackets A model for carbon flux diversion towards fructan biosynthesis in perennial ryegrass is shown in Fig. 6. The differentially expressed transcripts during cold acclimation potentially involved in carbon flux diversion towards fructan biosynthesis showed deduced amino acid sequence similarity to proteins such as glyceraldehyde-3-phosphate dehydrogenase (GAPDH) ( Fig. 7. The fructan content was increased around two-fold in 'Veyo' and around five-fold in 'Falster' during cold acclimation. Sucrose content was increased around five-fold in both genotypes.

Activation of protective mechanisms during cold acclimation
Transcriptome analysis showed the activation of numerous protective mechanisms in 'Veyo' and 'Falster' during cold acclimation. Moreover, a number of processes enriched in the sets of genes differentially expressed during cold acclimation were common to both 'Veyo' and 'Falster'. This is not surprising because both genotypes have shown acquisition of freezing tolerance during cold acclimation [15]. Consistent with the previous studies [2,4,37] our results also showed the induction genes involved in protective mechanisms such as redox homeostasis mechanisms, signal transduction, cell membrane stabilization, and carbohydrate metabolism during cold acclimation. Induction of cold-responsive genes such as cold-regulated, dehydrationresponsive, late-embryogenesis-abundant, and ice recrystallization inhibition homolog genes has been shown previously in perennial ryegrass cv. 'Caddyshack' during cold acclimation [28]. Our results showed the induction of these genes also in 'Veyo' and 'Falster' genotypes during cold acclimation. Analysis of differentially expressed transcripts during cold acclimation revealed that up-or downregulation of large number of transcripts in both 'Veyo' and 'Falster'. Down regulation of genes involved in mechanisms such as photosynthesis and respiration during cold acclimation, has been previously shown [28,38]. Our results also support the hypothesis that in addition to gene induction, gene repression might also play roles in cold acclimation [4].
Induction of genes encode plasma membrane associated proteins such as calcineurin B-like protein, and elicitor-responsive protein might play important roles in Ca 2+ -sensing and defense signaling [39,40] in both genotypes. The secretory carrier-associated membrane protein (SCAMP6) protein enriched in "Veyo" and "Falster" might be involved in the membrane trafficking of polysaccharides. The co-localization of NtSCAMP2 from Nicotiana tabacum and pectin in vesicular compartments has been shown [41], suggesting that NtSCAMP2 plays a role in intracellular vesicular-mediated transportation of polysaccharides. Induction of genes encode plasma membrane intrinsic protein aquaporins has been observed in both genotypes. These proteins might be involved in cellular hydraulic conductivity, transportation of solutes across membranes, and carbohydrate compartmentation [42,43] during cold acclimation. Enrichment of functional domains such as S_TKc, PP2Cc, HSF, cyclin, and FH during cold acclimation indicates the activation of processes such as phosphorylation and dephosphorylation of target enzymes, protein folding and protein translocation, cytoprotection, transcriptional regulation, and cytoskeleton rearrangements [44][45][46][47].
'Falster' genotype, adapted to cold climates showed a stronger transcriptional differentiation during cold acclimation suggestive of differences in the cold acclimation responses between these two perennial ryegrass genotypes adapted to different climatic conditions. Previous studies have demonstrated coordinated expression of functionally diverse FTs and the association of specific isoforms with the high sugar trait in perennial ryegrass  Expression profiles are only shown in cases where significant differential expression was identified during cold acclimation. The x-axis represents time during cold acclimation (d) and y-axis represents relative gene expression (arbitrary units) [20,48]. Our results show the induction of more transcripts with sequence homology to fructosyltransferase genes (FTs) in the 'Falster' genotype adapted to cold climates, compared to 'Veyo' in response to low-temperature stress. Some of them might represent different isoforms of FTs. Induction of specific and / or multiple isoforms of FTs might play roles in higher fold induction of fructan in cold adapted genotypes of perennial ryegrass.
The KEGG pathways 'Phenylpropanoid biosynthesis' and 'Circadian rhythm' were significantly enriched in the differentially expressed gene set of 'Veyo' originating from Italy. Phenylpropanoid metabolism produces a large array of secondary metabolites that are involved in protective mechanisms against abiotic stresses [49,50]. The observation that the down regulation of genes involved in these pathways may point to a low adaptability of this Italian variety to low temperature stresses. Our results indicate the up regulation of proteasome pathway in 'Falster' during cold acclimation. This may point to a high adaptability of 'Falster' to low temperature stresses. The proteasome pathway degrades misfolded proteins that are generated in response to environmental stresses [51].
Apart from the induction of genes involved in numerous protective mechanisms, post-transcriptional and post-translational mechanisms regulating the function of proteins [52] also play roles in abiotic stress tolerance. Our results showed divergence of gene expressionprofiles between "Veyo" and "Falster" genotypes during cold acclimation. The proteomic response to cold acclimation can be also varied between perennial ryegrass genotypes [53]. Further investigation of such differences between genotypes is important for a better understanding of biological processes involved in cold acclimation.

Coordinated gene expression directs carbon flux towards fructan biosynthesis
It has already been suggested that fructans play important roles in abiotic stress tolerance in temperate grasses [13,15]. Therefore, in this study we looked at the expression of genes involved in fructan metabolism. Genes encoding transcription factors such as MYB, bZIP, AP2/ EREBP, WRKY, and NAC induced in 'Veyo' and 'Falster' during cold acclimation might play important roles in coordinating genes involved in protective mechanisms [54][55][56]. Previous studies have shown the up-regulation of 1-SST and 6-SFT promoter-driven reporter genes by a R2R3 MYB transcription factor from wheat [57], suggesting that transcription factors might play important roles in coordinating fructan-related gene expression. Consistent with the previous studies [28,31], our results , sucrose (c), and fructan (d) contents before cold acclimation when the plants at 20°C, and after 17 days at 7°C are shown. Data represent the mean ± SE, obtained from three replicates of the analysis. Asterisks indicate significant differences (P < 0.05) compared to before cold acclimation also show the induction of genes encoding transcription factors potentially involved in coordination of fructanrelated gene expression in perennial ryegrass during cold acclimation. A large number of genes encoding kinases and phosphatases were also differentially expressed during cold acclimation in both 'Veyo' and 'Falster' (Additional files 4 and 5). Some kinases and phosphatases are involved in the sucrose-mediated induction of fructan synthesis in plants [25,26].

Conclusions
RNA-Seq enabled an untargeted profiling of the transcriptome of two perennial ryegrass genotypes adapted to different climatic conditions. The stronger transcriptional differentiation during cold acclimation in genotype 'Falster' explains its higher adaptability to cold climates compared to 'Veyo'. This study also describes the differentially expressed genes involved in fructan metabolism and other protective mechanisms during cold acclimation that can be useful for further research such as the functional analysis of genes, and comparative genomic studies.

Plant materials and growth conditions
Perennial ryegrass (Lolium perenne L.) variety 'Veyo' originating from Italy and adapted to warmer climates, and ecotype 'Falster' originating from Denmark and adapted to cold climates [30] were propagated vegetatively in a glasshouse. Following propagation, the plants were grown in a glasshouse (photoperiod of 9 h:15 h, light:dark) at~20°C for 12 weeks. The plants were cut (12 cm above the soil level) two times during the 12 weeks in glasshouse, and another time prior to transfer into a controlled environment chamber. The plants (at the preelongation stage) were then transferred to a growth chamber for 7 d (photoperiod of 9 h:15 h, light:dark), with 450 μmol photons m −2 s −1 light intensity at 20°C. The plants were cold-acclimated in a controlled growth room for 17 d (photoperiod of 9 h:15 h, light:dark), with 450 μmol photons m −2 s −1 light intensity and relative humidity of~70 % at 7°C as previously described [15].

Measurement of total fructan
The blades of fully developed leaves harvested in three biological replicates before and after cold acclimation were freeze-dried and ground to powder. Soluble carbohydrates were extracted from 100 mg of sample with 25 ml of 0.1 M acetate buffer (pH 5.0) for 1 h at 65°C. Extracts from each sample (2 ml) were hydrolysed using an equal volume of 0.074 M H 2 SO 4 for 70 min at 80°C. The amounts of glucose and fructose before and after acid hydrolysis were measured, using the hexokinase-phosphoglucose isomerase glucose-6-phosphate dehydrogenase system, by calculating the sucrose and fructan contents, as previously described [13].

RNA-Seq analysis
We have previously demonstrated an improved cold stress tolerance in both 'Veyo' and 'Falster' during the first 17 d of cold acclimation [15]. Therefore, in this study, same time points were selected for RNA-Seq analysis to elucidate the molecular mechanisms of cold acclimation. The blades of fully developed leaves were harvested in three biological replicates at the onset of the daily photoperiod on d 0, 9, 13, and 17 of cold acclimation at 7°C. The samples were frozen in liquid nitrogen immediately after harvesting and stored at -80°C. until the analysis of carbohydrates and gene expression levels was performed. The samples were ground separately in liquid nitrogen. Total RNA was extracted using RNeasy® Plant Mini Kit (Qiagen), and On-column DNAse I digestion (Qiagen) was performed to avoid genomic DNA contamination, according to the manufacturer's instructions. Total RNA quality was verified using an Agilent 2100 Bioanalyzer (RNA 6000 Nano Assay) and quantified by Quant-iT TM RiboGreen® RNA Reagent assay, according to the manufacturer's protocol. Total RNA samples were prepared for sequencing, using the TruSeq RNA Sample Preparation Kit (Illumina).
Illumina HiSeq 2000 platform was used for highthroughput sequencing (Beijing Genomics Institute, Hong Kong). Where the average quality score in a sliding window analysis (10 % of read length) fell below 20, poorquality bases towards the 3′ end of reads were trimmed using Sickle [61]. Reads shorter than 40 bp after trimming were discarded. De novo assembly of RNA-Seq data was carried out using the Trinity assembler, as previously described [62]. De novo transcriptomes were assembled for 'Veyo' , and for 'Falster' separately (minimum contig length of 200, and minimum kmer coverage of 2).
RNA-Seq Reads were aligned back to their respective assemblies using Bowtie as previously described [63]. The abundance estimates of transcripts separately for each time point (d 0, 9,13, and 17 of the cold acclimation) and genotype ('Veyo' and 'Falster') were calculated using RSEM [64]. Pairwise comparisons were carried out between all the selected time points and differentially expressed transcripts (P ≤ 0.05) were identified using edgeR as previously described [65].
RNA-Seq count data from each time point was normalized using trimmed mean of M-values (TMM) normalization [66] in edgeR to compute fragments per feature kilobase per million read mapped (FPKM) values. These variance stabilized data was used for clustering the transcripts. The transcripts showing differential expression at any time point during cold acclimation were clustered using K-means algorithm. Hierarchical cluster analysis was carried out separately on two sets of differentially expressed transcripts, a differentially expressed set from 'Veyo' , and a differentially expressed set from 'Falster'. Each set contained transcripts that were significantly differentially expressed (P ≤ 0.05) between any paired comparison for samples from 'Veyo' or 'Falster'. The differentially expressed gene sets were mapped onto Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway maps to identify networks differentially regulated between 'Veyo' and 'Falster'.
The differentially expressed transcripts were annotated using Blast2GO as previously described [31]. Functional annotation information was also assigned to each transcript by using each sequence as a query in a BLASTx against a UniProt database. Gene ontology annotation (GO) of the possible functions of unigenes was carried out using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) bioinformatic tool [32]. In order to compare the differentially expressed transcript sets identified in 'Veyo' and 'Falster' , BLASTn (E-value threshold 10 e-10 ) was used to map the transcripts back to a common reference L. perenne transcriptome [33] as previously described [31]. Principle component analysis (PCA) plot was generated based on variance stabilized expression data in order to separate transcription profiles according to genotype and treatment.
The RNA-Seq results were validated using quantitative reverse transcriptase-polymerase chain reaction (qRT-PCR) standard curve method for absolute quantification on a selection of genes potentially involved in fructan metabolism that were identified as differentially expressed during cold acclimation. PCR procedure and the primers used for the expression analysis of Lp1-SST, Lp1-FFT, Lp6G-FFT, and Lp1-FEH genes have been previously described [15]. PCR primers and the quantitative RT-PCR primers used for the expression analysis of LpFTL, LpCWI-1 and LpCWI-2 genes are shown (Additional files 16 and 17).

Availability of supporting data
Illumina sequence data are available from ArrayExpress with the accession number E-MTAB-2779.