Cold stress, such as chilling and freezing temperatures, is a major environmental factor that influences plant distribution, growth and development [1]. A few plants cope with cold stress by acquiring cold tolerance, a process termed cold acclimation. During this process, various biochemical and physiological changes occur and make the plant more tolerant under chilling or freezing conditions [2]. Cold acclimation involves the signal transduction of various genes responding to cold stress, transcriptional regulation and posttranslational regulation of transcription factors [1]. The induction of cold responsive (COR) genes have been shown to result in cold acclimation for alfalfa [3]. Many transcription factors involved in cold tolerance have been identified [4]. C-repeat binding factors (CBFs), also known as dehydration-responsive element-binding (DREB) proteins, are a group of transcription factors that are well-studied in plant cold stress. CBFs can bind to the promoters of COR genes and activate their expression [1]. The ectopic expression of CBF1 activates the expression of COR and induces cold tolerance in Arabidopsis [5]. Regulators upstream of CBF have also been identified, including ICE1, MYB15, CAMTA and HOS1 [1,6-8]. Although considerable studies have documented the molecular mechanism of plant cold tolerance, deep studies are necessary to highlight the transcriptional regulation of cold tolerance.

Kentucky bluegrass (Poa pratensis ) is a perennial, cool-season species of grass native to Europe and is widely used for making lawns in parks, gardens, golf courses and sports fields [9]. This plant is strongly resistant to cold stress and can survive at temperatures as low as -38°C. This adaptability suggests that Kentucky bluegrass has a distinct molecular mechanism that enables it to adapt to chilling or freezing stress.

In this study, we sequenced and analyzed the transcriptome of Kentucky bluegrass under normal and cold treatment conditions using RNA-seq, aiming to identify more transcription factors related to cold tolerance of Kentucky bluegrass. Differentially expressed genes (DEGs) under cold treatment conditions were analyzed. Genes belonging to several transcription factor families were identified as DEGs. The expression levels of a few selected genes were confirmed by quantitative RT-PCR. Our results provide useful information for understanding the Kentucky bluegrass response to cold stress and offer new information that will be useful for further studies.

Plant materials and treatments

Kentucky bluegrass (Poa pratensis cv. Merit) plants were grown in 30 cm diameter pods for three months under normal management. At least 30 plants were contained in each pod. When plants reached a height of 10 cm, they were subjected to cold treatment (-5°C) for 24 h in a growth chamber. During the cold treatment, plant leaves were sampled every 6 h for a total of 5 samples (0, 6, 12, 18 and 24 h) for each treatment. The samples without treatment (control, Ctrl) and those collected at 6 h after treatment (Cold 6 h) were used for RNAseq. All samples were immediately frozen in liquid N2 and stored at -80°C for further use. Relative electrical conductivity and soluble sugar content measurements were taken according to Zhang et al. [10]. Each plant pod was used as a replicate; three replicates were used for each treatment. For the RNA-seq, samples of these three replicates were mixed together for RNA extraction. In total, two samples (Cold 6 h and Ctrl) were used for RNA-seq.

RNA extraction and RNA-seq

Total RNA was extracted using the cetyl trimethyl ammonium bromide (CTAB) method [11]. Briefly, mRNA was isolated from total RNA using oligo(dT) beads and was cut into short fragments with a 10 × fragmentation buffer. First-strand cDNAs were synthesized with random primers using the short mRNA fragments as templates. Second-strand cDNAs were synthesized with DNA polymerase I (TaKaRa) and were purified using a QIAquick PCR Purification Kit (Qiagen). The purified cDNAs were subjected to end reparation and polyadenylation and were mixed with Solexa adapters. Suitable fragments (200-700 bp) were recovered from an agarose gel (Agencourt AMPure XP, Bechman) and amplified by PCR to generate cDNA library. The library was sequenced using an Illumina HiSeqTM 2500. The length of each reads was 125 bp.

Bioinformatics analysis

Low-complexity reads and reads containing adapter sequences were removed from the raw reads data to generate clean reads. A de novo assembly of the transcriptome was constructed from the resulting clean reads using the short-reads assembly program Trinity [12]. Reads with overlaps were assembled to form contigs, which were then assembled to unigenes. Unigenes were first aligned to the NCBI non-redundant (NR) protein database using BLASTx (E-value < 10-5). All assembled unigenes were also annotated using public protein databases, including NR [13], Swiss-Prot [14], Gene Ontology (GO) [15] and Clusters of Orthologous Groups (COG) [16].

Differentially expressed genes (DEGs) were analyzed using the method described by Li and Dewey [17], and the false discovery rate (FDR) was used to determine P-value thresholds via multiple testing [18]. The fragments per kilobase of transcript per million mapped reads (FPKM) method was used to calculate the rate of DEGs [19], with FDR ≤ 0.01 and a fold change value ≥ 2. The heat map for DEGs involved in cold stress was constructed using Cluster 3.0. All the raw data has been deposited into NCBI Sequence Read Archive (SRA) under accession number SRP075404.

qRT–PCR

First-strand cDNA was synthesized from 500 ng of total RNA using the M-MLV RTase cDNA Synthesis Kit (Cat#D6130, TaKaRa). The cDNA solution was diluted 10 times with H₂O and was subsequently used as templates for qRT-PCR assays. qRT-PCR was conducted as described by Tan et al. [20]. Specific primers for each gene were designed using Primer3 (http://frodo.wi.mit.edu/) and are listed in Supplementary Table S1. The actin gene, homologue of Brachypodium distachyon [21], which showed no differential expression in our study was used as an internal control. Three technical replications were conducted for each measurement.

Physiological characteristics of Kentucky bluegrass and RNA-seq

Soluble sugar content and electrical conductivity (EC) are regarded as important protection factors in plants under environmental stress [22]. In this study, we first measured the soluble sugar content and EC of Kentucky bluegrass plants under cold treatment (-5°C). The soluble sugar content started to increase after 6 h of cold treatment and continued to increase over the remaining treatment duration. Furthermore, the soluble sugar content was significantly increased after Cold 6 h (Figure 1). Moreover, the EC measurements showed the same trend as those of the soluble sugar contents (Figure 1). These results suggested that cold treatment resulted in physiological changes in the Kentucky bluegrass plants.

Figure 1: Soluble sugar content and relative electrical conductivity of Kentucky bluegrass during cold treatment. Kentucky bluegrass plants were grown in a growth chamber at -5°C for 24 h. The leaves were sampled every 6 h.

Control sample (Ctrl) and sample collected at 6 h after cold treatment (Cold 6 h) were used for RNA-seq. In total, we obtained 19 and 22 million reads from the Cold 6h and Ctrl samples, respectively (Table 1). All of the clean reads were assembled into 9,006,576 contigs, with an average contig length of 53.8 bp. The contigs were joined to generate 75,934 unigenes with an average length of 793.6 bp (Table 2). More than 63% of reads were mapped to unigenes (Table 1). To annotate the assembled unigenes, sequence similarity searches (Evalue < 10-5) were performed against public databases, including the NCBI NR, Swiss-Prot protein, GO, and COG databases. In total, 53,762 unigenes (70.8%) were successfully annotated in public databases (Table 3).

Table 1: Number of sequenced and mapped reads.

Table 2: Statistics of sequence assembly.

Table 3: Number of functional annotations for all unigenes in public databases.

Identification of DEGs between the Ctrl and Cold 6 h samples

To identify genes involved in cold tolerance, we compared the transcripts of control and Cold 6 h samples by differential expression analysis using EBSeq software [23]. The transcript abundance of each gene was estimated using FPKM. In total, 3,896 unigenes were differentially expressed (log2Ratio ≥ 2, FDR ≤ 0.01), out of which 2,410 unigenes were down-regulated and 1,486 unigenes were up-regulated in the Cold-6hrs samples (Figure 2). Out of these 3,896 unigenes, 3,207 were successfully annotated in public databases (Supplementary Table S2). GO analysis was employed to analyze the DEGs. Of the 3,207 DEGs identified, 1,944 were successfully categorized into GO groups. The largest numbers of DEGs were found with the phrases ‘cell part’ or ‘cell’ in the cellular component category, ‘binding’ or ‘catalytic activity’ in the molecular function category, and ‘cellular process’ or ‘metabolic processes’ in the biological process category (Figure 3).

Figure 2: Differentially expressed genes (DEGs) in cold treatment (Cold 6 h) and control samples. Significantly up- or down-regulated genes are marked in blue, and genes showing no significant differential expression are marked in red.

Figure 3: Gene Ontology (GO) assignment of all DEGs. The unigenes were mapped to three main categories: cellular component, molecular function and biological process. The right-hand y-axis indicates the number of annotated unigenes. The color of the bars corresponds to the color of numbers on right-hand y-axis.

COG analysis of DEGs

The DEGs were also subjected to COG analysis. Of the 3,207 DEGs, 979 were annotated in the COG database (Supplementary Figure S1). The largest numbers of DEGs were found in the ‘General function prediction only’ COG category, followed by the ‘Replication, recombination and repair’ and ‘Transcription’ categories.

Identification of DEGs involved in cold stress

We then analyzed the genes that were identified as closely related to cold stress in previous reports which also had differential expression levels in the Cold 6 h and Ctrl samples in this study. These genes included transcription factors, cold-induced proteins, antioxidant enzymes, osmoregulation proteins, and proteins related to Ca²⁺ and abscisic acid (ABA) signaling (Supplementary Table S3). We generated a heat map according to their transcription levels (Figure 4). Cold shock proteins (CSPs) are involved in the cold stress response [24]. Five CSP genes were up- or down-regulated by cold treatment in our data. Two heat shock proteins (HSPs) were also down-regulated by cold treatment [25]. The transcriptional regulation of cold tolerance has been widely reported. DREB and the transcription factor MYB were significantly affected by cold treatment. These transcription factors have been described as major factors involved in plant coldstress responses [6-7,26-27]. Additionally, some antioxidant enzymes and osmoregulation proteins were up-regulated by cold treatment in our data. These proteins included peroxidase (POD), water stress-induced protein and a sugar transporter. Under cold stress, reactive oxygen species (ROS) accumulate in plants, and osmotic pressure becomes imbalanced in cells [28]. The accumulation of antioxidant enzymes and osmoregulatory proteins may protect cold-stressed cell membranes.

Figure 4: Heat map of relative expression levels for DEGs related to cold stress. DEGs including transcription factors, cold-induced proteins, antioxidant enzymes, osmoregulation proteins, and proteins related to Ca2+ and abscisic acid (ABA) signaling were selected for constructing the heat map.

.Identification of transcription factors involved in cold stress

The transcriptional regulation of cold stress has been widely documented in model plants [1]. One of the most important transcription factors is CBF (also known as DREB), which belongs to the Apetala2/ethylene response factor family. CBF can bind to the promoters of COR (cold responsive) genes and activate their expression [1,4]. Zhuang et al. [29] has reported that PpCBF3 is induced by CT in Kentucky bluegrass and overexpression of PpCBF3 improves freezing tolerance of transgenic Arabidopsis. In our data, 27 ERFs or DREBs genes were identified, including those tested by qRTPCR, which were up- or down-regulated by cold treatment (Table 4). These results suggested that members of this transcription factor family act as transcription activators and also possibly as transcription repressors in response of cold stress. Because transcriptional regulation is pivotal in the cold acclimation of plants, more work should be conducted to map the regulatory networks among these transcription factors, including transcriptional regulation and protein interactions between transcription factors.

Table 4: DEGs related to transcription factors.

Heat stress transcription factors (HSFs) also showed differential expression after cold treatment (Table 4). HSP is involved in cold stress in plants [30], thus HSFs, as transcription factors also probably participate in this process. The association of NAC proteins and WRKY domain-containing transcription factors with cold stress has been previously reported in various plants [31-34]. We also identified differentially expressed NACs and WRKYs in our data (Table 4). These findings facilitate potential studies focusing on the interactions of different transcription factors in the regulation of cold acclimation. Two ARFs (auxin response factors) were also down-regulated by cold treatment in our data (Table 4); this down-regulation might be caused by the inhibition of hormone signaling by cold treatment. However, whether ARF is directly involved in cold stress in plants remains unclear. More work is required to elucidate the role of ARF in cold stress.

CBF genes have been identified as coding for important transcription factors in cold acclimation in plants. CBF genes regulate more than 170 COR genes in Arabidopsis [35]. Additionally, non-CBF genes, including HSFC1, ZAT12, ZF, ZAT10 and CZF1, also transcriptionally regulate COR expression [36]. These CBF and non- CBF transcription factors only regulate approximately one-third of the 2,000 COR genes [35]; the factors that regulate the other two-thirds of COR genes have yet to be identified. Therefore, the identification of other transcription factors that regulate CORs is essential for a deeper understanding of transcriptional regulation in response to cold stress. In addition to the CBFs, the identified transcription factors in our study (Table 4) may contribute to the regulation of these CORs.

qRT-PCR Validation of RNA-seq Results

To confirm the accuracy of the RNA-seq results, 10 DEGs including ERF and DREB transcription factors were selected for a qRT-PCRbased comparison of their expression levels between the Cold 6 h and Ctrl samples (Figure 5). Three ERF genes were tested. The expression of c116404 was promoted by cold treatment. The expression levels of c117307 and c47216 peaked at 12 h after treatment and subsequently decreased. Two DREB transcription factors showed opposite expression patterns: c115634 was induced, whereas c122138 was suppressed by cold treatment, suggesting their positive or negative regulatory role during cold treatment. The involvement of SWEET (sugars will eventually be exported transporter) protein in cold stress has been reported in various studies [37-39]. AtSWEET15 is upregulated by cold treatment in Arabidopsis [37]. In our data, one SWEET gene (c121279) was highly induced at 12 and 18 h after cold treatment, but another SWEET gene (c115130) was decreased during cold treatment, suggesting that these SWEET genes may act as different roles during cold treatment. HSPs respond not only to high temperature stress but also to cold stress [30]. Two HSPs (c112042 and c117918) had enhanced expression rates that subsequently decreased in our data (Figure 5). The chaperone functions of HSPs most likely maintain cellular homeostasis in plant cells responding to stress. Moreover, a CSP gene (c117281) was promoted at 6 h after cold treatment and this expression pattern is similar to the report by Kim et al. [24] in which an Arabidopsis AtCSP3 is up-regulated by cold treatment. Although we did not perform biological replicates for these two samples in RNA-seq, the expression of cold-related genes was consistent with the results of previous studies; moreover, the verification of DEGs by qRT-PCR in the following studies confirmed the accuracy of the RNA-seq results

Figure 5: Confirmation of gene expression by qRT-PCR. The numbers under x-axis indicate the time under cold treatment (Cold 6 h, -5°C). Error bars represent SD calculated from three biological replications.

We identified 3,896 unigenes that showed differential expression between control and cold-treated Kentucky bluegrass. Several transcription factor families were identified as DEGs. The data presented could be used as a valuable resource for studies on the transcriptional regulation of cold tolerance.

This work was supported by the Innovative Program for Fruit Tree Research of Liaoning Province, China (LNGSCYTX-13/14-4).

All the authors agree on the contents of the paper and declare that they have no conflict of interest.