A transcriptomic analysis of bermudagrass (Cynodon dactylon) provides novel insights into the basis of low temperature tolerance

Cold stress is regarded as a key factor limiting widespread use for bermudagrass (Cynodon dactylon). Therefore, to improve cold tolerance for bermudagrass, it is urgent to understand molecular mechanisms of bermudagrass response to cold stress. However, our knowledge about the molecular responses of this species to cold stress is largely unknown. The objective of this study was to characterize the transcriptomic response to low temperature in bermudagrass by using RNA-Seq platform. Ten cDNA libraries were generated from RNA samples of leaves from five different treatments in the cold-resistant (R) and the cold-sensitive (S) genotypes, including 4 °C cold acclimation (CA) for 24 h and 48 h, freezing (−5 °C) treatments for 4 h with or without prior CA, and controls. When subjected to cold acclimation, global gene expressions were initiated more quickly in the R genotype than those in the S genotype. The R genotype activated gene expression more effectively in response to freezing temperature after 48 h CA than the S genotype. The differentially expressed genes were identified as low temperature sensing and signaling-related genes, functional proteins and transcription factors, many of which were specifically or predominantly expressed in the R genotype under cold treatments, implying that these genes play important roles in the enhanced cold hardiness of bermudagrass. KEGG pathway enrichment analysis for DEGs revealed that photosynthesis, nitrogen metabolism and carbon fixation pathways play key roles in bermudagrass response to cold stress. The results of this study may contribute to our understanding the molecular mechanism underlying the responses of bermudagrass to cold stress, and also provide important clues for further study and in-depth characterization of cold-resistance breeding candidate genes in bermudagrass.

transcription factor [4]. Transcriptome analysis also showed that only 12 % of the cold responsive genes are controlled by CBFs [5], suggesting that there were CBFindependent components involved in cold signaling. For example, loss function of HOS9 gene encoding a homeobox transcription factor causes reduced freezing tolerance without changing the expression of CBFs and their target genes [6]. Although much progress has been made toward elucidating the molecular mechanisms of plant responses to cold stress, how plants sense low temperature signals remain unanswered. The recent findings support the hypothesis that plant cells can perceive cold stress and subsequently trigger the production of second messengers, such as Ca 2+ via membrane rigidification [7].
In recent years, the RNA-Seq has become a key technology for investigating transcriptome profiling among different species by de novo assembly or mapping. Besides, RNA-Seq is an efficient means to generate functional genomic data for non-model organisms or those with genome characteristics extremely difficult to whole-genome sequencing [8,9]. For instance, RNA-Seq has been successfully applied to characterize the transcriptomic response to low temperature in Chrysanthemum (Chrysanthemum morifolium), lily (Lilium lancifolium) and tea (Camellia sinensis) [10][11][12].
Bermudagrass [Cynodon dactylon (L). Pers.] is one of the most widely used warm-season turfgrass species for parks, lawns, and sport fields especially in golf courses [13,14]. Bermudagrass displays high tolerance to salt, drought and heat stresses, but is sensitive to cold stress [15,16]. Cold stress is a key factor limiting widespread use of bermudagrass. Thus, it is important to improve cold tolerance for bermudagrass. Although previous studies have identified several physiological and metabolic changes in bermudagrass after cold treatment, including the expression of genes encoding chitinase, dehydrin and antioxidant enzyme, protein synthesis, amino acid metabolism [15][16][17][18][19][20], the physiological and molecular mechanism of cold stress response in bermudagrass is largely unknown.
To date, few studies have been carried out to the transcriptional studies in bermudagrass. The transcriptomic responses of bermudagrass to low temperature using RNA-Seq have not been reported so far. In this study, the RNA-Seq platform based on Illumina NGS technology was used to investigate the transcriptomic response to low temperature by comparing the different transcriptome between two cold contrasting bermudagrass genotypes (Cold-resistant and -sensitive) subjected to periods of subzero temperature with or without a prior CA. Thus, the objectives of the present study were to (a) identify genes involved in response to chilling/freezing; (b) elucidate the molecular mechanisms of cold tolerance through transcriptomic analysis of the two genotypes differing in tolerance to cold stress; (c) gain a deep insight into the molecular basis of CA process in enhancing plant freezing tolerance.

Plant materials and growth conditions
The 128 bermudagrass accessions were planted in the plastic pots (15 cm diameter and 20 cm tall) filled with matrix (brown coal soil: sand 1:1). Each accession was repeated 3 times. The plants were treated with 4°C for 21 d, and the plants cultivated under 30/25°C (day/ night) were set as the control. Transpiration rate and growth rate of the plants were determined every week. The membership function method of fuzzy mathematics was analyzed using the phenotypic traits after a 21 d chilling treatment. The membership values of each accession were the index of cold tolerance. After the first round screening, 5 relatively cold-tolerant and 5 coldsensitive accessions were obtained, respectively. To further screen the relatively most cold-tolerant and coldsensitive genotypes, the 10 accessions were treated with −5°C for 4 h with or without cold acclimation. Finally, the most promising cold-tolerant (R) and -sensitive(S) bermudagrass genotypes were selected and further confirmed, respectively (Additional file 1).
The cold-tolerant (R) and -sensitive(S) bermudagrass plants were grown in plastic pots with a mix of sand and peat soil (1/1, v/v) in the greenhouse with natural sunlight, relative humidity of 87 %, and temperatures of 30/20°C (day/night). The plants in pots are ramets of the same clone, and the genetic background for these plants is uniform. After two months of establishment, plants were transferred to controlled-environment growth chambers (HP300GS-C; Ruihua Instrument, Wuhan, China), with a 14-h photoperiod, photosynthetically active radiation at 450 μmol m −2 s −1 in the canopy level with a day/night temperature of 30/20°C and 70 % humidity. Plants were fertilized three times a week with half-strength Hoagland's solution until dripping throughout the experiment in order to keep them close to field capacity.

Treatments and experimental design
When allowed to acclimate for 3 days at normal condition, plants were exposed to various cold treatments. The cold-tolerant and -sensitive genotypes were divided into two groups (Group I, II). Plants in Group I were placed in a freezing chamber set to 4°C for 48 h before being transferred to −5°C for 4 h, whereas plants in Group II without CA were directly incubated at −5°C for 4 h. The leaf samples for transcriptome sequencing were collected at 0 h (named CdR_0, CdS_0), 24 h (CdRCA_24, CdSCA_24) and 48 h (CdRCA_48, CdSCA_48) after 4°C treatment, −5°C for 4 h after prior CA (CdRCA_4, CdSCA_4), and −5°C for 4 h without prior CA (CdRNA_4, CdSNA_4), respectively. At each sampling time point, the leaves from three pots (three replicates) of each genotype were pooled together as one biological replicate and frozen immediately with liquid nitrogen, and stored at −80°C in preparation for RNA-Seq analysis. There were ten samples in total used for Illumina Genome Analyzer deep sequencing.

RNA preparation
Total RNA was isolated from the leaves using TRIzol reagent according to the manufacturer's protocol (Invitrogen, CA, USA). Then, RNA degradation and contamination was monitored on 1 % agarose gels. RNA purity was checked using the Nano Photometer spectrophotometer (IMPLEN, CA, USA). The RNA concentration was measured using Qubit RNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was evaluated using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Transcriptome sample preparation and sequencing
The total amount of 3 μg RNA per sample confirmed for RIN values above 8.0 was used as input material in constructing the sequencing library. The library was generated using Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to manufacturer's recommendations, and ten index codes were added to the sample for subsequent documentation. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in Illumina proprietary fragmentation buffer. First-strand cDNA was synthesized using random oligonucleotides and SuperScript II. Secondstrand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities and enzymes were removed. After adenylation of 3' ends of DNA fragments, Illumina PE adapter oligonucleotides were ligated to prepare for hybridization. To select cDNA fragments of preferentially 200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, MA, USA). DNA fragments with ligated adaptor molecules on both ends were selectively enriched using Illumina PCR Primer Cocktail in a 10 cycle PCR. Products were purified (AMPure XP system) and quantified using the Agilent high-sensitivity DNA assay on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded sample was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the vender's instructions. After cluster generation, the library preparation was sequenced on an Illumina Hiseq 2000 platform and 100 bp single-end reads were generated.

Bioinformatic analysis Quality control
The raw reads were processed by removing reads containing adapter, reads containing ploy-N and low quality reads, and then the clean data (clean reads) were obtained. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality.

Transcriptome assembly
The left files (read1 files) from all libraries/samples were pooled into one big left.fq file, and right files (read2 files) into one big right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity [21] with min_kmer_cov set to 2 by default and all other parameters set default.

Gene functional annotation
Gene function was annotated using the nucleotide (Nt) and protein (Nr, Pfam and Swiss-Prot) database, and assigned to functional categories in the KOG/COG, GO and KEGG database by searching BLASTx with an E value cutoff of 10 −5 .

Differential expression analysis
Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq (2010) R package. P-value was adjusted using q value [22]. q value < 0.005 & |log2(foldchange)| > 1 was set as the threshold for significant differential expression.

GO enrichment analysis
Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R packages based on Wallenius non-central hyper-geometric distribution [23], which can be adjusted for gene length bias in DEGs.

KEGG pathway enrichment analysis
KEGG [24] is a database resource for understanding highlevel functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecularlevel information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS [25] software to test the statistical enrichment of differential expression genes in KEGG pathways.

Validation of RNA-seq data by real-time quantitative PCR
To validate the expression of the candidate gene, real-time quantitative RT-PCR was employed by the method described previously by Chen et al. (2012Chen et al. ( , 2013 [26,27], and the CdACT2 gene was used as a quantitative control.

Transcriptome sequencing and assembly
To comprehensively survey the genes associated with cold stress response in bermudagrass, ten cDNA libraries were constructed from total RNA extracted from leaves of bermudagrass (cold-resistant and cold-sensitive genotypes) with different cold treatments. The libraries were sequenced using the Illumina HiSeq™ 2000 platform. An overview of the RNA-Seq reads derived from the ten libraries was presented in Because of the absence of reference genomic sequences, de novo assembly was employed to construct transcripts from these RNA-seq reads. Trinity software was used for de novo assembly of the Illumina reads, which has been proven to be efficient for de novo reconstruction of transcriptomes from RNA-Seq data [21,28]. A total of 326,435 contigs were obtained from the clean reads with a mean length of 1277 bp and length ranging from 201 bp to 20202 bp ( Table 2). Among the 326,435 contigs, 121,166 unigenes were obtained with an average length of 706 bp. The longest and shortest unigene was 20,202 bp and 201 bp, respectively (N50 was 1276 bp, N90 was 269 bp).

Gene annotation
The unigenes were annotated by searching against the seven public databases (

Differential expression genes (DEGs) analysis under various cold treatments
DEGs (q-value < 0.005 and |log2 (fold change)| >1) were defined as genes that were significantly enriched or depleted in one sample relative to the other sample. From the ten comparisons, including treatment R1 (CdRCA_24 vs CdR_0), R2 (CdRCA_48 vs CdR_0), R3 (CdRCA_4 vs CdR_0), R4 (CdRNA_4 vs CdR_0), R5 (CdRCA_4 vs CdRNA_4), S1 (CdSCA_24 vs CdS_0), S2 (CdSCA_48 vs  (Fig. 4). Further hierarchical clustering method was employed to observe the overall expression pattern of the differentially expressed genes (Fig. 5). The blue bands identify low gene expression quantity, and the red represent the high gene expression quantity. The results revealed that more DEGs were detected in comparison R1 than that in S1, suggesting that global gene expressions were initiated more quickly in R genotype than those in S genotype, when they were exposed to cold stress. In addition, more DEGs were identified in the comparisons R3 and S3, which underwent a prior cold acclimation (CA) for 48 h, as compared to the treatments which didn't undergo CA (R4 and S4), respectively (Figs. 4 and 5). It should be noted that the number of DEGs in R genotype is larger than that in S genotype, when they were subjected to freezing conditions (−5°C for 4 h) with or without CA. However, there were no obvious differences between the comparisons S3 (CdSCA_4 vs CdS_0) and S4 (CdSNA_4 vs CdS_0) from the hierarchical clustering analysis (Fig. 5). When comparing the R5 (CdRCA_4 vs CdRNA_4) and S5 (CdSCA_4 vs CdSNA_4) treatments, it was surprisingly found that R5 had 4315 DEGs (1717 upand 2598 down-regulated), whereas only 269 DEGs (127 up-and 142 down-regulated) were identified in S5 treatment. Further analysis using a venn diagram showed that both unique and overlapping sets of differentially expressed genes were detected at each treatment in both R and S genotypes (Fig. 6). Among these DEGs, 432 were categorized as commonly induced genes in R genotype comparisons, R1 (CdRCA_24 vs CdR_0), R2 (CdRCA_48 Average length 1277 706 Note: The N50 size is computed by sorting all transcripts from largest to smallest and by determining the minimum set of transcripts whose sizes total 50 % of the entire transcript and unigene was the same; N90 was counted in the similar way vs CdR_0), R3 (CdRCA_4 vs CdR_0) and R4 (CdRNA_4 vs CdR_0), while 367 were identified as overlap in four S genotype comparisons, S1 (CdSCA_24 vs CdS_0), S2 (CdSCA_48 vs CdS_0), S3 (CdSCA_4 vs CdS_0) and S4 (CdSNA_4 vs CdS_0) (Fig. 6).

Function annotation of DEGs using the KEGG database
Unigene KEGG annotation was aimed at DEGs from the above comparisons. In the R1 comparison, 1531 DEGs were assigned to the KEGG database involving 160 pathways; for R2, 1413 DEGs were assigned to 159 pathways; for R3, 1245 DEGs were assigned to 156 pathways; for R4, 914 DEGs were assigned to 125 pathways; for S1, 948 DEGs were assigned to 138 pathways; for S2, 2345 DEGs were assigned to 167 pathways; for S3, 510 DEGs were assigned to 118 pathways; and for S4, 461 DEGs were assigned to 120 pathways. The details of the KEGG classification of the above comparisons are presented in Additional file 4.

Genes involved in the response to low temperature
The Ca 2+ signaling components mainly included calciumbinding protein (CBP), calmodulin-like protein (CML), calcium-dependent protein kinase (CDPK), calcineurin Blike protein (CBL), CBL-interacting protein kinases (CIPK), and calmodulin-binding receptor like kinases (CBRLK) [29]. In the R1 comparison, there were 6 CML, 2 CBRLK, 3 calmodulin-binding protein, 2 Calciumbinding protein, 1 extracellular calcium sensing receptor, 3 CDPK, 1 CBL and 12 CIPK. The equivalent order for the R2 comparison was, respectively, one, four, three, two, one, three, one and nineteen; for R3 comparison, two, three, one, two, one, three, one and eleven; for R4 comparison, four, zero, two, one, one, zero, one and four; for S1 comparison, three, zero, two, two, zero, two, zero and seven; for S2 comparison, ten, three, six, one, one, three, one and thirteen; for S3 comparison, zero, one, one, two, zero, three, zero and four; for S4 comparison, there are only one CBP, one CDPK and three CIPK. Among these differential expression Ca 2+ signaling genes, the expression of unigene (comp148141_c0) encoding calcium binding protein was up-regulated in the comparisons R2, R3, R4, S1 and S3. The transcripts of CBP unigene (comp132952_c0) was induced in R2, S1, S2, S3 and S4. By contrast, another CBP gene expression was only induced in the cold-resistant bermudagrass genotype under cold treatment (comparisons R1 and R3). It is very interesting to find that one gene expression (comp151017_c0) encoding  One CDPK gene (comp156791_c0) transcripts were also accumulated in comparisons R1, R2, R3 and S2. The complete details of DEGs involved in Ca 2+ signalling pathway are presented in Additional file 5. The CBL-CIPK signaling networks have been proven to play important roles in response to a wide range of stimuli. Here, only two CBL genes were identified as DEGs, and both genes were up-regulated by cold treatment. Induction of expression of one CBL gene (comp151010_c0) was observed in the following comparisons R1, R2, R3 and S2. Besides, another one CBL gene (comp151988_c1) was induced in comparisons R4, showing that the gene may be involved in plant response to chilling stress without a prior CA. The number of differentially expressed CIPK genes was 46 and 27 in comparisons of cold-resistant and -sensitive genotypes of bermudagrass, respectively, revealing that more CIPK genes were involved in cold response in the cold-resistant genotype.
It was interestingly found that most of the identified CIPK genes were down-regulated by cold stress, while 7 genes identified in the S1 comparison were all up-regulated. These results revealed that expression profiles of CIPK genes were different in R and S genotypes under cold condition. The complete details of DEGs associated with CIPK are presented in Additional file 6. Similarly, DEGs associated with the MAPK cascade were twelve in comparisons of coldresistant genotype, while only seven related genes were detected in comparisons of cold-sensitive genotype. The complete details of DEGs associated with MAPK are presented in Table 4. One MAPKKK gene (comp155944_c1) was found to be down-regulated in R1, R2, R3 and S3 comparisons, implying that the gene may be specifically involved in the CA process. The expression of another MAPKKK gene (comp158986_c0) was induced in R2, R3 and R4 comparisons, and the induction folds were higher in R3 (5.26) than that in R4 (3.51)  comparisons, suggesting that the gene could be more effectively activated to respond to chilling treatment after CA process.
In the present study, members of various low temperatureresponsive transcription factor (TF) families were identified. The major TF families presented were AP2/ERF, bHLH, WRKY and NAC family. There are 7 and 6 cold upregulated genes associated with the NAC family identified in various comparisons in R and S genotypes, respectively (Table 5). Of these NAC genes, comp148886_c0 and comp150085_c0 were induced by low temperature in both R and S genotypes, but the induction folds by cold were higher in R genotype than that in S genotype.
As shown in Table 6, comp160681_c0 and comp160771_c0 encoding WRKY TF were up-regulated in the R1, R2, R3 and S2 comparisons, suggesting that these two WRKY proteins are involved in the CA process in both R and S genotypes, but specifically involved in response to freezing treatment in plants with prior exposure to CA in R genotype. Another WRKY gene (comp160978_c0) was found to be differentially transcribed in R2, R3, R4, S1 and S3, implying that it is involved in the CA process and freezing treatment in plants with prior exposure to CA in both R and S genotypes, and that this WRKY gene also plays essential roles in freezing treatment without prior CA process in the R genotype. Moreover, comp144527_c0 gene was differentially expressed in R2, R3, S1 and S3 comparisons. There are seven, ten, six and three up-regulated expressed genes encoding bHLH transcription factors identified in R1, R2, R3 and R4 comparisons of R genotype, respectively, whereas four, ten, two and two were detected in S1, S2, S3 and S4 comparisons of the S genotype (Table 7). Overall, the number of up-regulated expressed bHLH genes in R genotype was greater than that of S genotype. Of the identified bHLH genes, the log2 (fold change) of comp155113_c1 gene reached to 6.4 and 4.98 in R2 and S2 comparisons, respectively, but it was not detected as DEGs in other comparisons. In addition, the bHLH gene (comp151458_c0) expression was largely upregulated in R1, R2 and R3 comparisons, while it was only found to be increased in the S2 comparison (Table 7).
CBF TFs belong to the AP2/ERF (APETALA2/ethyleneresponsive factor) superfamily. In our present study, comp133037_c0 gene encoding CBF3 TFs was found to be specifically and highly expressed in R genotype (R2 and R3 comparisons), and the log2 folds were increased to 7.04 and 5.99, respectively. Two CBF genes (comp143318_c0 and comp155879_c0) were commonly expressed in R and S genotypes (Table 8). In our present study, we identified low temperature sensing and signaling related genes and transcription factors as DEGs under different cold treatments. In addition, various functional proteins, such as LEA proteins and dehydrins, also accumulated under cold conditions.

KEGG pathway enrichment analysis for DEGs
The top-five enriched pathways by DEGs in comparison R1 were photosynthesis, photosynthesis -antenna proteins, nitrogen metabolism, carbon fixation pathways in prokaryotes and carotenoid biosynthesis (Additional file 7). The rich factor for the above five pathways was 49.00, 46.00, 31.43, 28.57, and 30.43 %, respectively, while an equivalent order for the S1 comparison was 16.98, 19.    25 were induced and only one was down-regulated. In contrast, only 9 genes were found to be up-regulated in the S1 comparison. Likewise, there are 10 up-regulated expressed genes along with one down-regulated gene in the photosynthesis pathway in the R3 comparison, whereas it was found that only two and one were identified as up-and downregulated genes, respectively (Fig. 8). Moreover, it was unexpectedly discovered that 23 photosynthesis related genes were all down-regulated by freezing temperature without prior CA in R genotype.

Validate the DEGs by real-time RT-PCR analysis
To validate the expression data obtained from RNAsequencing, 10 DEGs were selected for confirmation by real-time RT-PCR analysis. The qRT-PCR results showed a strong correlation with the RNA-seqgenerated data (Pearson correlation coefficients r = 0.878; Additional file 8).

Global patterns of transcription in response to low temperature
The data available on the molecular basis of the bermudagrass response to low temperature is very limited. In recent years, the development of novel high-throughput sequencing has provided an opportunity to identify cold-related genes in different species by de novo assembly or mapping, thereby contributing to elaborate the molecular mechanism underlying the response to low temperature [8][9][10][11][12]. In the present study, ten bermudagrass cDNA libraries were con- Overall, the number of DEGs in the R genotype was larger than that in the S genotype under various cold treatments. The results from parallel comparisons R1 and S1 revealed that more DEGs  were detected in comparison R1 (3295) than that in S1 (1793), suggesting that global gene expressions were more quickly initiated in the R genotype than those in S genotype, when they were exposed to cold stress. From the hierarchical clustering analysis, it was found that the S genotype began to appear in a similar model like CdRCA_24 when they were exposed to cold for 48 h (CdSCA_48), further supporting the hypothesis that the S genotype triggered gene expression more slowly than that of the R genotype under cold stress. Besides, there were more DEGs identified in comparisons R3 and S3, which underwent a prior cold acclimation (CA) for 48 h, as compared to the treatments which didn't undergo CA (R4 and S4), respectively. Our results further provide additional evidence supporting that plants could obtain enhanced tolerance to freezing temperature when they undergo CA process, and that the acquired resistance may be attributed to a large alteration of global patterns of gene transcription in CA process. However, there was a big difference between R and S genotypes during the CA process as indicated that R5 (CdRCA_4 vs CdRNA_4) has 4315 DEGs (1717 up-and 2598 down-regulated), but only 269 DEGs (127 up-and 142 down-regulated) were identified in S5 (CdSCA_4 vs CdSNA_4) treatment. These results suggest that the R genotype could more effectively activate gene expression during CA process, and thereby better respond to freezing temperature, as compared to S genotype. It was speculated that the prerequisite for plants to obtain enhanced tolerance to freezing through a cold acclimation process is that the plant needs appropriate cold resistance levels.

Low temperature sensing and signaling genes
It has been shown that various abiotic stresses, including cold, can trigger intracellular changes in free Ca 2+ concentration, thereby generating the so-called Ca 2+ signature, which can be sensed by Ca 2+ sensors and then transduced through the interaction with their target protein to regulate the expression of stress-responsive genes [29][30][31][32][33]. There are Ca 2+ sensor proteins of three major classes: CDPK, CaM (CML) and CBL [29,31,32]. Here, six, one, two and four CMLs were identified as significant DEGs in R1, R2, R3 and R4 comparisons, respectively. Among these 13 CMLs, 9 were found to be down-regulated, indicating that some reduced CML genes expression may contribute to enhance plant tolerance to cold stress. Not surprisingly, it was revealed that over-expression of Arabidopsis CaM3 impairs cold induction of RD29A, KIN1 and KIN2 [34]. One CML gene (comp147675_c0) was downregulated in R1, R2, and R3 comparisons. Another CML gene (comp135890_c1) transcript was decreased in R1 and R3 comparisons, while it unexpectedly increased in R4 comparison. Likewise, the expression of comp148637_c0 and comp152137_c1 CML genes was significantly down-regulated in the R1 comparison, but induced in R4 comparison. These results suggest that CML family proteins may play different roles in the CA process and freezing response with or without a prior CA. Furthermore, no CML genes were identified in S3 and S4 comparison in the S genotype, implying that CML family proteins may function as important signaling responders in conferring bermudagrass tolerance to freezing temperature.
The Arabidopsis and rice genomes harbor 34 and 29 CDPK-encoding genes, respectively [35,36]. CDPKs have been identified as being involved in cold signaling. OsCPK7/OsCDPK13 is activated by cold treatment [37] and overexpression of either OsCPK7/OsCDPK13 or OsCPK13/OsCDPK7 improves cold tolerance in transgenic rice [37,38]. In the present study, there were many CDPKs identified as DEGs in each comparisons. One CDPK gene (comp156791_c0) transcript is accumulated in comparisons R1, R2, R3 and S2, suggesting that the gene may be not only involved in CA process in both R and S genotypes, but also exclusively involved in the freezing response through prior CA in the R genotype in bermudagrass.
The CBL proteins are characterized as a group of plant calcium sensors that could exclusively interact with CIPK proteins. The CBL-CIPK signaling components constitute a specific regulatory network of Ca 2+ signaling in plant cells [39][40][41]. Many CBLs and CIPKs have been identified as being involved in plant responses to various abiotic stresses. However, there are few reports on the CBL-CIPK involved in cold stress responses to date. Previous studies have revealed that AtCBL1 is involved in cold response [42,43]. CIPK3 and CIPK7 were reported to be involved in response to cold stress [39,42]. Here, two CBL genes were identified as DEGs. One CBL gene (comp151010_c0) was induced in the following comparisons, R1, R2, R3 and S2. Another CBL gene (comp151988_c1) was induced in R4 comparisons, suggesting that the gene may be involved in plant response to chilling stress without a prior CA. There are 46 and 27 DEGs encoding CIPK found in comparisons of cold-  resistant and -sensitive genotypes of bermudagrass, respectively. These differentially cold-regulated CBL-CIPK components may be useful for breeding coldresistant bermudagrass in the future.
Major classes of TF involved in the response to low temperature It has been well established that transcription factors (TFs) play important roles in response to different abiotic and biotic stresses. Here, members of various low temperatureresponsive transcription factor (TF) families were identified as DTGs in the treatments involved in a process of low temperature acclimation or freezing with or without CA process. NAC (NAM, ATAF, and CUC) is a plant specific transcription factor family with diverse roles in plant development and in response to abiotic stresses [44][45][46][47][48][49]. Hu et al. (2008) reported that over-expression of a stress-responsive NAC gene, SNAC2, increases rice tolerance to cold and salt [50]. Similarly, overexpression of TaNAC2 resulted in enhanced tolerances to salt, drought and freezing stresses in Arabidopsis [51]. More recently, Banana MaNAC1 was proven to be a direct target of MaICE1 and involved in cold stress through interacting with MaCBF1 [52]. Here, comp148886_c0 and comp150085_c0 were found to be induced by cold at higher folds in the R genotype than that in the S genotype, suggesting that the two genes may play an important function in conferring R genotype with enhanced cold tolerance, and should be the focus of future studies in bermudagrass.
WRKY TFs are a large family of regulators involved in various developmental and physiological processes, especially in response to diverse biotic and abiotic stresses. Recently, the results from high-throughput transcriptomic analyses have identified that 61 of the Populus WRKY genes were induced by various biotic and abiotic treatments, including cold [53]. Transgenic rice overexpressing OsWRKY76 led to drastically increased susceptibility to M. oryzae, but enhanced tolerance to cold stress [54]. In the present study, two WRKY genes, comp160681_c0 and comp160771_c0, were induced in the R1, R2, R3 and S2 comparisons, suggesting that these two WRKY proteins play essential roles in the CA process in both R and S genotypes, but specifically involved in the response to the freezing treatment in plants with prior exposure to CA in the R genotype, suggesting its key roles in conferring bermudagrass enhanced freezing tolerance in the R genotype after CA, compared to the S genotype. Moreover, one WRKY gene (comp160978_c0) was examined to be significantly upregulated in the R4 comparison, but not induced in the S4 comparison, indicating that its distinctive function was involved in the R genotype response to chilling without the CA process. These results provide new information that the identified WRKY genes from coldresistant bermudagrass may serve as a target gene for breeding new varieties in future.
To date, plant bHLH TFs have been demonstrated to function as transcription regulators involved in a diversity of biological processes, including flowering [55], trichome development [56], root hair formation and development [57,58], nodule vascular patterning [59] and the photoinduced signal transduction [60]. Furthermore, bHLH TFs are involved in the plant response to various abiotic stresses, including drought [61], cold [4,62,63] and iron deficiency [64]. However, although only a few bHLHs have been identified to be involved in cold tolerance mainly in model plant, cold responsive bHLHs needs further identification, and the underlying mechanisms need further elucidation. In recent years, using RNA-Seq and digital gene expression (DGE) technologies, low temperature induced bHLH genes have been identified in grape (Vitis amurensis and Vitis vinifera) [65], Chrysanthemum (Chrysanthemum morifolium) [11] and tea (Camellia sinensis) [10]. Here, we also identified cold induced bHLHs in bermudagrass. These genes may play important roles in the enhanced cold hardiness of bermudagrass and should be the focus of future studies in bermudagrass. CBF TFs have been best proven to play primary roles in response to cold stress. Six CBFs have been identified in Arabidopsis, and three of them, namely, CBF1/DREB1B, CBF2/DREB1C, and CBF3/DREB1A, have been proved to play a primary role in cold acclimation [66][67][68][69][70][71][72]. In our present study, CBFs specifically and highly expressed in R genotype were identified, and the results contributed to the understanding of the mechanism of bermudagrass response to cold.
Enriched KEGG pathway participated in response to low temperature Photosynthesis is regarded as a highly integrated and regulated process which is highly sensitive to environmental changes, because it needs to balance between the light energy absorbed by the photosystems and the energy consumed by metabolic sinks of the plant [73]. It is clear that optimum plant growth and development require a balance in the rates of source versus sink processes. However, cold stress can lead to an imbalance between the source of energy and the metabolic sink, thus requiring photosynthesis adjustments to maintain the balance of energy flow [74]. Previous studies have also reported that low-temperature modulation of the photosynthetic apparatus may be an important factor during the induction of freezing resistance in cereals [75]. As described above, our results revealed that the R genotype may better respond to chilling and freezing, with prior CA, through activating photosynthesis pathway related gene expression, in contrast to the S genotype.