The DX genotype has a better cold tolerance compared with GN genotype
To elucidate the physiological responses of both genotypes exposed to cold stress, biomass and several important physiological indices were measured. Cold stress relatively reduced the plant fresh weight of GN but not of DX (Additional file 1: Figure S1a). Concurrent with growth inhibition, there was a lower relative water content in GN than DX exposed to 24 h and 5 d of cold treatment (Additional file 1: Figure S1b). Malondialdehyde (MDA), a product of lipid peroxidation, considered as an indicator of oxidative injury under abiotic and biotic stress [18]. MDA content and electrolyte leakage levels increased (P < 0.05) in both plants under 24 h and 5 d of cold stress, while a smaller increase in these values was observed in DX (Additional file 1: Figure S1c, d), indicating lower oxidative damage in DX under cold stress. DX accumulated (P < 0.05) higher proline than GN in both control and stressed conditions (Additional file 1: Figure S1e). These results suggest that the DX genotype was more cold-tolerant than GN genotype.
Transcriptome assembly and annotation
To further elucidate the molecular mechanisms underlying differential cold tolerance in DX and GN, twenty-four cDNA libraries were created from three independent biological samples of four treatments (0, 3, 24 h, and 5d) and sequenced using the Illumina HiSeq™ 4000 sequencing platform. More than 65 M high-quality pair-end reads were retrieved after trimming for each RNA-seq sample. Sequencing was done on the Illumina platform generating paired end reads of 100 bp each. Cleaned reads were de novo assembled using Trinity software: a total of 200,520 and 181,331 transcripts with a N50 of 1809 and 1777 bp were generated in DX and GN genotypes, respectively (Additional file 2: Table S1). This indicates a high quality assembly. The size distribution of transcripts was shown in Additional file 3: Figure S2.
To validate and annotate the of assembled unigenes, using E-value < 1e-5, they were blast searched against seven public databases, including NR, NT, COG, GO, KEGG, Swissprot, and Interpro protein database. In total, 164,827 (82.2 %) and 150,545 (83.02 %) unigenes in the DX and GN assemblies were found in at least one of these databases (Additional file 4: Table S2). Overall, the unigene sequences were most similar to gene sequences from Hordeum vulgare subsp. vulgare, Aegilops tauschii, Triticum urartu, Brachypodium distachyon via BLASTx matches (Additional file 5: Figure S3).
To characterize the functional classifications of annotated unigenes, GO and KEGG analyses were performed to access the distributions of functional categories. A total of 62,654 (31.29 %) and 55,886 (30.82 %) unigenes were annotated in GO for DX and GN, respectively, and classified into 57 functional groups, including 23 groups in biological process, 17 in molecular function, and 17 in cellular components (Additional file 6: Figure S4). In total, 82,192 (40.99 %) unigeness in DX and 76,364 (42.11 %) in GN were assigned to 20 KEGG pathways. These unigenes were mainly involved in ‘Translation’, ‘Lipid metabolism’, and ‘Transport and catabolism’ (Additional file 7: Figure S5).
Gene expression under cold stress
Differential gene expression was analyzed relative to a control (0 h) grown under control conditions. To characterize the differentially expressed genes (DEGs) over the 5 d of cold exposure, a total of 5436 and 4323 DEGs (P-value ≥ 0.8 and |log2 (fold change)| >2) were identified and analyzed for the DX and GN genotypes, respectively (Fig. 1a; Additional file 8: Table S3). Most of the cold regulated DEGs are late-response genes. In DX, 2413 (44.4 %) unigenes were induced/repressed exclusively at 5 d of cold treatment, while only 1030 (18.9 %) and 569 (10.9 %) unigenes were specific at 3 or 24 h of cold treatment. A similar change in DEGs over time was observed in GN, starting from 1580 DEGs after 3 h to 2078 DEGs after 5 d of cold stress (Fig. 1a). For all comparisons in pairs, we identified 5350 DEGs in DX and 4237 in GN, and only 86 DEGs were commonly regulated for both genotypes (Fig. 1b), indicating that the cold-induced transcriptomic responses were largely genotype specific. A total of 4074 DEGs (P-value ≥ 0.8 and |log2 (fold change)| >2) were identified between the treated and control DX and GN plants (Additional file 9: Table S4). A total of 94 and 76 DEGs in the DX and GN genotypes represented 52 common response genes (8 persistently upregulated) at all time points. While 2451 genes were specifically regulated for each time point in DX and GN, illustrating a high species-specific transcriptomic plasticity in response to cold stress (Fig. 1c; Additional file 10: Table S5; Additional file 11: Table S6; Additional file 12: Table S7; Additional file 13: Table S8).
A hierarchical cluster analysis (HCL) was performed to characterize the expression patterns of the common expressed DEGs after cold exposure. A large proportion of DEGs in DX and GN followed a similar expression pattern after 24 h and 5 d of cold treatment, which clustered together. In contrast, expression profile of 3 h-exposed plants clustered separately, this indicated obvious differences in the global gene expression patterns at early and later phases of cold stress in the two genotypes (Fig. 1d).
Pathways enrichment analysis of DEGs
Using the KEGG database, pathways displaying significant changes (Q value ≤ 0.05) in response to the cold treatment were identified in both genotypes (Additional file 14: Table S9). In the 3 h cold treatment sample, 21 and 16 KEGG pathways were significantly enriched in DX and GN, respectively. Among the top 5 enriched pathways, ‘mRNA surveillance pathway (ko03015)’, ‘circadian rhythm – plant (ko04712)’, and ‘RNA transport (ko03013)’ were common regulated in response to 3 h cold stress in both genotypes. The intermediate and late phases of the cold response were characterized by protective response through modulating cellular metabolic homeostasis. In the 24 h cold treatment sample the enriched pathway were similar to those at 5 d in both genotypes. In DX, ‘phenylpropanoid biosynthesis (ko00940)’, ‘cyanoamino acid metabolism (ko00460)’, and ‘galactose metabolism (ko00052)’ were the top three enriched pathways after 24 h and 5 d of cold stress. In GN, ‘glutathione metabolism (ko00480)’, ‘flavonoid biosynthesis (ko00941)’, ‘biosynthesis of secondary metabolites (ko01110)’, ‘phenylpropanoid biosynthesis (ko00940)’, and ‘flavone and flavonol biosynthesis (ko00944)’ were the most significantly enriched.
Early gene-expression changes in response to cold stress
To gain insight into the functional categories of common and genotype-specific DEGs induced by 3 h of cold stress, GO term annotations were analyzed. For genes that were detected in the two genotypes at 3 h, the top enriched GO categories were invovled in ‘transcription factors (TFs)’ and ‘signal transduction’. Numerous genotype-specific TFs and signaling molecules were detected during early cold stress, indicating that DX and GN activated the downstream defense response through different signal transduction and transcriptional regulation pathways.
Transcription factor plays a key role in the regulation of upstream cold signal transduction which are capable of activating a cascade of downstream gene transcript [3]. A total of 118 genes encoding TFs belonging to 11 families were differentially expressed in both genotypes after 3 h of cold exposure (Additional file 15: Table S11). These TFs were asssociated with response to abiotic and biotic stresses and the regulation of developmental processes, including AP2-EREBP, basic helix–loop–helix (bHLH), MYB-related, MADS, NAC, basic leucine zipper (bZIP), cysteine-2/histidine-2 zinc finger proteins (C2C2)-Dof, and heat shock transcription factors (HSF). CBFs were the most represented subfamily with 40 members, including five DX-specific and 19 GN-specific genes. Additionally, other TFs such as bZIP, bHLH, NAC were cold-regulated in the two plants (Fig. 2). MADS, Sigma70-like, and Alfin-like subfamliy TFs were uniquely expressed in DX (Additional file 15: Table S11).
When suffering from unfavorable conditions, plants could trigger multiple stress-responsive signal transduction pathways that activate gene transcription and the downstream physiological adaptation [19]. Two DX-specific Ca2+ signaling-related genes, glutamate receptor 2.7 (GLR2.7) and calmodulin-binding transcription activator 4 (CMTA4), were induced during 3 h of cold stress and two GN-specific, glutamate receptor 3.4 (GLR3.4) and CBL-interacting protein kinase 7 (CIPK7) were induced during this time frame. Additionally, the up-regulation of calmodulin-like protein 3 (CML3) was higher in DX, suggesting CML3 might contribute to higher cold tolerance in DX (Additional file 16: Table S12). Protein phosphorylation and dephosphorylation have been implicated in cold signal transduction, including 75 genes annotated as encoding protein kinases and phosphatases in both genotypes (Additional file 17: Table S12). Ten genes related to phytosulfokine receptor 2 were over-accumulated in both genotypes, indicating that phytosulfokine receptor signaling may play an important role in regulating the early cold signalings in E. nutans.
We also observed metabolic changes after 3 h of cold stress, indicated by the GO term ‘biosynthetic process’ (GO:0009058). For example, early changes involved in secondary metabolism processes, such as the induction of genes encoding flavonoid and lignin biosynthesis (CL3209.Contig3_All, CL8645.Contig8_All, Unigene11053_All, Unigene39201_All, etc) were detected after 3 h (Additional file 10: Table S5; Additional file 11: Table S6; Additional file 12; Table S7). Furthermore, lipid metabolism changed during cold stress, indicated by up-regulation of genes coding wax, suberin (CL14480.Contig22_All, CL1381.Contig2_All, CL2649.Contig4_All, etc), and fatty acid biosynthesis (Unigene82701_All, CL11835.Contig5_All). Genes coding sucrose: sucrose 1-fructosyltransferase (1-SST) indicated the need for the regulation of fructan biosynthesis. Moreover, the early response encompassed the accumulation of numberous cold responsive protein-coding genes (CL10923.Contig3_All, CL6160.Contig2_All, and Unigene38511_All, etc.).
In order to further identify specific metabolic pathways differentially affected by early cold stress in DX and GN leaves, a KEGG enrichment analysis (Q value ≤ 0.05) was performed. Several different metabolic pathways were affected by early cold stress in both genotypes, including mRNA surveillance pathway (ko03015)’, ‘circadian rhythm – plant (ko04712)’, ‘RNA transport (ko03013)’, and ‘biosynthesis of secondary metabolites (ko01110)’ (Additional file 18: Table S10). In the class of biosynthesis of secondary metabolites, most enzymes of stilbenoid, diarylheptanoid and gingerol biosynthesis (ko00945) and flavonoid biosynthesis (ko00941) were specifically induced in both genotypes. In addtion, seven genes of phenylpropanoid biosynthesis (ko00940) were specifically induced in DX in response to early cold stress. Specifically, most genes of starch and sucrose metabolism (ko00500), galactose metabolism (ko00052), and cutin, suberine and wax biosynthesis (ko00073) were induced in DX during the early cold stress (Additional file 18: Table S10).
DEGs in response to intermediate and late phases of the cold stress
The intermediate and late phases of the cold response were characterized by increased lipid metabolism, oxidative stress, carbohydrate metabolism, secondary metabolism, and photosynthetic process.
‘Lipid metabolism process’ (GO:0006629) and ‘lipid transport’ (GO:0006869) were significantly induced in both genotypes during the intermediate and late stress treatment. A total of 99 lipid metabolism-related genes were identified, most of them involved in oxylipin, wax, cutin, suberin, sterol, and fatty acid biosynthetic and metabolic processes (Fig. 3; Additional file 16: Table S13). Genotype-specific regulation was prominent, as different genes—albeit with similar functions—were affected in the two genotypes. Genes encoding wax biosynthesis were strongly up-regulated at 24 h and 5 d in DX. Three oxylipin biosynthesis-related genes exhibited higher expression in DX, indicating oxylipin might contribute to better cold tolerance in DX. Notably, an approximately 10 fold increase in the gene encoding fatty acid desaturase (CL3006.Contig17_All) was observed after cold treatment at 5 d in DX (Additional file 16: Table S13).
Cold stress resulted in a higher reactive oxygen species (ROS) accumulation and lower antioxidant enzyme activities in GN compared with DX (Additional file 19: Figure S6), which might contribute to lower oxidative damage in DX. To further characterize the differences in ROS-related gene expression between tolerant and sensitive genotypes, cold-responsive genes related to ROS were analyzed. In all, 125 genes coding ROS producing and detoxification, including 29 in DX, 31 in GN and 65 in both genotypes (Additional file 20: Table S14), were identified as having significantly different expression profiles after late cold treatment. Consistent with higher production of ROS in GN, several genes coding respiratory burst oxidase protein and polyamine oxidase (CL20776.Contig4_All, CL20776.Contig46_All, Unigene5932_All), invovled in the generation of ROS, were more strongly induced by the late cold stress in GN than in DX (Fig. 4). Eighty-three DEGs were identified as encoding enzymes associated with ROS scavenging, including peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX), monodehydroasorbate reductase (MDAR), glutathione S-transferase (GST), glutathione peroxidase (GPX), polyphenol oxidase (PPO), glutaredoxin, thioredoxin. Among the antioxidant enzymes we identified, GST and POD play the most important roles in scavenging ROS in both plants. Moreover, the induction of several genes involved in polyamine and tocopherol biosynthesis in both genotypes, were also noteworthy (Fig. 4; Additional file 20: Table S14).
Carbohydrate metabolism was affected by cold stress in both genotypes (Additional file 21: Table S15). Fructan plays important roles in response to cold stress, corroborated by the induction of 1-SST and differential expression of fructan 6-exohydrolase (6-FEH) in both plants (Fig. 5a, b). The fructan content was increased around two-fold in GN and five-fold in DX during 5 d of cold treatment. Consistent with increased concentrations in sucrose, fructose, and raffinose, enzymes such as sucrose synthase 1, beta-amylase, alpha-amylase isozyme C, and galactinol synthase 2 (Fig. 5; Additional file 21: Table S15), were induced in both plants. Additionally, cold stress increased the abundance of alpha, alpha-trehalose-phosphate synthase (TPS6) transcript by 5.2-fold in DX at 5 d under cold stress.
The ‘secondary metabolism’ category also showed genotype-specificity in both plants. In all, 110 genes encoding secondary metabolism biosynthesis, including 44 DX specificity, 34 GN specificity and 32 commonly expressed (Additional file 22: Table S16). Flavonoids and stilbenes have ROS-scavenging activity that protects against oxidative damage and controls ROS levels under abiotic stresses [20]. Most of the flavonoid biosynthesis pathway responded late in the cold, with higher expression in DX than in GN. In the stilbenes biosynthesis pathway, one putative O-methyltransferase 2 (OMT2; Unigene44864_All) was upregulated by 9.6-fold in DX and trans-resveratrol di-O-methyltransferase (ROMT; CL13661.Contig2_All) was induced by 7.9-fold in GN. In DX tropane alkaloid, caffeine, menthol and paclitaxel biosynthesis were induced uniquely, this contrasted with GN in which benzylisoquinoline alkaloid biosynthesis was exclusively depressed. Lignin metabolism was significatly enriched in both genotypes. Altered pathways included, for example, the up-regulation of lignin catabolism and down-regulation of ligin biosynthesis in GN (Fig. 6).
Numerous genes encoding components of photosystem I (PSI) and PSII declined significantly during the intermediate and late phases of cold treatment. Over-represented GO terms describing sets of declining genes after 24 h and 5 d encompassed: ‘photosynthesis’ (GO:0015979), chloroplast thylakoid membrane (GO:0009535), and photosynthesis, light harvesting (GO:0009765) (Additional file 23: Table S17). Additionally, six genes encoding chlorophyll a biosynthesis significantly decreased in both genotypes, including protochlorophyllide reductase (POR) and glutaminyl-tRNA synthetase (HEMA). Indeed, chlorophyll content had a marked reduction in both plants under cold stress, showing a greater reduction in GN compared to DX (Fig. 7a). Transcriptome analysis showed that PSII was inhibited more severely in GN (Additional file 23: Table S17). To confirm this, PAM was used to determine the maximum quantum efficiency of PSII (Fv/Fm), apparent electron transport rate (ETR), and non-photochemical quenching (NPQ). During cold stress, Fv/Fm and ETR decreased, whereas NPQ increased in both plants. DX maintained higher Fv/Fm, ETR, and lower NPQ than GN under cold stress (Fig. 7).
Genes involved in specific metabolic pathways and cellular response processes in both genotypes were further analysed using KEGG enrichment analysis. The intermediate and late phases of the cold response were characterized by protective response through modulating cellular metabolic homeostasis, which is reflected by the largest numbers of metabolic pathways (ko01100) in both genotypes in response to the prolonged cold stress. A variety of different biochemical pathways were affected by cold stress in both plants, including cartbohydratre metabolism, energy generation, lipid metabolism, and secondary metabolism (Additional file 18: Table S10). In the category of lipid metabolism, genes of sterol biosynthetic process (ko00100), alpha-linolenic acid metabolism (ko00592), and most cutin, suberine and wax biosynthesis (ko00073) were induced in DX during the late cold stress, while genes of alpha-linolenic acid metabolism pathway in GN were down-regulated. Induction of glycolysis/ gluconeogenesis pathways (ko00010) in DX leaves is essential to avoid energy starvation caused by reduced photosynthesis under late cold stress. Flavonoid biosynthesis (ko00941) and stilbenoid, diarylheptanoid and gingerol biosynthesis (ko00945) were significantly enriched in the DX and GN leaves exposed the prolonged cold stress. Majority genes involved in these processes were up-regulated in both genotypes, except genes of flavonoid biosynthesis in 24 h-treated GN leaves. The transcriptomes of the DX and GN leaves further displayed an enrichment of the category protein processing in endoplasmic reticulum (ko04141) in response to late cold stress (Additional file 18: Table S10).
ABA metabolism and signalling in response to cold stress
ABA is an important stress hormone that mediates abiotic stress responses in plants [21]. The genes involved in the ABA response increased over time in both genotypes. This is supported by the over-represented GO term ‘response to abscisic acid stimulus’ (GO:0009737) and ‘abscisic acid mediated signaling pathway’ (GO:0009738) seen after 5 d of cold treatment, especially in DX genotype (Additional file 24: Table S18). Six ABA synthesis-related genes, encoding zeaxanthin epoxidase and a probable aldehyde oxidase, were induced in both genotypes. ABA catabolism contained genes abscisic acid 8'-hydroxylase 3 (CYP707A7) and abscisate beta-glucosyltransferase (AOG) were exclusively induced in GN, while abscisic acid 8'-hydroxylase 1 (CYP707A5) declined only in DX. We further meaured the ABA content in both plants. ABA levels immediately increased in both genotypes in response to cold, with DX showing this effect to a greater extent (Fig. 8).
Dehydrins function as candidate genes regulating cold tolerance
Numerous COR genes were induced during the cold stress in both genotypes (Additional file 25: Table S19). Most COR genes were expressed at common levels of expression, indicating these COR genes were conserved in regulating down-stream defense responses in both genotypes. Four genes encoding dehydrin (2, 3, and RAB15) and cold shock protein (CSP1) over-accumulated only in DX. While two late embryogenesis abundant protein-related genes (LEA-like, LEA1) were induced exclusively in GN. DHN4 was induced in DX but repressed in GN at the later stages of cold stress. Genes encoding LEA14-A, cold-regulated 413 plasma membrane protein (COR413PM), dehydrin (COR410), and cold-responsive protein (COR14a) had higher transcriptional abundance in DX compared with GN, indicating those CORs may contribute to the enhanced cold tolerance in DX. In contrast, DHN5 and LEA3 were much stronger over-represented in GN than in DX (Additional file 25: Table S19).
Co-expression network and regulatory interactions
In order to further clarify the regulatory network induced by cold stress in both genotypes, a co-expression network analysis was performed using Pearson correlation thresholds [22]. In total, 98 and 140 genes formed a cold-responsive co-expression network in DX and GN, respectively (Additional file 26: Figure S7; Additional file 27: Table S20). These networks were directly connected with each other via 1632 edges in DX and 1288 edges in GN, suggesting that the coexpressed genes are most likely coregulated. Twenty-six of these coexpressed genes (22 in DX, 4 in GN) were strongly interconnected, with each gene having more than 50 edges (Additional file 26: Figure S7). These genes were therefore defined as hub genes. These hub genes were directly connected to each other via 428 edges in DX and 12 edges in GN, forming a highly interconnected subnetwork. Twenty-two hub genes (84.6 %) were annotated as COR (COR413PM1, COR410, DHN5, etc.), reflecting their central roles in the acquisition of cold tolerance.
Validation of gene expression profiles by qRT-PCR
To confirm the accuracy and reproducibility of this Illumina RNA-seq result, 26 DEGs including 10 hub genes were chosen for qRT-PCR. These genes were from various functional categories include signaling, transcript factor, metabolism, and cold-responsive protein. The primer sequences, FPKM and qRT-PCR result are listed in Additional file 28: Table S21 and Additional file 29: Figure S8. All 26 genes in the two genotypes showed the same expression patterns in the qRT-PCR assays as in the RNA-Seq data (Additional file 29: Figure S8). This independent evaluation revealed the reliability of the RNA-seq data.