Rapid and dynamic alternative splicing impacts the Arabidopsis cold response transcriptome

Background Plants have adapted to tolerate and survive constantly changing environmental conditions by re-programming gene expression. The scale of the contribution of alternative splicing (AS) to stress responses has been underestimated due to limitations in RNA-seq analysis programs and poor representation of AS transcripts in plant databases. Significantly, the dynamics of the AS response have not been investigated but this is now possible with accurate transcript quantification programs and AtRTD2, a new, comprehensive transcriptome for Arabidopsis. Results Using ultra-deep RNA-sequencing of a time-course of Arabidopsis thaliana plants exposed to cold treatment, we identified 8,949 genes with altered expression of which 2,442 showed significant differential alternative splicing (DAS) and 1,647 genes were regulated only at the level of AS (DAS-only). The high temporal resolution demonstrated the rapid induction of both transcription and AS resulting in coincident waves of differential expression (transcription) and differential alternative splicing in the first 6-9 hours of cold. The differentially expressed and DAS gene sets were largely non-overlapping, each comprising thousands of genes. The dynamic analysis of AS identified genes with rapid and sensitive AS within 3 h of transfer to the cold (early AS genes), which were enriched for splicing and transcription factors. A detailed investigation of the novel cold-response DAS-only gene, U2B”-LIKE, suggested that it regulates AS and is required for tolerance to freezing. Conclusions Our data indicate that transcription and AS are the major regulators of transcriptome reprogramming that together govern the physiological and survival responses of plants to low temperature.


50
Background 51 Plants adapt to and survive adverse environmental conditions by reprogramming their 52 transcriptome and proteome. Low temperatures negatively affect growth and development 53 but, in general, plants from temperate climatic regions can tolerate chilling temperatures (0-54 15°C), and can increase their freezing tolerance by prior exposure to low, non-freezing 55 temperatures through a process called cold acclimation. Cold acclimation involves multiple 56 physiological, biochemical and molecular changes to reduce growth, modify membrane 57 properties, and produce the proteins and metabolites required to protect cell integrity 58 throughout the period of freezing exposure [1][2][3]. In Arabidopsis thaliana, these cold- genes [1,3,6,7]. Increased expression of these genes during cold acclimation increases 68 freezing tolerance, which can also be achieved by ectopic expression of the CBFs in the 69 absence of cold acclimation [1][2][3]6]. Besides the CBF regulon, other transcriptional 70 pathways are required for the activation of cascades of gene expression which together 71 drive the cold response, acclimation and plant survival [8]. Plants also adapt to daily and 72 seasonal changes in temperature, and the circadian clock provides a mechanism by which 73 they anticipate the predictable diel cycles of temperature and light/dark to optimise gene 74 expression and consequent physiology [9]. The circadian clock also regulates cold-75 biological repeats and time-of-day variation in expression to be taken into account in the 157 statistical analysis to produce more accurate and robust results. It is important to note that in 158 the time-series analysis, we compared gene and transcript abundances between equivalent 159 time-points at 20°C and those in day 1 and day 4 at 4°C to remove the effects of time-of-day 160 so that the changes detected are due to reduction in temperature (contrast groups in 161 Additional file 1: Figures S1 and S2). We firstly analysed differential expression at the gene 162 level (DE). A gene was considered differentially expressed if it showed a log2-fold change ≥1 163 (≥2-fold change) in expression in at least two consecutive contrast groups (adjusted p<0.01) 164 ( Fig. 2a; Additional file 1: Figure S1; Additional file 2: Table S1). Using these stringent 165 criteria, we identified a total of 7,302 genes which were significantly differentially expressed 166 in response to low temperature when compared to 20°C levels. Of these, 48.2% were up-167 regulated and 51.8% down-regulated (Additional file 1: Figure S3; Additional file 3: Table  168 S2). 169 Secondly, we used the transcript-level data to identify genes that were differentially 170 alternatively spliced (DAS) between contrast groups. F-tests were carried out to examine the 171 consistency of expression changes among the transcripts and gene (see Methods) to detect 172 DAS genes. Criteria for genes being DAS were that in at least two consecutive contrast 173 groups 1) at least one of the transcripts differed significantly from the gene with an adjusted 174 p-value <0.01, and 2) at least one of the transcripts of the gene showed a Δ Percent Spliced 175 (ΔPS) of ≥0.1 (to keep genes where there is a significant change coming from transcript(s) 176 with large differences in their relative abundance) ( Fig. 2a; Additional file 1: Figure S2; 177 Additional file 2: Table S1). We identified 2,442 DAS genes (Fig. 2a) of which 795 were also 178 DE genes (regulated by both transcription and AS) and 1,647 genes that were not DE 179 (regulated by AS only). Thus, of the total of 8,949 genes, which exhibited significantly altered 180 levels of differential gene and/or transcript level expression, 27.3% were differentially 181 alternatively spliced, consistent with widespread DAS in response to cold ( Fig. 2a; Additional 182 file 1: Figure S2; Additional file 2: Table S1). At any particular time-point, between ca. 600-183 3,700 genes were DE and ca. 400-1,450 genes were DAS when compared to 20°C levels 184 (Additional file 1: Figure S3). 185 To pinpoint the individual transcripts responsible for a gene being identified as a DAS 186 gene, a differential transcript usage (DTU) analysis was carried out by testing the expression 187 changes of every transcript against the weighted average of all the other transcripts of the 188 gene using a t-test (Additional file 1: Figure S2). DTU transcripts were identified as those 189 which differed from the gene level and which showed a ΔPS of ≥0.1 in at least two 190 consecutive contrast groups with an adjusted p-value of <0.01. In total, 34% (4,105) of 191 expressed transcripts of DAS genes were classed as DTU (Additional file 2: Table S1) of 192 which ~70% were protein-coding and ~30% contained premature termination codons (PTC). 193  We next explored differences between immediate responses upon transfer to low 224 temperature and the response to prolonged cold acclimation by comparing which genes 225 were DE and DAS in day 1 and/or day 4 after transfer to 4 o C. Around 50% (3,573 genes) 226 and 60% (1,440) of DE and DAS genes, respectively, were common to both days with the 227 remainder being either unique to day 1 or day 4 (Additional file 1: Figure S4; Additional file 2: 228 Table S1). Thus, changes in gene-level expression and AS occurred throughout the cold 229 period: either transiently (occurring in day 1 at 4°C and returning to 20°C levels by day 4), 230 persisting throughout the cold period (occurring in day 1 at 4°C and remaining at day 4), or 231 occurring later, only in day 4. We propose that these patterns of gene-level expression and 232 AS reflect different contributions to low temperature perception, initial cold responses and 233 physiological acclimation to cold and freezing temperatures. 234

DE and DAS analyses identify novel cold response genes and AS of cold regulators 236
Previous analyses of differential gene expression in wild-type Arabidopsis plants exposed to 237 cold used microarrays [12,33,34] and, more recently, 7,35] (Additional files 4 238 and 5: Table S3 and S4, respectively). There was a substantial overlap between the cold 239 response DE genes identified in those studies and the DE genes identified here (Fig. 2b). 240 Critically, we identified an additional 1,708 novel DE cold response genes (Additional file 6: 241 Table S5). As expected, we showed cold induction of CBFs and selected COR genes (none 242 of which undergo AS) (Additional file 1: Figure S5a and b). However, the first significant 243 increase in CBF expression (CBF2) was detected between 3 and 6 h after onset of cold 244 treatment (applied at dusk) whereas expression of CBFs has been detected within less than 245 1 h and peaked at around 1-2 h in other studies conducted in constant light [36]. The 246 differences in timing of expression of the CBF and COR genes seen here will partially reflect 247 variation in the age of plants tested and the experimental conditions (summarised in 248 Additional file 7: Supplementary Notes) that here included the application of cold at dusk with 249 gradual reduction in temperature at the beginning of the cold treatment (13.5°C at 30 min, 250 8°C at 1 h, 5.5°C at 90 min and 4°C at 2 h; Additional file 1: Figure S6). On the other hand, 251 the differences in the detection of DE genes (Fig. 2b) most likely reflect the quality and depth 252 of our data, including comparisons that remove the effects of the time-of-day variation, a 253 much higher number of time-points, and robust statistical analysis (Additional file 4: Table  254 S3; Additional file 7: Supplementary Notes). 255 Differential alternative splicing in response to cold was analysed previously using an 256 algorithm to extract individual probe hybridisation data from microarrays of Arabidopsis 257 seedlings exposed to 4°C [12]. Comparison with TAIR9 transcripts identified over 200 DAS 258 transcripts although only half of those tested experimentally were validated [12]. While 259 demonstrating that AS occurred in response to cold, the limited resolution of this approach 260 only detected a fraction of DAS genes and transcripts. In comparison, the analysis here 261 identified 2,442 DAS genes and 4,105 DTU transcripts. In particular, we identified 1,500 262 novel cold response DAS genes of which 1,331 displayed regulation only by AS (DAS-only) 263  Table S7). Cluster 12 shows highly 301 increased gene level expression in the first 3 h of cold treatment reflecting the 302 reconfiguration of membranes in response cold to maintain fluidity and protect against 303 subsequent freeze damage [3]. 304 For the DAS genes, the most enriched functional terms were related to mRNA 305 splicing ( Fig. 2c; Additional file 8: Table S6). One hundred and sixty-six (7%) DAS genes 306 were RNA-binding proteins, spliceosomal proteins or SFs. Hierarchical clustering of the DTU AS isoform transcripts at different times during Day 1 at 4°C while cluster 4 (n = 326) shows 320 late up-regulation of transcripts on the fourth day at 4°C. Clear gain in rhythmic expression 321 of AS transcripts upon cold is seen in cluster 7 (n = 258). Cluster 8 (n = 233) includes 322 transcripts with increased expression within the first 3 h of cold treatment. The z-score scale 323 represents mean-subtracted regularized log-transformed TPMs. The coloured bars above 324 the heat map indicate whether samples were exposed to light (coloured) or dark (black) in 325 the 3 h before sampling. 326 327

Rapid cold-induced changes in AS accompany the major transcriptional responses 328
The high temporal resolution of the time-course allowed us to examine the relative dynamics  Table S8). The speed of induction is highlighted 337 by 648 and 388 genes showing significant DE and DAS after only 3 h of cold (T10) (Fig. 4a,  338 b), respectively, despite the gradual reduction in temperature which takes 2 h to reach 4°C 339 (Additional file 1: Figure S6). Notably, three-quarters (76.5%; 849) of the DAS genes 340 detected within 9 h of cold ( Fig. 1a, T10-T12) also had large AS changes with at least one 341 transcript having ΔPS >0.2 (Additional file 11: Table S9). A further indicator of the speed of 342 AS response and the sensitivity of some AS to reductions in temperature was demonstrated 343 by identifying those DAS genes which show isoform switches (IS), where the relative 344 abundance of different isoforms is reversed in response to cold (Additional file 1: Figures 345 S11 and S12). We developed the Time-Series Isoform Switch (TSIS) program to identify ISs 346 in time-series data [45]. Using TSIS, a total of 892 significant (p<0.001) ISs that involved two 347 abundant transcript isoforms (see Methods) were identified in 475 unique DAS genes (Fig. 348 5a; Additional file 12: Table S10). The ISs involved different types of AS events and 77% 349 involved at least one protein-coding transcript isoform (Fig. 5a, b; Additional file 12: Table  350 S10). These either generated isoforms which coded for different protein variants, or where 351 AS occurred in the 5' or 3'UTR, the transcript isoforms coded for the same protein, or one of 352 the transcripts was non-protein-coding (e.g. PTC-containing). TSIS determines the two time-353 points between which a significant isoform switch occurs and, consistent with the rapid 354 changes in AS, the majority (57%) occurred between 0-6 h following transfer to 4°C ( Fig. 5a; 355 Additional file 12: Table S10). Thus, immediately in response to lowering temperature, there 356 are waves of transcriptional and AS activity involving thousands of genes. Tukey t-tests were performed to compare each temperature reduction results against 20°C 390 control. Significant differences are labelled with asterisks (*, p<0.05; **, p<0.01; ***, 391 p<0.001). 392 393

Cold-regulated expression and AS of transcription and splicing factors 394
The waves of differential expression and AS in response to cold likely reflect regulation by 395 transcription factors (TFs) and splicing factors/RNA-binding proteins (SF-RBPs). Differential 396 expression of TFs in response to cold has been well documented (see Background). Here, 397 532 of the 2,534 TFs predicted in Arabidopsis (Additional file 13:  (Table 1). Thus, many TF and SF-RBP genes were regulated by AS in response to lower 403 temperatures. The majority have not previously been associated with the cold response and 404 represent putative novel cold response factors (Additional file 14: Table S12). We next 405 identified the TF and SF-RBP genes with the fastest (0-6 h after onset of cold) and largest 406 changes in expression and AS (log2 fold change ≥1.5 (equivalent to 2.83-fold change) for DE 407 genes and ΔPS >0.25 for at least one transcript in DAS genes) (Additional file 14: Table  408 S12). Fifty-nine TF and 47 SF-RBP DAS genes were identified as "early AS" genes. The TF 409 genes included a high proportion of circadian clock genes as well as genes associated with 410 abiotic stress, flowering time and hormone responses (Additional file 15: Table S13). The 411 SF-RBP genes included serine-arginine-rich (SR) and heterogeneous ribonucleoprotein 412 particle (hnRNP) protein genes known to be involved in stress responses and regulation of 413 the clock. For many of the early AS genes, the AS changes were either only observed in day 414 1 at 4°C, or persisted through the cold treatment, and many involved isoform switches 415 (Additional file 15: Table S13). Thus, both transcription and AS determine expression levels 416 of TFs and SFs and many of these genes were regulated rapidly in response to reduction in 417 temperature. 418 419

Speed and sensitivity of AS responses to small reductions in temperature 424
The rapid and large changes in AS suggest that many AS events are sensitive to relatively 425 small changes in temperature. To investigate further, we examined the effect on AS of early 426 AS genes when the temperature was lowered from 20°C to 18°C, 16°C, 12°C, 8°C and 4°C 427 for 12 h (Fig. 5c-e; Additional file 1: Figure S13). We observed that the relative levels of AS 428 transcript isoforms were dependent on the temperature. Three of the six genes analysed 429 showed significant changes in the level of at least one AS isoform with only a 2°C reduction 430 in temperature (to 18°C) while others were affected by 4°C or 8°C reductions ( Fig. 5e; 431 Additional file 1: Figure S13). We then examined the speed and sensitivity of AS by taking 432 multiple samples between 0 and 3 h after the onset of cold (

U2B"-LIKE is regulated at the AS level and is required for freezing tolerance 438
Many of the early AS genes, including TFs and SFs, showed large and rapid changes in AS 439 which alter the levels of protein-coding transcripts (Additional file 15: Table S13). We 440 hypothesised that such significant changes in AS in response to low temperature are 441 important in the overall cold acclimation process of the plant and lead to improved ability to 442 tolerate freezing conditions after acclimation. In support of this, four of the early AS genes 443 have previously been shown to be required for cold acclimation and tolerance to freezing:  (Table 2). To 445 examine whether other early AS genes may be involved in cold acclimation, we selected the 446 SF-RBP gene U2B"-LIKE because it was a novel DAS-only gene with an adaptive 447 expression pattern. We isolated a knockout mutant of the U2B"-LIKE gene (AT1G06960; Fig.  448 6a; Additional file 1: Figure S15). U2B"-LIKE has two main AS transcripts, the fully spliced 449 protein-coding mRNA and an isoform with retention of intron 4 (I4R) (P1 and P2, respectively 450 - Fig. 6a). In wild-type plants, the protein-coding P1 transcript isoform showed rhythmic 451 expression at 20°C, loss of rhythm during day 1 at 4°C, maintaining a high level of 452 expression throughout the remaining cold treatment (Fig. 6a). In freezing tolerance tests 453 conducted at -8.0°C and -8.5°C, the u2b"-like mutant plants showed greater sensitivity to 454 freezing; u2b"-like did not survive freezing at -8.5°C after cold acclimation while wild-type 455 plants recovered (Fig. 6b). u2b"-like mutant and WT plants both recovered at -8.0°C. 456 Differential sensitivity of the mutant was confirmed in quantitative electrolyte leakage 457 analyses; leaf tissue of u2b"-like suffered significantly increased ion leakage (cellular 458 damage) at -10°C than wild type plants (Fig. 6c) indicating that expression of U2B"-LIKE is 459 required for cold acclimation and freezing tolerance.  Arabidopsis contains two U2B"-related genes: U2B" (AT2G30260) and U2B''-LIKE 480 (AT1G06960). The two proteins are very similar: 80% identical and 90% similar at the 481 protein level (Additional file 1: Figure S16). U2B" is an U2snRNP-specific protein which 482 binds, along with U2A', to stem-loop IV of U2snRNA in both plants and human (Additional 483 file 7: Supplementary Notes). In the u2b"-like mutant, there was no expression of U2B"-LIKE 484 but expression of the U2B" paralogue (which was neither DE nor DAS in cold) (Additional file 485 1: Figure S15) was detected, suggesting that U2B" protein could not compensate for the lack 486 of U2B"-LIKE in the u2B"-like mutant and therefore that they had functionally diverged. To 487 investigate whether U2B"-LIKE affected AS regulation, we then compared AS patterns of 41 488 genes (including 34 DAS or DE+DAS genes identified here) in wild-type and u2b"-like 489 mutant plants. Five genes showed significantly different AS (p<0.05 and >10% difference in 490 splicing ratio between the mutant and wild-type) (Additional file 1: Figure S17; Additional file 491 16: Table S14). These included decreased levels of fully spliced, protein-coding transcripts 492 of PIF7 (AT5G61270; Additional file 1: Figure S17) which along with TOC1 and PHYB 493 represses expression of CBFs [40] and HOMOLOGUE OF HY5 (HYH), a clock input gene. 494 Thus, U2B"-LIKE is one splicing factor that contributes to correct splicing of PIF7 linking 495 U2B"-LIKE-dependent AS to regulation of the major cold response pathway. Therefore, the 496 freezing sensitivity of the u2b"-like mutant may be due to altered AS and expression of 497 specific genes required for cold acclimation (Additional file 7: Supplementary Notes). we demonstrate the dynamic contribution of AS by the rapid cold-induced wave of AS 506 activity accompanying the transcriptional response (Fig. 4) and the sensitivity of AS of some 507 genes to small reductions in temperature (Fig. 5). We also significantly demonstrate the 508 extent of AS by showing that over 2,400 genes are regulated by AS in response to cold with 509 over 1,600 regulated only at the AS level (Fig. 2). The massive changes in expression and 510 AS involved thousands of genes reflecting activation of both transcription and splicing 511 pathways and networks. The speed and extent of the cold-induced AS suggest that AS, 512 along with the transcriptional response, is a major driver of transcriptome reprogramming for 513 cold acclimation and freezing tolerance. 514 With over 2,400 genes regulated by AS, multiple different mechanisms are likely to 515 control the splicing decisions. Reduction in temperature to 4°C is expected to reduce the 516 rate of biochemical reactions and potentially affect transcription and splicing. We observed 517 that the vast majority of introns in the pre-mRNAs of all the cold-expressed genes are 518 efficiently spliced throughout the cold treatment. Therefore, low temperature does not cause 519 a general defect in splicing reflecting the ability of temperate plants to grow in a wide range 520 of fluctuating temperatures. Nevertheless, low temperatures may directly affect AS 521 regulation. For example, in mammals, secondary structures in pre-mRNAs affect splice site 522 selection [48] and cooling could stabilise such structures. Similarly, splicing is largely co-523 transcriptional and slower rates of RNA polymerase II (PolII) elongation promote selection of 524 alternative splice sites [49]. Both of these mechanisms will undoubtedly be involved in the 525 cold-induced AS changes of some of the genes seen here. However, the sensitivity of AS to 526 reductions in temperature of only a few degrees and clear rhythmic expression profiles of AS 527 transcript isoforms in plants exposed to constant 4°C temperature for four days (e.g. cluster 528 7 in Fig. 3 and Additional file 1: Figure S10) argue against such mechanisms being widely 529 responsible for the cold-induced AS changes observed here. Local or global DNA 530 methylation and chromatin modifications can also affect the rate of PolII elongation or help to 531 recruit SFs to affect splice site selection [49]. In plants, epigenetic regulation is responsible 532 for suppression of FLC by vernalisation in the seasonal response to cold [50]. Furthermore, 533 altered histone 3 lysine 36 tri-methylation (H3K36me3) was recently shown to affect some 534

AS events induced by higher ambient temperatures within 24 h [14]. Alongside dynamic 535
changes in histone marks at specific stress-induced genes [51,52] it is likely that some of 536 the cold-induced AS here reflects local epigenetic changes. 537 We showed that the levels of hundreds of TF and SF-RBP gene transcripts changed 538 in response to cold at both the transcriptional and AS levels. Therefore, splicing decisions in 539 the physiological response to low temperature are most likely controlled by altered 540 abundance, activity or localisation of SFs or other RNA-interacting proteins [11, 18, 19, 24, 541 53]. In particular, we identified TF and SF-RBP genes with large and rapid changes in AS. 542 Most of the early AS transcription factor genes were regulated only by AS and, therefore, 543 had not been identified previously as cold response transcription factors. Nevertheless, the 544 rapid cold-induced changes in the AS of some known cold response genes: CAMTA3 which 545 activates the CBFs [54], and VRN2 and SUF4 which are involved in vernalisation and 546 silencing of FLC [43, 50], have not been described previously and our results introduce AS 547 as a novel component in their regulation. It will be interesting to address the function of the 548 novel AS-regulated TFs and the function of AS of these and known cold response TFs on 549 cold acclimation and vernalisation in future experiments. 550 In contrast, the early AS SF-RBP genes included SR and hnRNP protein genes 551 known to respond to changes in temperature (e.g. SR30, RS40, GRP8, SR45A, PTB1, 552 RBP25 etc.) [11,23,24]. Many SF-RBP genes are regulated by AS-NMD and the rapid 553 induction of AS in these and other early AS genes affects the abundance of protein-coding 554 transcripts and presumably of the splicing factors themselves to alter AS of downstream 555 targets. Various spliceosomal and snRNP protein genes are also among the early AS genes. 556 These include GEMIN2 (snRNP assembly) which is cold-induced, involved in regulation of 557 the circadian clock, and enhances U1snRNP assembly to compensate for reduced 558 functionality of U1snRNP at low temperatures [22]. Interestingly, a number of U1snRNP core 559 and associated protein genes (U1-70k, LUC7B, LUC7RL, PRP39A, RBM25) [55, 56] 560 respond rapidly to cold via AS. The early AS genes also include two wound-induced RNA-561 binding proteins, UBA2a and UBA2c [57] and may also therefore be involved in the cold 562 response. Three LAMMER kinase genes (AFC1, AFC2 and AME3) which regulate SR 563 proteins via phosphorylation showed changes in their expression due to AS suggesting that 564 lower temperatures affect activation/deactivation of specific splicing factors which are targets 565 of these kinases [47,58]. In addition, over 20 putative RNA-binding proteins, kinases and 566 RNA helicases with little or no known function are among the novel early AS genes. Four of 567 the early AS SF-RBP genes (RCF1, STA1, GEMIN2 and AME3; Table 2) have been shown 568 to be involved in freezing tolerance [22,47] and we provide initial evidence for another early 569 AS gene, U2B"-LIKE, being involved in freezing tolerance and acclimation ( Fig. 6; Additional 570 file 7: Supplementary Notes). Our results identify over 100 splicing and transcription 571 regulatory genes, whose expression is rapidly and drastically altered by AS in response to 572 cooling. Future work will address the function of these putative regulators and specific 573 transcript isoforms in cold acclimation. 574 The speed of change of AS may be one of the earliest responses to cooling. We 575 showed significant AS within only 40-60 minutes of cooling and with subtle reductions in 576 temperature of as little as 2°C (Fig. 5). Similar responses are seen in mammals where 577 neuronal stimulation and rapid changes of intracellular sterols also activate splicing/AS 578 within minutes without de novo transcription or protein production [59,60]. In addition, a 1°C 579 change in body temperature activated a program of AS changes within 30 min which 580 involved temperature-sensitive phosphorylation of SR proteins and AS of the U2-associated 581 factor, U2AF26 [61]. Therefore, cold-induced AS programs may similarly involve rapid 582 phosphorylation/dephosphorylation of SFs such as SR proteins to modulate AS [61,62].

Plant Material and Growth Conditions 654
Arabidopsis thaliana Col-0 seeds were surface sterilized, stratified in the dark at 4°C for 4 655 days, and grown hydroponically to maturity (5 weeks) in Microclima environment-controlled 656 cabinets (Snijders Scientific), maintaining 20°C, 55% relative humidity, and 12 h light (150 657 μE m −2 s 1 ):12 h dark as described previously [72,73]. 10-13 Arabidopsis rosettes were 658 harvested and pooled at each sampling time. Harvesting occurred every 3 h over the last 24 659 h at 20°C, and on days 1 and 4 after transfer to 4°C giving 26 time-points in the time-series 660 (Fig. 1a). Day 1 at 4°C represents the "transition" from 20°C to 4°C when plants first begin to 661 experience the temperature decrease; Day 4 at 4°C represents "acclimation" where plants 662 have been exposed to 4°C for 4 days (Fig. 1a). Three biological replicates were generated 663 for each time-point in separate experiments (78 samples in total). The same growth cabinet 664 was used for all repeats to eliminate the potential effects of minor changes in light intensities 665 and light quality on gene expression. Additionally, to avoid interference in the experiment 666 from possible microclimates within the growth cabinet, 26 trays for each time-point were 667 placed in a randomised fashion. The switch to 4°C from 20°C was initiated at dusk. In a 668 temperature reduction, the cabinet used here typically takes 2 h to reach 4°C air 669 temperature (Additional file 1: Figure S6). Tissue was rapidly frozen in liquid N2 and stored at 670 −80°C until isolation of RNA and preparation of cDNA. 671

RNA Extraction 673
Total RNA was extracted from Arabidopsis tissue using RNeasy Plant Mini kit (Qiagen), 674 followed by either on-column DNase treatment (for HR RT-PCR, see below), or the TURBO 675 DNA-free™ kit (Ambion) (for library preparation and qRT-PCR, see below).  Table S15). 687 Residual adaptor sequences at both 5' and 3' ends were removed from raw reads using 688 cutadapt version 1.4.2 (https://pypi.python.org/pypi/cutadapt/1.4.2) with quality score 689 threshold set at 20 and minimum length of the trimmed read kept at 20. The "--paired-output" 690 option was turned on to keep the two paired read files synchronized and avoid unpaired 691 reads. The sequencing files before and after the trimming were examined using fastQC 692 version 0.10.0. 693 694

Quantification of transcripts and AS 695
Arabidopsis transcript expression from our RNA-seq experiment was carried out using 696 Salmon version 0.82 [26] in conjunction with AtRTD2-QUASI augmented by 8 genes that 697 were not originally present [27]. For indexing, we used the quasi mapping mode to build an 698 auxiliary k-mer hash over k-mers of length 31 (--type quasi -k 31). For quantification, the 699 option to correct for the sequence specific bias ("--seqBias") was turned on. The number of 700 bootstraps was set to 30 and all other parameters were on default settings. AtRTD2-QUASI 701 is a modification of AtRTD2, a high quality reference transcript dataset for Arabidopsis 702  Table S16. 707 708

Differential gene expression (DE) and AS (DAS) analysis of the RNA-seq data 709
To carry out differential expression analysis, transcript quantification results generated by 710 Salmon were processed and refined in successive steps (Additional file 1: Figure S18). First, 711 transcript and gene read counts were generated from TPM data correcting for possible gene 712 length variations across samples using tximport version 0.99.2 R package with the option 713 "lengthScaledTPM" [74]. Second, read count data from sequencing replicates were summed 714 for each biological sample. Third, genes and transcripts that are expressed at very low levels 715 were removed from downstream analysis. The definition of a low expressed gene and 716 transcript was determined by analysing mean-variance relationships [30]. The expected 717 decreasing trend between the means and variances was observed in our data when 718 removing transcripts that did not have ≥ 1 counts per million (CPM) in 3 or more samples out 719 of 78, which provided an optimal filtering for low expression transcripts. At the gene level, if 720 any transcript passed the expression level filtering step, the gene was included as an 721 expressed gene and then the normalisation factor, which accounted for the raw library size, 722 was estimated using the weighted trimmed Mean of M values method using edgeR version change of expression of the other transcripts were tested using a t-test for the same 18 743 contrasts using the DiffSplice function (Additional file 1: Figure S2). DTU transcripts which 744 coded for protein isoforms were determined from transcript translations of AtRTD2 [27]. 745 Genes were significantly DE at the gene level if they had at least two contrast groups 746 at consecutive time-points with adjusted p<0.01 and ≥2-fold change in expression in each 747 contrast group (Additional file 1: Figure S1). Genes/transcripts with significant DAS/DTU, 748 had at least two consecutive contrast groups with adjusted p<0.01 and with these contrast 749 groups having at least one transcript with ≥10% change in expression (Additional file 1: 750 Figure S2). Gene functional annotation was performed using R package 751 RDAVIDWebService version 1.8.0 [78][79][80]. The possibility of a gene and transcript being 752 identified by accident as cold responsive by the statistical method was tested. Time-points 753 T1 and T9 (Fig. 1a) are virtually identical as they both represent dusk samples at 20°C, the 754 only difference being they are 24 hours apart, such that few DE or DAS genes and DTU 755 transcripts were expected when comparing these time-points. Indeed, no significant DE or 756 DAS gene, nor DTU transcript, was identified between T1 and T9. This suggests our 757 statistical method to select cold-responsive genes and transcripts is conservative and 758 controls the number of false positives. 759 760 Identification of isoform switches 761 5,317 high abundance transcripts, whose average expression accounts for >20% of total 762 gene expression at at least one time-point, were selected from DAS gene transcripts for the 763 isoform switch analysis using the TSIS R package which is a tool to detect significant 764 transcript isoform switches in time-series data [45]. Switches between any two time-points 765 were identified by using the default parameters in which i) the probability of switch (i.e. the 766 frequency of samples reversing their relative abundance at the switches) was set to >0.5; ii) 767 the sum of the average differences of the two isoforms in both intervals before and after the 768 switch point were set at ΔTPM>1; iii) the significance of the differences between the 769 switched isoform abundances before and after the switch was set to p<0.001; and iv) both 770 intervals before and after switch must consist of at least 2 consecutive time-points in order to 771 detect long lasting switches. SUPPA version 2.1 [46] was then used to identify the specific 772 AS events (e.g. intron retention, alternative 3' or 5' splice site selection, exon skip) that 773 distinguished the pair of switch transcript isoforms. 774 775

High-resolution (HR) RT-PCR 787
HR RT-PCR reactions were conducted as described previously [82]. Gene-specific primer 788 pairs were used for analysing the expression and alternative splicing of different genes 789 (Additional file 19: Table S17). For each primer pair, the forward primer was labelled with 6-790 carboxyfluorescein (FAM). cDNA was synthesised from 4 μg of total RNA using the Sprint 791

RT Complete -Double PrePrimed kit following manufacturer's instructions (Clontech 792
Laboratories, Takara Bio Company, USA). The PCR reaction usually contained 3 μL of 793 diluted cDNA (1:10) as a template, 0.1 μL of each of the forward and reverse primers (100 794 mM), 2 μL of 10 X PCR Buffer, 0.2 μL of Taq Polymerase (5U/μL, Roche), 1 μL of 10 mM 795 dNTPs (Invitrogen, Life Technologies) and RNase-free water (Qiagen) up to a final volume 796 of 20 μL. For each reaction, an initial step at 94 o C for 2 min was used followed by 24-26 797 cycles of 1) denaturation at 94 o C for 15 sec, 2) annealing at 50 o C for 30 sec and 3) 798 elongation at 70 o C for either 1 min (for fragments smaller than 1000 bp) or 1.5 min (for 799 fragments between 1000-1200 bp) and a final extension cycle of 10 min at 70 o C. 800 To separate the RT-PCR products, 1.5 μL of PCR product was mixed with 8.5 μL of Hi-Di TM 801 formamide (Applied Biosystems) and 0.01 μL of GeneScan TM 500 LIZ TM dye or 0.04 μL of 802 GeneScanTM 1200 LIZ TM dye size standard and run on a 48-capillary ABI 3730 DNA 803 Analyser (Applied Biosystems, Life Technologies). PCR products were separated to single 804 base-pair resolution and the intensity of fluorescence was measured and used for 805 quantification in Relative Fluorescent Units (RFU). The different PCR products and their 806 peak levels of expression were calculated using the Genemapper® software (Applied 807 Biosystems, Life Technologies). 808 809 Identification and characterisation of the u2b''-like mutant 810 cDNA was synthesised as described above for HR RT-PCR. PCR was performed using 811 cDNAs and GoTaq Green DNA polymerase (Promega) following manufacturer's instructions. 812 Primer sequences are provided in Additional file 19: Table S17. 813 814

Freezing and electrolyte leakage assay 815
Cold acclimated plants were assessed for damage after freezing conditions. Sterilised seeds 816 were sown on MS-agar plates and after 7 days seedlings were transferred to peat plugs for 817 growth in 12:12 light-dark cycles, 150 to 200 μE/m 2 /s at 20°C for 4 weeks. Plants were then 818 transferred to 5°C, 10:14 LD cycles (150 μE/m 2 /s) for ca. 14 d after which they were used in 819 either a qualitative or quantitative assay. In the qualitative assay, cold-acclimated plants 820 were transferred at dusk to either -8.0°C or -8.5°C for 24 h, then transferred to 5°C, 10:14 821  Description of data: Gene lists of known Arabidopsis cold-response genes in Table S5.   Table S8a. Gene lists of 7,302 DE genes organised by the time-point 940 at which they first become significantly differentially expressed upon cold. Table S8b. Gene 941 descriptions of the 648 genes which are first detected as significantly differentially expressed 942 at 3 hours after the start of cold treatment (T10) in Table S8a. Table S8c. Gene descriptions 943 of the 1452 genes which are first detected as significantly differentially expressed at 6 hours 944 after the start of cold treatment (T11) in Table S8a. Table S8d. Gene lists of 2,442 DAS 945 genes organised by the time-point at which they first become significantly differentially 946 alternatively spliced upon cold. Table S8e. Gene descriptions of the 388 genes which are 947 first detected as significantly differentially alternatively spliced at 3 hours after the start of 948 cold treatment (T10) in Table S8d. Table S8f. Gene descriptions of the 768 genes which are 949 first detected as significantly differentially alternatively spliced at 6 hours after the start of 950 cold treatment (T11) in Table S8d Table S10a. Isoform switch analysis of transcripts from DAS genes 963 using the TSIS package . Table S10b. AS events identified with SUPPA 964 (Alamancos et al., 2015) for each isoform switch. Table S10c. Gene lists of DAS genes with 965 Isoform Switches (IS) between 0-3 h cold (T9-T10), 3-6 h cold (T10-T11) and 0-6 h cold (T9-966 T11 at the alternative splicing level (DAS) and TF genes with the earliest and largest changes in 979 AS. Table S12b. Summary of 166 SF-RBP genes differentially regulated by cold at the 980 alternative splicing level (DAS) and SF-RBP genes with the earliest and largest changes in 981 AS. Table S12c. TF and SF-RBP DAS genes with the largest and most rapid changes in AS. 982 Table S12d. TF and SF-RBP DE genes with the largest and most rapid changes in 983 expression. Table S12e. Summary of functions/classifications of TF genes showing the 984 largest and most rapid changes in gene expression and/or AS in cold treatment. Table S12f. 985 Classifications of SF-RBPs showing the largest and most rapid changes in gene expression 986 and/or AS in cold treatment. 987