Through vernalization, plants achieve flowering competence by sensing prolonged cold exposure (constant exposure approximately 2-5 °C). During this process, plants initiate defense responses to endure cold conditions. Here, we conducted transcriptome analysis of Arabidopsis plants subjected to prolonged cold exposure (6 weeks) to explore the physiological dynamics of vernalization and uncover the relationship between vernalization and cold stress.
Time-lag initiation of the two pathways and weighted gene co-expression network analysis (WGCNA) revealed that vernalization is independent of cold acclimation. Moreover, WGCNA revealed three major networks involving ethylene and jasmonic acid response, cold acclimation, and chromatin modification in response to prolonged cold exposure. Finally, throughout vernalization, the cold stress response is regulated via an alternative splicing-mediated mechanism.
These findings illustrate a comprehensive picture of cold stress- and vernalization-mediated global changes in Arabidopsis.
Plants are sessile organisms that passively sense environmental signals, such as temperature. Vernalization is a process through which plants achieve flowering following prolonged cold exposure. Typically, winter annual accessions of Arabidopsis require several weeks of cold exposure before flowering, whereas its rapid-flowering accessions do not require vernalization. The difference between the two accessions is determined by the expression of the dominant allele of FRIGIDA (FRI) [1, 2]. FRI encodes a scaffold protein that functions as an activator of FLOWERING LOCUS C (FLC). FLC encodes a MADS-BOX transcription factor that functions as a suppressor of the floral integrators FT and SUPPRESSOR OF OVEREXPRESSION OF CO 1 (SOC1) .
Before vernalization, FLC chromatin is in an active transcription state during vegetative growth. FRI forms a complex with FRI-LIKE 1 (FRL1), FRI ESSENTIAL 1 (FES1), SUPPRESSOR OF FRI 4 (SUF4), and FLC EXPRESSOR (FLX) to recruit other transcription factors and chromatin modifiers, ultimately activating FLC . Upon cold exposure, FLC suppression is initiated via upregulation of the noncoding FLC antisense transcript COORAIR . Subsequent suppression is realized by dynamic replacement of active markers via trimethylation of lysine 27 of histone H3 (H3K27me3) by Polycomb repressive complex 2 (PRC2) . PRC2 is conserved in both animals and plants Homologs of the Drosophila H3K27 methyltransferase E(z), namely CURLY LEAF (CLF) and SWINGER (SWN), together with VERNALIZATION 2 (VRN2), FERTILIZATION INDEPENDENT ENDOSPERM (FIE), and SUPPRESSOR OF IRA 1 (MSI1) constitute the core components of PRC2 [7, 8]. VIN3, encoding a chromatin-remodeling plant homeodomain (PHD) finger protein, forms a heterodimer with the paralog VIN3-like 1 (VIL1)/VERNALIZATION 5 (VRN5) and joins the PRC2 core to serve as the cold-specific PHD (VIN3)-PRC2. Moreover, Polycomb partners VAL1 and VAL2 serve as epigenome readers that recognize the cis-regulatory element at the FLC locus to recruit histone deacetylases HDA9 and PRC2. The former catalyzes H3K27ac deacetylation to H3K27, and the latter catalyzes H3K27 trimethylation to H3K27me3, synergistically inhibiting FLC expression to promote flowering . As the temperature increases (22 °C), LIKE HETEROCHROMATIN PROTEIN 1 (LHP1), and VRN2 recognize H3K27me3, thus stably maintaining FLC suppression. Multiple genes alter the chromatin structure of the FLC locus to inhibit FLC expression [10, 11].
The entire process of vernalization requires at least 1 month of cold exposure for enabling plants to saturate the vernalization response. During this period, plants must also initiate cold acclimation to endure harsh environments. Cold acclimation is a rapid adaptive response that enables plants to acquire freezing tolerance. In most cases, 1-2 d of exposure to low but non-freezing temperatures is sufficient for plants to acquire cold acclimation . A. thaliana achieves maximum freezing tolerance after 7 d of exposure to temperatures as low as 2 °C [13, 14]. When cold-acclimated plants are treated with non-acclimation-inducing temperatures, they lose their freezing tolerance. This is known as deacclimation . The required duration for deacclimation is much shorter. For instance, 24 h is sufficient to deacclimatize Arabidopsis . There are two basic types of cold-response pathways that enable plants to achieve cold acclimation, namely CBF-dependent and CBF-independent pathways. ICE1-CBF-COR is the core model of the CBF-dependent pathway. Inducer of CBF expression 1 (ICE1) is an MYC-like basic helix–loop–helix transcription factor that can bind to the MYC cis-acting elements in the CBF promoter [16, 17]. HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENE1 (HOS1), SAP AND MIZ1 DOMAIN-CONTAINING LIGASE 1 (SIZ1), and OPEN STOMATA 1 (OST1) regulate ICE1 through ubiquitination, sumoylation, and phosphorylation, respectively, thus participating in the CBF-dependent pathway [18,19,20]. The COR genes comprise four gene families, namely low temperature-inducible (LTI), cold-inducible (KIN), responsive to desiccation (RD), and early dehydration-inducible (ERD) genes . Of note, 10% of all COR genes are regulated by CBF1, CBF2, and CBF3 (or DREB1b, DREB1c, and DREB1a, respectively [22, 23]. At low but non-freezing temperatures, ICE1 directly binds to the CBF promoter to activate its expression, further enhancing COR expression to enhance cold tolerance [24,25,26].
Phytohormones are essential in the regulation of freezing tolerance. JA positively regulates the ICE–CBF pathway to enhance freezing tolerance in Arabidopsis . Blocking JA biosynthesis and signaling produces hypersensitivity to freezing tolerance . The ethylene signaling pathway negatively regulates freezing tolerance. EIN3 suppresses the expression of CBFs and Type-A ARR genes by directly binding to their promoters [24, 28].
The relationship between cold stress response and vernalization in plants remains controversial. In 2004, Sung and Amasino  reported that one distinction between cold acclimation and vernalization is the time lag between the two pathways (approximately 10 d) with Col-FRI seedlings. As such, vernalization occurs approximately 10 d later than cold acclimation, based on the induction time of VIN3 [12, 29]. As the first gene of the vernalization pathway, VIN3 can be detected within 1 day of cold treatment in the background of rapid-flowering ecotypes, whereas in Col-FRI, VIN3 can only be induced by cold exposure longer than 2 weeks [29, 30]. In 2009, Seo et al.  reported that CBFs could induce FLC expression under intermittent cold (0–6 h), thereby delaying flowering. Vernalization could override this effect by inhibiting FLC expression, demonstrating a complicated relationship between these two pathways. HOS1 is a negative regulator of cold responses. The demonstration that HOS1 can regulate FLC expression through chromatin remodeling under cold temperatures providing new insights into the crosstalk between the two pathways . However, Bond et al.  illustrated the independence of cold acclimation and vernalization by showing that none of the key components of cold acclimation signaling, such as ICE1 and HOS1, play a role in VIN3 induction. To this end, in the present study, we conducted transcriptome analysis to explore the relationship between vernalization and cold acclimation.
Transcriptional dynamics of Vernalization in Arabidopsis
With the aim of profiling the whole picture of transcriptional dynamics during vernalization and cold stress response in Arabidopsis, we conducted transcriptome analysis of plants exposed to 42 d of cold required to saturate the vernalization response in Col-FRI . We harvested whole plants of FRI-Col (Col-0 with a functional FRI allele) and set up eight sampling time points, four of which (0 d, 14 d, 28 d, and 42 d) were designed to explore the effects of vernalization and the remaining four (0.5 h, 1 d, 29 d, and 30 d) to explore the effects of cold stress. Of note, the 29 d samples were subjected to 1 d of deacclimation (22 °C) and the 30 d samples were subjected to 1 d of reacclimation (Fig. 1a). Here, we focused on the cold exposure phase of vernalization, and 1 d was sufficient to deacclimate but not to devernalize Arabidopsis. The BGISEQ-500 was used to detect differentially expressed genes (DEGs). The correlation heatmap showed that the T14d, T28d, T42d, T0h, and T29d transcriptomes were similar to one another (Fig. 1b). A total of 31,744 DEGs were identified. At the initiation stage of cold exposure, the expression of only 2709 genes was altered within 0.5 h, and the expression level in these samples was almost one-fourth of that in the 1 d samples (Fig. 1c). In addition, the highest number of DEGs was detected in the 1 d samples, indicating that 1 d was the most drastic response time. Genes in the 30 d samples exhibited no significant changes, even though they experienced 1 d of cold exposure after recovery, which may be explained by the gain of cold acclimation as the plants had already been exposed to cold for 28 d. The number of DEGs was similar in the 14, 28, and 42 d samples, suggesting that plants maintained a high level of response to prolonged cold exposure (Fig. 1c). Among the short-term cold exposure treatments, the T0hVST0.5 h, T0dVST1d, and T29dVST30d samples shared only 982 DEGs. Among long-term cold exposure treatments, the T0hVST14d, T0VS28d, and T0VS42d samples shared 5651 DEGs. These results indicate that short-term cold response is more flexible, whereas long-term cold response is more stable (Fig. 1d).
Time-course analysis was conducted by clustering all genes from different time points to investigate their expression dynamics (Fig. 2). Genes in clusters 1 and 2 showed a rapid response to cold within 0.5 h and 1 d, respectively (Fig. 2a, b). Gene Ontology (GO) analysis indicated that genes involved in response within 0.5 h were sensitive to stress and were enriched in wounding response, defense response, and ethylene-activated signaling pathways (Fig. 2a). Genes involved in response within 1 d were associated with ribosomal assembly and protein translation (Fig. 2b). Genes in clusters 3 and 4 showed similar expression patterns (Fig. 2c, d), exhibiting a rapid response to temperature change at 0.5 h and 29 d, respectively (Fig. 2c, d). GO analysis indicated that the upregulated genes in these two clusters were enriched in DNA repair and RNA modification, perhaps because such a lasting response can induce DNA damage and enhance replication (Fig. 2c, d). Genes in clusters 5 and 6 were downregulated during cold exposure and were sensitive to increases in temperature (Fig. 2e, f). GO analysis revealed that genes in cluster 5 were enriched in the reductive pentose phosphate cycle, redox processes, and cytokinin response, while those in cluster 6 were enriched in cell division and cell cycle. These results indicate that the cell cycle, energy consumption, and oxidative activity are effective at relatively high temperatures, but are suppressed at low temperatures (Fig. 2e, f). Genes in cluster 7 maintained high expression after 1 d of cold exposure (Fig. 2g) and were primarily enriched in phosphorylation-related processes, including the MKK and CIPK9 signaling pathways as well as intracellular protein transport (Fig. 2g). This result suggests that phosphorylation plays an important role in the response to long-term temperature changes. Genes in cluster 8 were enriched in cold response and cold acclimation (Fig. 2h); the expression of genes involved in cold response peaked at 14 d and stably dropped thereafter, whereas that of genes involved in cold acclimation peaked at 42 d. Genes in cluster 9 were involved in the regulation of flower development (Fig. 2i). In addition, genes in cluster 8 were sensitive to temperature change, but those in cluster 9 showed no such activity. Collectively, these expression patterns indicate the primary relationships and differences between short- and long-term responses to cold (Fig. 2h, i) (Table S1).
Relationship between cold acclimation and Vernalization
Cold acclimation and vernalization were likely initiated with a time lag of approximately 10 d. To explore the overlaps and interactions between cold acclimation and vernalization, we focused on genes involved in the relevant pathways (Table S2). The expression of CBF1, which is the key gene involved in cold stress response, peaked within 0.5 h and remained stable thereafter (Fig. 3a). The expression of other genes known to be involved in the cold response also showed a typical upward trend. ICE1, COR15A, COR15B, COR47, COR413PM1, COR413IM1, LTI30, LTI65, LTI78, ERD2, ERD3, ERD4, ERD7, ERD10, ERD14, KIN1and KIN2 were significantly upregulated within 1 d (Fig. 3f) (Table S2). Conversely, the expression of FLC, which is the key gene involved in vernalization, was initially suppressed at 14 d, along with the induction of VIN3, which is considered the first gene activated in the vernalization pathway (Fig. 3a). In addition, NTL8, which was recently shown to upregulate VIN3 under long-term cold , showed similar expression to VIN3. The expression of the PRC2 genes showed a typical upward trend. VAL1 expression showed an obvious upward trend during vernalization, while VAL2 expression dropped to the normal level after a slight increase (Table S2). Notably, the expression of almost all genes in the vernalization pathway was altered to regulate FLC expression after 14 d (Fig. 3e). Box plots showed that the expression patterns of genes involved in the two pathways were distinct in terms of the time point of their change (1 d for cold acclimation and 14 d for vernalization) (Fig. 3g, h). We confirmed this observation using qPCR, and the patterns of FLC, CBF1, and VIN3 expression were consistent between qPCR and RNA-seq (Fig. 3b, c, d).
HOS1 is a crosstalk gene between cold stress and vernalization, which regulates FLC expression under intermittent cold conditions and physically interacts with ICE1 to mediate its ubiquitination [18, 32]. The constant increase in FLC expression and decrease in ICE1 expression within 1 d may be attributed to this function of HOS1. However, the upward trend of FLC expression and the downward trend of ICE1 expression did not last long. After 1 d of cold exposure, the expression of these two genes showed converse trends, indicating that such interactive activity would not affect the timing of initiation of the two pathways. Therefore, consistent with the time-lag notion, rapid response genes are not involved in quantitative response, that is, cold acclimation and vernalization are independent from the perspective of overlapping regulatory genes.
Network analysis of cold stress and Vernalization based on weighted gene co-expression network analysis (WGCNA)
Furthermore, WGCNA was performed to explore the relationship between cold stress and vernalization. The analysis yielded 28 modules (Fig. 4a and b; Table S3). Three core networks were identified as the key models of cold response, and genes with known functions were selected. We designated these models as the ERF, LTI, and SCC networks based on hub genes with high module membership (MM) values (Table S4). Genes with MM values over 0.98, are listed in Table 1. Genes in the ERF network belonged to cluster 1 and were involved in an instant cold response. Consistent with the GO_P analysis of cluster 1, the ERF network included ethylene and jasmonic acid (JA) response factors, with the hub genes ERF11, ERF104, JAZ7, and WRKY40 (MM > 0.99) (Fig. 4c). ETHYLENE RESPONSIVE ELEMENT-BINDING FACTORs (ERFs) are the core ethylene response factors, and five ERFs, namely ERF1, ERF2, ERF4, ERF5, and ERF6, were also involved in this network. JASMONATE-ZIM-DOMAIN PROTEINs (JAZs) play pivotal roles in JA signaling, and five JAZs, namely JAZ1, JAZ5, JAZ7, JAZ8, and JAZ10, were involved in this network (Fig. 4c). Additionally, the ethylene biosynthesis genes ACS6 and ACS11 and the JA biosynthesis gene AOC3 functioned at the relative outer circle of the ERF network, suggesting that ethylene and JA are the key hormones for short-term cold response. ZAT6, ZAT7, and ZAT12 of the C2H2 ZINC FINGER TRANSCRIPTION FACTOR (ZFP) family were also part of the ERF network, and these genes were primarily involved in the stress response (Fig. 4c).
The LTI network involved typical cold acclimation genes from cluster 8. The expression levels of genes in T14d, T28d, T30d, and T42d samples were significantly different from those at T0h (Fig. 4g) (Table S4). The hub genes of this network included LOW TEMPERATURE INDUCED 30 (LTI30), LATE EMBRYOGENESIS ABUNDANT 14 (LEA14), and A. thaliana ALPHA CARBONIC ANHYDRASE 8 (ACA8) (MM > 0.98) (Fig. 4d). Notably, IAA1 and IAA29 were also part of this network, indicating that auxin is involved in cold acclimation.
The SCC network was the largest model detected in the present study. Genes in this network were first downregulated at T0.5h and T1d, and then upregulated at T14d, T28d, and T42d, indicating that this network is likely involved in response to long-term cold (Fig. 4h; Table S4). SCC2, SCC3, and CUL4 (MM > 0.99) were the hub genes of the SCC network (Fig. 4e). CUL4 is one of the CULLIN (CUL) RING UBIQUITIN LIGASEs (CRLs), which are involved in substrate ubiquitination . SCC2 and SCC3 are essential for maintaining centromere cohesion during anaphase I. Many genes from the SCC network are involved in chromatin modification. These included HOS1, the trimethylase ATX1, LD, and FLD from the autonomous flowering pathway, SNL1 related to deacetylation, and others [32, 38,39,40,41] (Fig. 4e). These three networks represent three major parts of the cold response. The first is an instantaneous response, particularly at T0.5h, which is related to JA/ethylene signaling. The second is a rapid and lasting response that is sensitive to 28-30d temperature change, which is related to cold acclimation pathway. The third is a stable response featuring a consistent expression level at T14d, T28d, and T42d, which is related to chromatin modification (Fig. 4f, g, h). Notably, the VIN3 and FLC as core genes of vernalization were not involved in any module of WGCNA, indicating that the cold response is independent of vernalization.
Alternative splicing mediation during Vernalization
Alternative splicing is a ubiquitous co-transcriptional RNA modification through which multiple transcripts can be generated from a single gene. Temperature is closely associated with alternative splicing. Several mechanisms of alternative splicing have been reported, including skipped exons (SE) (a specific exon is excluded from mature mRNA), mutually exclusive exons (MXEs) (choice between two constitutive exons), alternative 3′/5 splicing sites (A3SS/A5SS) (distinct 3′ or 5′ splicing sites are generated in the resulting isoforms), and retained introns (RIs) . RI is the predominant form of alternative splicing in plants and generates transcripts with premature termination codons (PTCs), thus leading to nonsense mRNA decay (NMD) . A total of 1540 differential alternative splicing (DAS) genes were identified, accounting for 4.85% of all DEGs. Overall, the proportion of RIs decreased during cold exposure, while MXEs appeared after 1 d. The proportion of A3SS also significantly increased (Fig. 5a). These results suggest that plants attempt to alter splicing patterns to cope with environmental cues more efficiently.
MAF1 (FLM) and MAF2 are regulated by alternative splicing. They are homologous to FLC and serve as floral repressors during vernalization. However, their expression was upregulated during cold exposure and downregulated following recovery at 22 °C, as expected (Figure S1). Therefore, MAF1 and MAF2 are involved in an alternative splicing-mediated thermosensory pathway rather than vernalization. FLM regulates flowering through change in the FLM-β/δ ratio when plants experience temperature fluctuations between 16 °C and 27 °C . At 4 °C, FLM-β was the only upregulated transcript, while FLM-δ expression was relatively stable (Fig. 5b), indicating that FLM suppresses flowering in an FLM-β-dependent and FLM-δ-independent manner under cold but non-freezing conditions. MAF2 also responds to cold by altering the expression of the sole transcript MAF2var1 .
Regarding specific response pathways, the alternative splicing mechanism significantly affected every known step during cold signaling (Table S5). The cold response starts with signal transduction. The Ca2+-permeable cyclic nucleotide-gated channel 5 (CNGC5), CNGC6, and phytochromes PHYA and PHYB are the most upstream genes that respond to temperature signals in an alternative splicing-dependent manner (Fig. 5b) [46, 47]. CBL-interacting protein kinase 3 (CIPK3) is a crucial kinase for Ca2+ signal transduction, and PHYTOCHROME-INTERACTING FACTOR 4 (PIF4) and PIF7 are important factors for thermosensory regulation via phyB; all these factors are regulated through alternative splicing [48,49,50]. The overall expression levels of CIPK3, PIF4, and PIF7 did not change significantly, but their functional transcripts were differentially expressed (Fig. 5b). BES1, EIN4, and SPY are essential regulators of brassinosteroid, ethylene, and GA signaling; ICE1, COR15A, and COR27 are the core genes of the cold acclimation pathway; and ATX1, ATX2, and FCA are the chromatin or RNA regulators. All these genes are regulated through alternative splicing (Fig. 5b) [51,52,53,54,55,56]. Finally, we selected two vernalization genes, namely VRN5 and VAL2, whose overall expression was not upregulated, but the proportion of their sole transcripts was altered (Fig. 5b, S1).
We classified these analyzed genes into three types: change in the proportion of different transcripts, induction of functional transcripts, and synergistic function of all transcripts (Fig. 5b). CNGC5, ICE1, EIN4, PHYB, PIF4, RAB8, ATX1, VRN5, and VAL2 belonged to type one, in which the proportions of different transcripts were altered. This type of gene may have antagonistic functions. FLM, MAF2, CIPK3, COR27, and ATX2 belonged to type two, whose transcript expression was distinctly increased; CNGC6, COR15A, BES1, PHYA, PIF7, and FCA belonged to type three, which functioned synergistically. Therefore, alternative splicing is involved in the intricate responses to temperature changes.
Relationship between cold acclimation and Vernalization
Temperature is an essential environmental stimulus that significantly affects plant growth and development. Plants use various strategies to cope with different durations of cold. Vernalization is a quantitative cold-sensing process that enables plants to flower in warm springs. Simultaneously, plants must also initiate a stress response to cold during vernalization. Cold acclimation enables plants to rapidly adapt to cold temperatures. A time lag of approximately 10 d between the initiation of the two pathways has been proposed , and whether these two successive pathways interact under the same cold treatment should be elucidated. The expression of VIN3 cannot be detected before T14 d. Additionally, a recent study showed that NTL8, an upstream regulator of VIN3, does not function similarly to other genes involved in cold response, such as COR47 , whose expression pattern was consistent with that of VIN3 (Fig. 3c). FLC was upregulated by cold at T1d. This effect was overridden by vernalization following VIN3 induction, as FLC expression continued to decrease. Core genes of cold stress response CBF1 and CBF2 were upregulated at T0.5h (Fig. 3). Thus, the aforementioned time-lag notion was confirmed by the CBFs and VIN3 induction times (Fig. 3). In addition, other genes involved in the cold signaling pathway, such as CAMTA3, MPK3, MPK6, CIPK7, and OST1, also show increased expression to mediate cold tolerance [57, 58] (Table S5). The independence between cold acclimation and vernalization has been illustrated by the demonstration that core components of the cold-responsive pathway, such as ICE1, HOS1, and IP3-related CVP2, have no effect on VIN3 expression . The WGCNA results provided complementary verification of the detection of CBFs, ICE1, HOS1, and CVP2 LIKE 1(CVL1) in modules, while VIN3 and FLC were co-expressed with no genes related to the cold-response pathway.
Sensors of cold acclimation and Vernalization
Several recent studies have focused on temperature sensors. COLD1 acts as a cold sensor at 2–4 °C and confers chilling tolerance in Japonica rice . COLD1 is homologous to AtGTG1 and AtGTG2 and functions together with RGA1 to mediate the cold-induced influx of extracellular Ca2+ in rice . The expression levels of GTG1 and GTG2 decreased within 1 d and returned to normal levels after 14 d (Figure S2). CNGCs are nonspecific cation channels, putatively downstream of GTG1, and the expression levels of CNGCs and GTGs were similar (Figure S2) . PHYB acts as a thermal sensor at 12–27 °C, sensing temperature change through thermal reversion . At the transcriptional level, PHYB senses low temperature via downregulation for 1 d, followed by upregulation under long-term cold conditions (Figure S2). PIFs, with known binding sites, are downstream of PHYB, and their expression pattern is similar to that of PHYB (Figure S2). According to Jung et al. , ELF3 responds to increasing temperatures within 22–27 °C by reversibly forming liquid droplets. Interestingly, the expression pattern of ELF3 resembled that of GTGs or PHYB—decreased within 24 h of cold treatment and increased thereafter, and the three reported thermal sensors showed similar expression changes when treated at 4 °C (Figure S2). Therefore, low temperatures may have a similar effect on these thermal sensors at the transcriptional level.
Whether sensors of cold stress and vernalization are the same remains debatable. Suppression of FLC expression involves VIN3-dependent and independent pathways . Temperature-dependent growth has been reported as a long-term thermosensor for VIN3-dependent repression of FLC . Reduced NTL8 protein accumulation due to slow growth at low temperatures is the major cause of VIN3 accumulation . At the transcriptional level, the expression pattern of NTL8 resembled that of VIN3, being upregulated at 14 d and downregulated at 29 d (Fig. 4f). The effect of cold treatment on NTL8 mRNA induction in Col-FRI background was more obvious than that in the Col background . Physical disruption of the FLC-containing gene loop is considered as a sensor of the VIN3-independent pathway of FLC repression during the vernalization process [64, 65]. The transcription activity of FLC reportedly drops to its lowest level within 1 week. However, FLC showed a short-term HOS1-dependent upward trend in response to cold at 1 d (Fig. 3a). This trend can be maintained at 10 d of cold exposure, as recently demonstrated by another transcriptomic analysis . One of the supporting pieces of evidence is that UPTREAM OF FLC (UFC), which is located 4.7 kb upstream of FLC, is suppressed in parallel with FLC during vernalization . Indeed, the expression of UFC was also upregulated at 1 d and then downregulated thereafter.
Alternative splicing mediation during cold exposure
Temperature is considered a key factor in the regulation of alternative splicing. Temperature-dependent alternative splicing is associated with various temperature-related processes, such as hibernation in mammals and cold acclimation in plants and fish [67,68,69]. Calixto et al.  found that plants undergo rapid and dynamic alternative splicing in response to short-term cold conditions. With the aim of exploring alternative splicing mediation under long-term cold conditions, we assessed the overall trends of change and regulation of individual genes during 42 d of cold treatment. Alternative splicing was active throughout cold exposure. Many splicing regulation-related genes showed increased expression during cold treatment (Figure S3). Overall, the proportion of the five mechanisms was altered during cold exposure to generate additional functional transcripts, thereby enhancing cold tolerance more efficiently (Fig. 5a). Specific genes whose overall expression did not change may be greatly upregulated as a single transcript. Typically, genes involved in cold acclimation and vernalization respond to cold via three major mechanisms: change in the proportion of different transcripts, induction of functional transcripts, and synergistic function of all transcripts, which enable plants to respond to environmental cues more intricately and efficiently (Fig. 5, S1).
In conclusion, using transcriptomic analysis of the entire vernalization process, we uncovered the independence of cold acclimation and vernalization and further revealed the response networks involved in prolonged cold exposure of plants. Plants tend to alter splicing patterns (decreasing the proportion of RIs) in response to long-term cold, and alternative splicing mechanisms may regulate the entire process of cold response. Thermal sensors exhibit similar expression patterns under non-freezing but cold conditions, and they may sense non-freezing cold at the transcriptional level.
Materials and methods
Plant material and growth conditions
Arabidopsis FRI-Col (Col-0 with a functional FRI allele) was used in this experiment. A functional FRI locus from the Sf2 line was introgressed into Col to construct FRI-Col by the R. Amasino lab , and the FRI-Col seeds were provided by Dr. Yuehui He (Shanghai Center for Plant Stress Biology, Chinese Academy of Sciences). Permissions for using these materials were obtained from the Chinese Academy of Sciences. Seeds were surface-sterilized with 75% ethyl alcohol for 1 min, followed by 10% sodium hypochlorite for 15 min, washed six times with sterile water, and stratified at 4 °C for 2 d before being sown on 1/2MS medium. Seedlings were incubated in growth chambers at 22 °C under a 16-8 h day-night period for approximately 2 wk. and then three biological replicates of samples were harvested when they bore four true leaves (T0h). The remaining samples were transferred to an 8-16 h day-night period at 4 °C for 30 min (T0.5h), 1 d (T1d), 14 d (T14d), or 28 d (T28d), and finally harvested. The 29 d samples were subjected to 1 d of recovery at 22 °C (T29d), and the 30 d samples were exposed to cold for 1 d following recovery (T30d). The 42 d samples were subjected to an additional 14 d of treatment at 4 °C under short-day conditions and then harvested (T42d). All the samples were harvested at zeitgeber time 7 (ZT 0 is light-on) except that T0.5h was at ZT 7.5.
RNA extraction and RNA-seq library construction
Total RNA was extracted using RNAiso Plus (TaKaRa) and processed for mRNA enrichment and rRNA removal. The mRNA was enriched with a polyA tail using magnetic beads with OligodT. Fragmentation buffer was added to break the mRNA into short fragments at high temperature, and the fragmented mRNA was used as the template to synthesize the first-strand cDNA. Next, the second-strand cDNA was synthesized, and the recovered cDNA was purified using a commercial kit. Next, the base “A” was added to the 3 ends of the cDNA, and a linker was connected. The size of the fragment was determined, and the fragments were subjected to PCR amplification. The quality of the constructed library was checked, and the libraries were sequenced.
DEGs and DAS genes
High-throughput sequencing was performed using the BGISEQ-500 platform. After several data processing steps (including removal of adaptor sequences, null reads, and low-quality reads), pure reads were obtained from the original sequence. After obtaining clean reads, HISAT was used to align the clean reads to the reference genome (GCF_000001735.4_TAIR10.1). Bowtie2  was used to align clean reads to the reference gene sequences, and RSEM  was used to calculate the gene expression level of each sample.
rMATS was used to detect DAS genes between different samples and splicing events of the samples. rMATS was used for the analysis of DAS genes based on RNA-seq data. It uses the rMATS statistical model to quantify the expression of alternative splicing events in different samples and then calculates the P value with a likelihood-ratio test to indicate whether the two groups of samples are in IncLevel (Inclusion Level), which uses the Benjamini–Hochberg algorithm for correcting the P values to obtain the false discovery rate.
To assess similarities in expression patterns among the groups, we analyzed the transcriptome profiles of biological replicates. The log2-normalized FPKM values of gene expression were input into the WGCNA package in R  to generate gene networks. A standard process was used to minimize the noise. An adjacency matrix was constructed using a soft threshold power of 9. Networks were identified using a dynamic tree-cut algorithm with a minimum cluster size of 30 and merging threshold of 0.25. Hub genes were identified based on their eigengene connectivity (KME) . Networks were visualized using Cytoscape v3.5.1 (https://cytoscape.org/).
qRT-PCR of DEGs
The reliability of DEGs or transcripts identified through RNA-seq was evaluated through qRT-PCR of FLC, VIN3, CBF1, CBF2, MAF1, and SOC1. An Eppendorf Mastercycler Ep RealPlex 2S (Hamburg, Germany) fluorescence quantifier was used. Reactions were performed using the 2× SYBR Green qPCR Master Mix (Bimake), following the manufacturer’s instructions. The specific reaction system contained 10.0 μL of SYBR® Premix Ex Taq™ II, 1.0 μL each of 10 μM forward and reverse primers, 5.0 μL of cDNA template, and 3.0 μL of ddH2O. The reaction conditions were as follows: 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s, 55 °C for 15 s, and 72 °C for 20 s. Finally, gene expression levels were quantified using the 2-ΔΔCT method.
Gazzani S, Gendall AR, Lister C, Dean C. Analysis of the molecular basis of flowering time variation in Arabidopsis accessions. Plant Physiol. 2003;132(2):1107–14 https://doi.org/10.1104/pp.103.021212.
Michaels SD, He Y, Scortecci KC, Amasino RM. Attenuation of FLOWERING LOCUS C activity as a mechanism for the evolution of summer-annual flowering behavior in Arabidopsis. Proc Natl Acad Sci U S A. 2003;100(17):10102–7 https://doi.org/10.1073/pnas.1531467100.
Searle I, He Y, Turck F, et al. The transcription factor FLC confers a flowering response to vernalization by repressing meristem competence and systemic signaling in Arabidopsis. Genes Dev. 2006;20(7):898–912 https://doi.org/10.1101/gad.373506.
Choi K, Kim J, Hwang HJ, et al. The FRIGIDA complex activates transcription of FLC, a strong flowering repressor in Arabidopsis, by recruiting chromatin modification factors. Plant Cell. 2011;23(1):289–303 https://doi.org/10.1105/tpc.110.075911.
Wood CC, Robertson M, Tanner G, Peacock WJ, Dennis ES, Helliwell CA. The Arabidopsis thaliana vernalization response requires a polycomb-like protein complex that also includes VERNALIZATION INSENSITIVE 3. Proc Natl Acad Sci U S A. 2006;103(39):14631–6 https://doi.org/10.1073/pnas.0606385103.
Mylne JS, Barrett L, Tessadori F, et al. LHP1, the Arabidopsis homologue of HETEROCHROMATIN PROTEIN1, is required for epigenetic silencing of FLC. Proc Natl Acad Sci U S A. 2006;103(13):5012–7 https://doi.org/10.1073/pnas.0507427103.
Oono Y, Seki M, Satou M, et al. Monitoring expression profiles of Arabidopsis genes during cold acclimation and deacclimation using DNA microarrays. Funct Integr Genomics. 2006;6(3):212–34 https://doi.org/10.1007/s10142-005-0014-z.
Kim YS, Lee M, Lee JH, Lee HJ, Park CM. The unified ICE-CBF pathway provides a transcriptional feedback control of freezing tolerance during cold acclimation in Arabidopsis. Plant Mol Biol. 2015;89(1-2):187–201 https://doi.org/10.1007/s11103-015-0365-3.
Dong CH, Agarwal M, Zhang Y, Xie Q, Zhu JK. The negative regulator of plant cold responses, HOS1, is a RING E3 ligase that mediates the ubiquitination and degradation of ICE1. Proc Natl Acad Sci U S A. 2006;103(21):8281–6 https://doi.org/10.1073/pnas.0602874103.
Miura K, Jin JB, Lee J, et al. SIZ1-mediated sumoylation of ICE1 controls CBF3/DREB1A expression and freezing tolerance in Arabidopsis. Plant Cell. 2007;19(4):1403–14 https://doi.org/10.1105/tpc.106.048397.
Mustilli AC, Merlot S, Vavasseur A, Fenzi F, Giraudat J. Arabidopsis OST1 protein kinase mediates the regulation of stomatal aperture by abscisic acid and acts upstream of reactive oxygen species production. Plant Cell. 2002;14(12):3089–99 https://doi.org/10.1105/tpc.007906.
Kidokoro S, Yoneda K, Takasaki H, Takahashi F, Shinozaki K, Yamaguchi-Shinozaki K. Different cold-signaling pathways function in the responses to rapid and gradual decreases in temperature. Plant Cell. 2017;29(4):760–74 https://doi.org/10.1105/tpc.16.00669.
Shi Y, Tian S, Hou L, et al. Ethylene signaling negatively regulates freezing tolerance by repressing expression of CBF and type-a ARR genes in Arabidopsis. Plant Cell. 2012;24(6):2578–95 https://doi.org/10.1105/tpc.112.098640.
Jia Y, Ding Y, Shi Y, Zhang X, Gong Z, Yang S. The cbfs triple mutants reveal the essential functions of CBFs in cold acclimation and allow the definition of CBF regulons in Arabidopsis. New Phytol. 2016;212(2):345–53 https://doi.org/10.1111/nph.14088.
Hu Y, Jiang L, Wang F, Yu D. Jasmonate regulates the inducer of cbf expression-C-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. Plant Cell. 2013;25(8):2907–24 https://doi.org/10.1105/tpc.113.112631.
Hu Y, Jiang Y, Han X, Wang H, Pan J, Yu D. Jasmonate regulates leaf senescence and tolerance to cold stress: crosstalk with other phytohormones. J Exp Bot. 2017;68(6):1361–9 https://doi.org/10.1093/jxb/erx004.
Bond DM, Wilson IW, Dennis ES, Pogson BJ, Jean FE. VERNALIZATION INSENSITIVE 3 (VIN3) is required for the response of Arabidopsis thaliana seedlings exposed to low oxygen conditions. Plant J. 2009;59(4):576–87 https://doi.org/10.1111/j.1365-313X.2009.03891.x.
Seo E, Lee H, Jeon J, et al. Crosstalk between cold response and flowering in Arabidopsis is mediated through the flowering-time gene SOC1 and its upstream negative regulator FLC. Plant Cell. 2009;21(10):3185–97 https://doi.org/10.1105/tpc.108.063883.
Méndez-Vigo B, Picó FX, Ramiro M, Martínez-Zapater JM, Alonso-Blanco C. Altitudinal and climatic adaptation is mediated by flowering traits and FRI, FLC, and PHYC genes in Arabidopsis. Plant Physiol. 2011;157(4):1942–55 https://doi.org/10.1104/pp.111.183426.
Zhao Y, Antoniou-Kourounioti RL, Calder G, Dean C, Howard M. Temperature-dependent growth contributes to long-term cold sensing [published correction appears in nature. 2020 Sep;585(7824):E8]. Nature. 2020;583(7818):825–9 https://doi.org/10.1038/s41586-020-2485-4.
Kim S, Choi K, Park C, Hwang HJ, Lee I. SUPPRESSOR OF FRIGIDA4, encoding a C2H2-type zinc finger protein, represses flowering by transcriptional activation of Arabidopsis FLOWERING LOCUS C. Plant Cell. 2006;18(11):2985–98 https://doi.org/10.1105/tpc.106.045179.
Wang Z, Cao H, Sun Y, et al. Arabidopsis paired amphipathic helix proteins SNL1 and SNL2 redundantly regulate primary seed dormancy via abscisic acid-ethylene antagonism mediated by histone deacetylation. Plant Cell. 2013;25(1):149–66 https://doi.org/10.1105/tpc.112.108191.
Monteuuis G, Wong JJL, Bailey CG, Schmitz U, Rasko JEJ. The changing paradigm of intron retention: regulation, ramifications and recipes. Nucleic Acids Res. 2019;47(22):11497–513 https://doi.org/10.1093/nar/gkz1068.
Airoldi CA, McKay M, Davies B. MAF2 is regulated by temperature-dependent splicing and represses flowering at low temperatures in parallel with FLM. PLoS One. 2015;10(5):e0126516 Published 2015 May 8. https://doi.org/10.1371/journal.pone.0126516.
Fischer C, DeFalco TA, Karia P, et al. Calmodulin as a Ca2+−sensing subunit of Arabidopsis cyclic nucleotide-Gated Channel complexes. Plant Cell Physiol. 2017;58(7):1208–21 https://doi.org/10.1093/pcp/pcx052.
Kim KN, Cheong YH, Grant JJ, Pandey GK, Luan S. CIPK3, a calcium sensor-associated protein kinase that regulates abscisic acid and cold signal transduction in Arabidopsis. Plant Cell. 2003;15(2):411–23 https://doi.org/10.1105/tpc.006858.
Kumar SV, Lucyshyn D, Jaeger KE, et al. Transcription factor PIF4 controls the thermosensory activation of flowering. Nature. 2012;484(7393):242–5 Published 2012 Mar 21. https://doi.org/10.1038/nature10928.
Galvāo VC, Fiorucci AS, Trevisan M, et al. PIF transcription factors link a neighbor threat cue to accelerated reproduction in Arabidopsis. Nat Commun. 2019;10(1):4005 Published 2019 Sep 5. https://doi.org/10.1038/s41467-019-11882-7.
Antoniou-Kourounioti RL, Hepworth J, Heckmann A, et al. Temperature sensing is distributed throughout the regulatory network that controls FLC Epigenetic silencing in vernalization. Cell Syst. 2018;7(6):643–655.e9 https://doi.org/10.1016/j.cels.2018.10.011.
Horii Y, Shiina T, Uehara S, et al. Hypothermia induces changes in the alternative splicing pattern of cold-inducible RNA-binding protein transcripts in a non-hibernator, the mouse. Biomed Res. 2019;40(4):153–61 https://doi.org/10.2220/biomedres.40.153.
Lee I, Amasino RM. Effect of Vernalization, photoperiod, and light quality on the flowering phenotype of Arabidopsis plants containing the FRIGIDA gene. Plant Physiol. 1995;108(1):157–62 https://doi.org/10.1104/pp.108.1.157.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323 Published 2011 Aug 4. https://doi.org/10.1186/1471-2105-12-323.
We thank Dr. Yuehui He (National Key Laboratory of Plant Molecular Genetics & Shanghai Center for Plant Stress Biology, The Chinese Academy of Sciences) for providing the seeds of Col-FRI.
This work was supported by grants from the National Natural Science Foundation of China (31930100, 31872146), and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
Authors and Affiliations
State Key Laboratory of Crop Genetics and Germplasm Enhancement, Key Laboratory of Landscaping, Ministry of Agriculture and Rural Affairs, Key Laboratory of Biology of Ornamental Plants in East China, National Forestry and Grassland Administration, College of Horticulture, Nanjing Agricultural University, Nanjing, 210095, China
All authors have read and approved the manuscript. FL, QH and JJ designed the experiments. QH and JJ prepared samples for RNA-seq analysis. FL analyzed the data. FL and JJ wrote the manuscript. FC and JJ revised the manuscript.
The authors declare that Professor R. Amasino has given the permission for Dr. Yuehui He to provide seeds for this experiment, and permissions for using these materials was obtained from the Chinese Academy of Sciences.
Consent for publication
The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Heatmap showing the expression pattern of selected DEGs. Gene expression data were normalized to Log2(FPKM+1); red and blue represent up- and downregulated genes, respectively.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.