Dynamic changes in the transcriptome landscape of Arabidopsis thaliana in response to cold stress

Abstract

Plants must reprogram gene expression to adapt constantly changing environmental temperatures. With the increased occurrence of extremely low temperatures, the negative effects on plants, especially on growth and development, from cold stress are becoming more and more serious. In this research, strand-specific RNA sequencing (ssRNA-seq) was used to explore the dynamic changes in the transcriptome landscape of Arabidopsis thaliana exposed to cold temperatures (4°C) at different times. In total, 7,623 differentially expressed genes (DEGs) exhibited dynamic temporal changes during the cold treatments. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that the DEGs were enriched in cold response, secondary metabolic processes, photosynthesis, glucosinolate biosynthesis, and plant hormone signal transduction pathways. Meanwhile, long non-coding RNAs (lncRNAs) were identified after the assembly of the transcripts, from which 247 differentially expressed lncRNAs (DElncRNAs) and their potential target genes were predicted. 3,621 differentially alternatively spliced (DAS) genes related to RNA splicing and spliceosome were identified, indicating enhanced transcriptome complexity due to the alternative splicing (AS) in the cold. In addition, 739 cold-regulated transcription factors (TFs) belonging to 52 gene families were identified as well. This research analyzed the dynamic changes of the transcriptome landscape in response to cold stress, which reveals more complete transcriptional patterns during short- and long-term cold treatment and provides new insights into functional studies of that how plants are affected by cold stress.

Keywords: Arabidopsis thaliana; alternative splicing; cold stress; dynamic change; long non-coding RNA; ssRNA-seq; transcriptome landscape.

FIGURE 1

The overall analysis of the transcriptomes of FRI-Col in different times of cold treatment. (A) Pearson’s correlation coefficient of gene expression values of all samples. ***, P < 0.001. (B) Hierarchical clustering analysis (HCA) of gene expression values of all samples. (C) t-Distributed Stochastic Neighbor Embedding (t-SNE) of gene expression values of all samples. (D) Principal component analysis (PCA) of gene expression values of all samples. Axes show principal components (PC) 1, PC2 which explained 60.2% of the variance. Each color represents one time point and three biological replicates were performed for each time point in separate experiments (15 samples in total) in (A–D).

FIGURE 2

Changes of the transcriptome during cold treatment in Arabidopsis. (A) Bar graph showing total numbers of differentially upregulated (red) and downregulated (blue) genes at each pairwise time point. (B) Hierarchical clustering and heat map of Arabidopsis DEGs and Key GO Terms. The groups, Z-score values, and Clusters were assigned different colors. The gene number of the different clusters was shown on circle size. Bar plot of –log10 transformed FDR values are shown. BP, biological process; MF, molecular function; CC, cellular component.

FIGURE 3

Alternative splicing of the transcriptomes of FRI-Col in different times of cold treatment. (A) Identified differentially alternatively spliced (DAS) genes in the Arabidopsis transcriptome during cold treatment. (B) Summary table of AS events detected by the IsofromSwitchAnalyzeR. Exon skipping (ES), Mutually exclusive exon and multiple Exon Skipping (MXE), Alternative 5’ splice site (A5), Alternative 3’ splice site (A3), and intron retention (IR). (C) Stacked bar charts showing the AS events percentage at each pairwise time point. Line charts showing the ES, IR, A5, and A3 events in (D–G) respectively.

FIGURE 4

Expression dynamics reveal a temporal ordering of biological processes during cold treatment. (Top) Hierarchical clustering and heat map of Mfuzz cluster core expression profiles. (Bottom) GO enrichment analysis with clusterProfiler.

FIGURE 5

Construction and identification of the co-expression modules associated with the timepoints, based on the FPKM values of DEG and DAS genes. (A) WGCNA cluster dendrogram showing co-expression modules identified in 15 samples of 5 timepoints. Modules are assigned by distinct colors, respectively. (B) The number of DEGs and DElncRNAs in different Modules. (C) Correlations between modules and cold treatment timepoints. The correlation is estimated by the Pearson correlation coefficient method. (D) Predicted regulatory networks based on DElncRNAs and their potential targeted genes in DEG and DAS genes. Different colors represent different modules and the big dots represent lncRNAs. The detailed information was showing in Supplementary Table 16.

FIGURE 6

The regulatory networks of differently expressed TFs during cold treatment. (A) Bar plot showing the distribution of the top 20 differently expressed TFs in total DEG and DAS genes. (B) Bar plot showing the distribution of the top 5 of differently expressed TFs in DEG and DAS genes at each time point. (C) Enrichment analysis of TF DNA binding motifs within the promoter regions of DEG and DASs genes at each time point. (D) Predicted regulatory networks based on enriched TFs and their potential targeted genes in DEG and DAS genes. Different colors represent different TFs and the big dots represent TFs. The detailed information was showing in Supplementary Table 18.