- © 2018 The Authors.
Abstract
Plants have adapted to tolerate and survive constantly changing environmental conditions by reprogramming gene expression The dynamics of the contribution of alternative splicing (AS) to stress responses are unknown. RNA-sequencing of a time-series of Arabidopsis thaliana plants exposed to cold determines the timing of significant AS changes. This shows a massive and rapid AS response with coincident waves of transcriptional and AS activity occurring in the first few hours of temperature reduction and further AS throughout the cold. In particular, hundreds of genes showed changes in expression due to rapidly occurring AS in response to cold (“early AS” genes); these included numerous novel cold-responsive transcription factors and splicing factors/RNA binding proteins regulated only by AS. The speed and sensitivity to small temperature changes of AS of some of these genes suggest that fine-tuning expression via AS pathways contributes to the thermo-plasticity of expression. Four early AS splicing regulatory genes have been shown previously to be required for freezing tolerance and acclimation; we provide evidence of a fifth gene, U2B”-LIKE. Such factors likely drive cascades of AS of downstream genes that, alongside transcription, modulate transcriptome reprogramming that together govern the physiological and survival responses of plants to low temperature.
INTRODUCTION
Plants adapt to and survive adverse environmental conditions by reprogramming their transcriptome and proteome. Low temperatures negatively affect growth and development, but, in general, plants from temperate climatic regions can tolerate chilling temperatures (0–15°C) and can increase their freezing tolerance by prior exposure to low, nonfreezing temperatures through a process called cold acclimation. Cold acclimation involves multiple physiological, biochemical, and molecular changes that reduce growth, modify membrane properties, and produce the proteins and metabolites required to protect cell integrity throughout the period of freezing exposure (Thomashow, 2010; Knight and Knight, 2012; Zhu, 2016). In Arabidopsis thaliana, these cold-responsive changes reflect complex reprogramming of gene expression involving chromatin modification, transcription, posttranscriptional processing, posttranslational modification, and protein turnover (Thomashow, 2010; Knight and Knight, 2012; Barrero-Gil and Salinas, 2013; Kim et al., 2015; Zhu, 2016). Previous genome-wide microarray and RNA-seq studies have mostly focused on the cold-induced transcriptional response and thousands of differentially expressed genes have been reported. The best studied is the signaling cascade involving the C-REPEAT BINDING FACTOR/DEHYDRATION RESPONSE ELEMENT BINDING PROTEIN (CBF/DREB) transcription factors (CBF1-3), which recognize the C-repeat/dehydration responsive element motifs in the promoters of their target genes (Thomashow, 2010; Knight and Knight, 2012; Jia et al., 2016; Zhao et al., 2016). Increased expression of these genes during cold acclimation increases freezing tolerance, which can also be achieved by ectopic expression of the CBFs in the absence of cold acclimation (Thomashow, 2010; Knight and Knight, 2012; Jia et al., 2016; Zhu, 2016). Besides the CBF regulon, other transcriptional pathways are required for the activation of cascades of gene expression that together drive the cold response, acclimation, and plant survival (Park et al., 2015). Plants also adapt to daily and seasonal changes in temperature, and the circadian clock provides a mechanism by which they anticipate the predictable diel cycles of temperature and light/dark to optimize gene expression and consequent physiology (Harmer, 2009). The circadian clock also regulates cold-responsive gene expression through a process called gating where the magnitude of changes in gene expression depends on the time of day that the stimulus is applied (Harmer, 2009; Thomashow, 2010).
Alternative splicing (AS) has been linked to stress responses (Mastrangelo et al., 2012; Leviatan et al., 2013; Staiger and Brown, 2013; Hartmann et al., 2016; Klepikova et al., 2016; Li et al., 2016; Pajoro et al., 2017; Laloum et al., 2018), but little is known about its global contribution or dynamics. AS is a regulated process that produces different mRNA transcripts from precursor mRNAs (pre-mRNAs) of a single gene (Fu and Ares, 2014; Lee and Rio, 2015). The selection of splice sites is determined by sequence elements in the pre-mRNA that are recognized by transacting splicing factors (SFs), which recruit the spliceosome for intron removal. The concentration, localization, and activity of these factors in different cells and conditions defines splice site choice to generate different alternatively spliced isoforms and splicing patterns (Fu and Ares, 2014; Lee and Rio, 2015). AS regulates the level of expression of genes by generating nonproductive transcript isoforms that are targeted for nonsense-mediated decay (NMD) or impacts the proteome by generating protein variants with altered sequence and function. In higher plants, extensive AS variation has been observed in different cell types, organs, and developmental stages (Wang et al., 2014; Klepikova et al., 2016; Li et al., 2016; Vaneechoutte et al., 2017) and has been implicated in a wide range of physiological processes including responses to stress (Mastrangelo et al., 2012; Leviatan et al., 2013; Staiger and Brown, 2013; Hartmann et al., 2016; Klepikova et al., 2016; Li et al., 2016; Pajoro et al., 2017). Around 60% of Arabidopsis intron-containing genes undergo AS (Marquez et al., 2012). SFs are required for normal growth and development including control of flowering time, regulation of the circadian clock, and stress responses, suggesting that regulated AS of downstream targets is essential (Zhang and Mount, 2009; Reddy et al., 2013; Staiger and Brown, 2013; Schlaen et al., 2015; Cheng et al., 2017; Laloum et al., 2018). The importance of AS to the cold response has been demonstrated by altered cold sensitivity or tolerance when SFs are misexpressed (Reddy et al., 2013; Staiger and Brown, 2013; Laloum et al., 2018). These studies strongly suggest that AS networks are central coordinators of the cold response. However, virtually nothing is known about the extent and timing of the contribution of AS or how transcription and AS determine the dynamic changes in the transcriptome required for response to cooling, acclimation, freezing tolerance, and survival.
RNA-sequencing (RNA-seq) provides the means to understand how AS impacts gene expression and transcriptome reprogramming. However, the quality of AS analysis depends on accurate quantification of transcripts. Transcript assembly from short reads is often inaccurate due to misassembly of transcripts and nonassembly of bona fide transcripts both of which impact the accuracy of transcript quantification (Steijger et al., 2013; Hayer et al., 2015; Zhang et al., 2017). In addition, accuracy also decreases with increasing numbers of isoforms in a gene (Hayer et al., 2015). Rapid, nonalignment programs, Salmon or kallisto, have high accuracy of transcript quantification (Bray et al., 2016; Patro et al., 2017) but require well annotated and comprehensive transcriptomes, which are often limited in plant species. To improve the accuracy of transcript abundance and AS in RNA-seq analyses in Arabidopsis, we developed a comprehensive Reference Transcript Data Set (AtRTD2; Zhang et al., 2017) for use with Salmon or kallisto. AtRTD2 contains over 82,000 transcripts giving greatly increased coverage of AS transcript isoforms than currently available in the Arabidopsis TAIR10 and Araport11 transcriptome data sets. Importantly, we demonstrated the increased accuracy of AS measurements using Salmon/AtRTD2 and the negative impact of missing transcripts (Zhang et al., 2017). Although AS has been detected in many Arabidopsis genes in response to stress conditions (Filichkin et al., 2010; Staiger and Brown, 2013; Ding et al., 2014), the true scale and the dynamics of the AS response need to be addressed.
High-resolution RNA-seq time courses provide insights into dynamic changes in expression during, for example, stomatal or apical meristem development (Adrian et al., 2015; Klepikova et al., 2015). The value of RNA-seq time-series analysis to understand relative levels of gene expression and AS is illustrated in development and differentiation in animal systems (Cieply et al., 2016; Gloss et al., 2017; White et al., 2017). Here, we exploit time-series RNA-seq to examine the dynamics of differential AS alongside differential expression in response to lowering temperatures taking full account of the time-of-day variations in gene expression due to photoperiod and circadian control. The increased resolution from ultra-deep RNA-seq of multiple time points identified 8949 genes that were differentially expressed (DE) at the gene level and differentially alternatively spliced (DAS). These included 1647 genes that were regulated only at the AS level, the majority of which had not been identified as cold responsive by other methods. The high temporal resolution of both gene expression and AS shows the rapid and dynamic induction of changes of both total gene expression and AS of thousands of genes in the first few hours after the onset of cold. Response to low temperature is thus controlled by genome-wide changes in both transcription and AS. Early cold response genes include specific splicing and transcription factors that undergo rapid AS that is sensitive to small changes in temperature. Our data suggest a mechanism whereby dynamic changes in AS of splicing regulators contribute to the changes in the transcriptome that condition the developmental and physiological processes required by plants to respond to constantly changing environmental conditions.
RESULTS
Almost Half of the Expressed Genes Show Cold-Induced Differential Expression and/or Differential AS
To examine changes in gene expression and AS in response to low temperature, we performed deep RNA-seq on a diel time series of rosettes of 5-week-old Arabidopsis Col-0 plants grown at 20°C and transferred to 4°C (Figure 1A). Rosettes were sampled at 3 h intervals for the last day at 20°C, the first day at 4°C, and the fourth day at 4°C as described in Figure 1A (see Methods). The experiment was repeated separately with different batches of plants on three occasions giving three biological repeats per time point. We generated over 360 million 100-bp paired end reads for each of the 26 time points and quantified transcript abundance in transcripts per million (TPM) using Salmon (Patro et al., 2017) and AtRTD2-QUASI as reference transcriptome (Zhang et al., 2017), allowing us to determine patterns of expression at both the gene level (sum of all transcript abundances of a gene) and at the individual transcript isoform level (Supplemental Figures 1 and 2). Principal component analysis of the gene level expression data from across the 26 time points showed that temperature (32.4% of variance) and time of day (18.9% of variance) were the overwhelming drivers of gene expression (Figure 1B). To analyze the time-series data at gene and transcript levels to obtain DE, DAS, and differential transcript usage (DTU) results, we developed an analysis pipeline that exploited the general linear models available in limma (Law et al., 2014, 2016; Ritchie et al., 2015; see Methods), which allowed biological repeats and time-of-day variation in expression to be taken into account in the statistical analysis to produce more accurate and robust results. It is important to note that in the time-series analysis, we compared gene and transcript abundances between equivalent time points at 20°C and those in day 1 and day 4 at 4°C to remove the effects of time of day so that the changes detected are due to reduction in temperature (contrast groups in Supplemental Figures 1 and 2). We first analyzed differential expression at the gene level (DE). A gene was considered differentially expressed if it showed a log2 fold change ≥1 (≥2-fold change) in expression in at least two consecutive contrast groups (adjusted P < 0.01; Figure 2A; Supplemental Figure 1 and Supplemental Data Set 1). Using these stringent criteria, we identified a total of 7302 genes that were significantly differentially expressed in response to low temperature when compared with 20°C levels. Of these, 48.2% were upregulated and 51.8% downregulated (Supplemental Figure 3 and Supplemental Data Set 2).
Analyses of the Arabidopsis Response to Low Temperature in Diel Conditions.
(A) Schematic representation of sampling strategy. Time points of sampling are marked by vertical colored lines and labeled from 1 to 26. Five-week-old Arabidopsis rosettes were harvested every 3 h over a 24-h period at 20°C (vertical red lines). At dusk, the temperature was gradually reduced to 4°C, and harvesting continued during the first day at 4°C (vertical blue lines) and the fourth day at 4°C (vertical green lines). Black boxes, 12 h dark; white boxes, 12 h light.
(B) Principal component analysis of the Arabidopsis rosette transcriptome before and after a shift from 20°C to 4°C. Each data point (T1–T26) refers to one time point and represents the average gene expression (n = 3) from the RNA-seq data. The data points are connected by arrows in chronological order: black for 3 h of dark and yellow for 3 h of light. The dotted gray line joining T17 and T18 represents days 2 and 3 at 4°C.
DE and DAS Analyses of Arabidopsis Response to Low Temperature.
(A) Flow chart showing the distribution of the 8949 DE (blue) and DAS (red) genes (Supplemental Data Set 1). The DE and DAS gene sets are largely different with only 795 (11.26%) in common (overlap between blue and red circles in [B]).
(B) Euler diagram of DE (left) and DAS (right) genes identified here and compared with known cold response DE/DAS genes (dashed circle). Information on known cold response DE/DAS genes can be found in Supplemental Data Sets 3 and 4, respectively.
(C) Most significantly enriched GO terms for DE (shades of blue) and DAS (shades of red) genes. Bar plots of -log10 transformed FDR values are shown. BP, biological process; MF, molecular function; CC, cellular component.
(D) Number and percentage of DE genes that are significantly DE in both day 1 and day 4 at 4°C and unique to day 1 and day 4.
(E) Number and percentage of DAS genes that are significantly DAS in both day 1 and day 4 at 4°C and unique to day 1 and day 4.
Second, we used the transcript-level data to identify genes that were DAS between contrast groups. We examined the consistency of expression changes among the transcripts and the corresponding genes (see Methods; Supplemental Figure 2) to detect DAS genes. Criteria for genes being DAS were that in at least two consecutive contrast groups (1) at least one of the transcripts differed significantly from the gene with an adjusted P value of <0.01, and (2) at least one of the transcripts of the gene showed a Δ percent spliced (ΔPS) of ≥0.1 [to keep only genes where there is a significant change coming from transcript(s) with large differences in their relative abundance] (Figure 2A; Supplemental Figure 2 and Supplemental Data Set 1). We identified 2442 DAS genes (Figure 2A), of which 795 were also DE genes (regulated by both transcription and AS) and 1647 genes that were not DE (regulated by AS only). Thus, of the 8949 genes that exhibited significantly altered levels of differential gene and/or transcript level expression, 27.3% were differentially alternatively spliced, consistent with widespread DAS in response to cold (Figure 2A; Supplemental Figure 2 and Supplemental Data Set 1). At any particular time point, between ∼600 and 3700 genes were DE and ∼400 to 1450 genes were DAS when compared with 20°C levels (Supplemental Figure 3).
To pinpoint the individual transcripts responsible for a gene being identified as a DAS gene, a DTU analysis was performed by testing the expression changes of every transcript against the weighted average of all the other transcripts of the gene using a t test (Supplemental Figure 2). DTU transcripts were identified as those that differed from the gene level and showed a ΔPS of ≥0.1 in at least two consecutive contrast groups with an adjusted P value of <0.01. In total, 34% (4105) of expressed transcripts of DAS genes were classed as DTU (Supplemental Data Set 1) of which ∼70% were protein-coding and ∼30% contained premature termination codons (PTCs). Time points T1 and T9 (Figure 1A) are virtually identical as they both represent dusk samples at 20°C, the only difference being they are 24 h apart. This contrast group acted as a control for the analysis and no significant DE or DAS gene, nor DTU transcript, was identified in the comparison between T1 and T9.
We next explored differences between immediate responses upon transfer to low temperature and the response to prolonged cold acclimation by comparing which genes were DE and DAS in day 1 and/or day 4 after transfer to 4°C. Around 50% (3573 genes) and 60% (1440) of DE and DAS genes, respectively, were common to both days with the remainder being either unique to day 1 or day 4 (Figures 2D and 2E; Supplemental Data Set 1). Thus, changes in gene-level expression and AS occurred throughout the cold period: either transiently (occurring in day 1 at 4°C and returning to 20°C levels by day 4), persisting throughout the cold period (occurring in day 1 at 4°C and remaining at day 4), or occurring later, only in day 4. We propose that these patterns of gene-level expression and AS reflect different contributions to low temperature perception, initial cold responses, and physiological acclimation to cold and freezing temperatures.
DE and DAS Analyses Identify Novel Cold Response Genes and AS of Cold Regulators
Previous analyses of differential gene expression in wild-type Arabidopsis plants exposed to cold used microarrays (Vogel et al., 2005; Carvallo et al., 2011; Leviatan et al., 2013) and, more recently, RNA-seq at the gene level (Gehan et al., 2015; Jia et al., 2016; Zhao et al., 2016) (Supplemental Data Sets 3 and 4, respectively). There was a substantial overlap between the cold response DE genes identified in those studies and the DE genes identified here (Figure 2B). Crucially, we identified an additional 1708 novel DE cold response genes (Supplemental Data Set 5). We also analyzed the 3594 genes identified in previous studies but which were not differentially expressed here. Almost one-quarter did not meet our expression level criteria (>1 count per million in at least 3 of the 78 time-point replicates) and three-quarters did not meet the stringent criteria of at least two consecutive contrast groups with a log2 fold change of >1 and adjusted P value of <0.01.
As expected, we showed cold induction of CBFs and selected COLD-REGULATED (COR) genes (none of which undergo AS) (Figures 3A and 3B). The first significant increase in expression of the CBF genes was detected between 3 and 6 h after the onset of cold treatment (applied at dusk). Expression of CBFs has been detected within <1 h and peaked at around 1 to 2 h in other studies conducted in constant light (Gilmour et al., 1998). Here, all three CBFs had significant fold-change increases between 0 and 3 h consistent with early induction but P values were above 0.01. The relative timing of expression of other transcription factors reflected previous studies. For example, some of the earliest “first wave” transcription factors (TFs; Park et al., 2015), ZAT10, ERF11, WRKY33, and CRF3, were DE 3 h after onset of cold treatment (see Supplemental Data Set 13B). Also consistent is ZAT12, which coregulates the CBF pathway and negatively regulates CBF expression (Vogel et al., 2005; Park et al., 2015) and is significantly DE at 3 to 6 h of cold along with the CBFs. CBF2 negatively regulates CBF1 and 3 (Novillo et al., 2004) and the timing of the increase in expression of CBF2 may indicate such regulation (Figure 3A). Differences in timing of expression of cold response genes seen here likely reflect differences in experimental conditions (Supplemental Data Set 3) and stringent statistical analyses as well as gating of cold responsive expression by the circadian clock (Fowler et al., 2005). The key experimental differences were the age of plants (5-week-old versus 10- to 14-d-old seedlings) and the application of cold at dusk with gradual reduction in temperature at the beginning of the cold treatment (13.5°C at 30 min, 8°C at 1 h, 5.5°C at 90 min, and 4°C at 2 h; Supplemental Figure 4). Application of cold at dusk instead of in the light should avoid additional effects due to redox sensing-signaling associated with the photosynthetic apparatus (Hüner et al., 2013).
Expression Profiles of CBFs, Selected COR Genes, and Known Cold Response Genes.
(A) Single-exon and single-transcript CBF genes.
(B) Total gene expression level of three COR regulatory genes: COR78 and single-transcript genes KIN1 and COR47.
(C) Expression of gene and protein-coding transcripts of RCF1 (AT1G20920). The P1, P3, P4, P5, P6, and P7 transcript isoforms of RCF1 differ at their 3′UTR sequences.
(D) PIF7 (AT5G61270) has protein-coding P1 and intron-retained P2 as main expressed transcripts.
(E) The novel DAS gene PHYB (AT2G18790) has two main transcripts, protein-coding P1 and intron-retained JS1.
(F) The novel DAS gene SUF4 (AT1G30970) has two main transcripts, protein-coding P1 and intron-retained s2.
(C) to (F) Rapid and significant isoform switches detected by TSIS (Guo et al., 2017) are labeled with red circles. For clarity, transcripts whose expression was below 4 TPM at all time points were not included in the graphs. Error bars indicate se of the mean, n = 3 biological replicates × 3 sequencing replicates.
Differential alternative splicing in response to cold was analyzed previously using an algorithm to extract individual probe hybridization data from microarrays of Arabidopsis seedlings exposed to 4°C (Leviatan et al., 2013). The study only identified around 200 DAS transcripts and only half of those tested experimentally were validated (Leviatan et al., 2013). While demonstrating that AS occurred in response to cold, the limited resolution of this approach only detected a fraction of DAS genes and transcripts. In comparison, the analysis here identified 2442 DAS genes and 4105 DTU transcripts. In particular, we identified 1500 novel cold response DAS genes of which 1331 displayed regulation only by AS (DAS-only) (Figure 2B; Supplemental Data Set 5). Among the DAS genes identified here and that had previous evidence of involvement in the cold response, we observed dynamic changes in AS in, for example, REGULATOR OF CBF EXPRESSION1 (RCF1), PHYTOCHROME-INTERACTING FACTOR7 (PIF7), PHYTOCHROME B (PHYB), and SUPPRESSOR OF FRIGIDA4 (SUF4) (Figures 3C to 3F). RCF1 encodes a cold-inducible RNA helicase required for cold tolerance (Guan et al., 2013). It produces transcript isoforms that differ by AS of introns in the 3′ untranslated region (UTR), which may cause particular isoforms to be retained in the nucleus or trigger NMD to regulate RCF1 expression at different temperatures. PIF7 is a transcriptional repressor of the circadian expression of CBFs and is involved in photoperiodic control of the CBF cold acclimation pathway; its activity is regulated by the clock component, TOC1, and PHYB (Kidokoro et al., 2009). PIF7 shows temperature-dependent AS altering the relative levels of the fully spliced, protein-coding transcript and a non-protein-coding transcript with retention of intron 1. PHYB is a photoreceptor required for photomorphogenesis and is involved in the interaction of light- and cold-signaling pathways (Lee and Thomashow, 2012; Wang et al., 2016). The two PHYB transcript isoforms differ by removal or retention of intron 3, which alters the C-terminal sequences of the predicted proteins by the presence/absence of the histidine kinase related domain required for full functionality (Krall and Reed, 2000; Müller et al., 2009). Finally, SUF4 is required for delayed flowering in winter-annual Arabidopsis and is involved in chromatin modification of FLOWERING LOCUS C (FLC) (Ding et al., 2013). SUF4 AS transcripts differ by splicing/retention of intron 6 to alter the C-terminal sequences (Kim and Michaels, 2006), and here, we demonstrate rapid cold-induced changes in AS (Figure 3F). Thus, in addition to the previously observed cold-induced changes in gene level expression, we demonstrate that low temperature-dependent AS is an important further layer of regulation that controls the abundance of cold response genes and transcripts.
DE and DAS Genes Have Different Predicted Functions
The DE and DAS gene sets were largely different with an overlap of only 795 genes (Figure 2A). The most significantly enriched Gene Ontology (GO) terms for DE genes correlated with known physiological and molecular events in the cold response such as response to various stresses and ribosome production (Knight and Knight, 2012; Beine-Golovchuk et al., 2018) (Figure 2C; Supplemental Data Set 6). Hierarchical clustering of total gene expression levels of DE genes revealed transient, adaptive, and late expression profiles in response to cold and regulation between light and dark (Figure 4; Supplemental Figure 5 and Supplemental Data Set 7). Transcript expression profiles of individual DE genes showed similar responses (Supplemental Figure 6). Functional annotation of individual DE gene clusters was associated with response to cold, abscisic acid, and drought (cluster 6), decreased photosynthesis (cluster 10), increased ribosome production (cluster 3), and membrane and lipid biosynthesis (cluster 12) (Figure 4; Supplemental Figure 5 and Supplemental Data Set 7). Cluster 12 shows highly increased gene level expression in the first 3 h of cold treatment reflecting the reconfiguration of membranes in response to cold to maintain fluidity and protect against subsequent freeze damage (Knight and Knight, 2012).
Hierarchical Clustering and Heat Map of Arabidopsis DE Cold-Responsive Genes and Key GO Terms.
DE genes show segregation into 12 coexpressed modules. For simplicity, genes that do not fall into any cluster have been removed from the heat map (n = 492). Cluster 3 (n = 762) includes genes enriched for “structural constituent of ribosome” and is upregulated mostly in the fourth day at 4°C. Cluster 9 (n = 1140) includes genes enriched for “response to auxin” and is downregulated upon cold, whereas cluster 12 (n = 571) includes genes enriched for “plasma membrane” and is upregulated within the first 3 h of cold treatment. Full results of GO enrichment analysis of heat map DE clusters are shown in Supplemental Data Set 7. The z-score scale represents mean-subtracted regularized log-transformed TPMs. BP (red bars), biological process; CC (blue bars), cellular component; MF (green bars), molecular function. The colored bars above the heat map indicate whether samples were exposed to light (colored) or dark (black) in the 3 h before sampling.
For the DAS genes, the most enriched functional terms were related to mRNA splicing (Figure 2C; Supplemental Data Set 6). One hundred and sixty-six (7%) DAS genes were RNA binding proteins, spliceosomal proteins, or SFs. Hierarchical clustering of the DTU transcripts and expression profiles of individual DAS genes also showed transient, adaptive, and late expression patterns (Figure 5; Supplemental Figures 7 to 9). Functional annotation of the genes in the individual DTU clusters showed enrichment of terms involved in the plasma membrane and signal transduction (cluster 8, P < 0.001) as well as regulation of transcription (cluster 3, P < 0.001), RNA splicing (cluster 5, P < 0.0001), and chromatin binding (cluster 1, P < 0.0001) (Figure 5).
Heat Map of DTU Transcripts from DAS Genes.
DTU transcripts from DAS genes show segregation into 10 coexpressed clusters. For simplicity, transcripts that do not fall into any cluster have been removed from the heat map (n = 36). Clusters 1 and 2 show transcripts downregulated upon cold. Clusters 3, 5, and 10 show clear transient changes in AS isoform transcripts at different times during day 1 at 4°C, while cluster 4 (n = 326) shows late upregulation of transcripts on the fourth day at 4°C. Clear gain in rhythmic expression of AS transcripts upon cold is seen in cluster 7 (n = 258). Cluster 8 (n = 233) includes transcripts with increased expression within the first 3 h of cold treatment. The z-score scale represents mean-subtracted regularized log-transformed TPMs. The colored bars above the heat map indicate whether samples were exposed to light (colored) or dark (black) in the 3 h before sampling.
Rapid Cold-Induced Changes in AS Accompany the Major Transcriptional Responses
The high temporal resolution of the time course allowed us to examine the relative dynamics of the DE and DAS changes. The statistical model used in this analysis allowed us to determine precisely at which time point each DE and DAS gene first showed a significant change (start time point), along with the magnitude and duration of that change. The dynamics of the changes of DE and DAS genes was compared by plotting the distribution of start time points (Figures 6A and 6B). DE and DAS genes peaked at 6 to 9 h after onset of cold (Figure 1A, T11 and T12) and 6 h after onset of cold (Figure 1A, T11), respectively. Furthermore, 62.2% (1520) of the DAS genes and 47.6% (3473) of the DE genes were detected within the first 9 h of low temperature (Figures 6A and 6B; Supplemental Data Set 8). The speed of induction is highlighted by 648 and 388 genes showing significant DE and DAS after only 3 h of cold (T10) (Figures 6A and 6B), respectively, despite the gradual reduction in temperature, which takes 2 h to reach 4°C (Supplemental Figure 4). Notably, three-quarters (76.5%; 849) of the DAS genes detected within 9 h of cold (Figure 1A, T10-T12) also had large AS changes with at least one transcript having ΔPS >0.2 (Supplemental Data Set 9). A further indicator of the speed of AS response and the sensitivity of some AS to reductions in temperature was demonstrated by identifying those DAS genes that show isoform switches (ISs), where the relative abundance of different isoforms is reversed in response to cold (Supplemental Figures 8 and 9). We developed the Time-Series Isoform Switch (TSIS) program to identify ISs in time-series data (Guo et al., 2017). Using TSIS, a total of 892 significant (P < 0.001) ISs that involved two abundant transcript isoforms (see Methods) were identified in 475 unique DAS genes (Figure 7A; Supplemental Data Set 10). The ISs involved different types of AS events and 77% involved at least one protein-coding transcript isoform (Figures 7A and 7B; Supplemental Data Set 10). These either generated isoforms that encoded different protein variants, or where AS occurred in the 5′ or 3′UTR, the transcript isoforms encoded the same protein, or one of the transcripts was non-protein-coding (e.g., PTC-containing). TSIS determines the two time points between which a significant isoform switch occurs and, consistent with the rapid changes in AS, the majority (57%) occurred between 0 and 6 h following transfer to 4°C (Figure 7A; Supplemental Data Set 10). Thus, immediately in response to lowering temperature, there are waves of transcriptional and AS activity involving thousands of genes.
Rapid Changes in DE and DAS Genes in Response to Cold.
Histograms and density plots of the time points at which the 7302 DE (A) and 2442 DAS (B) genes first become significantly different in day 1 and day 4 at 4°C compared with 20°C. The genes that first show significant differences after longer exposure to cold (day 4 at 4°C) represent 20.45% of DE and 14.73% of DAS genes. Each gene is represented only once in each histogram (left y axis). The estimated density line of the number of genes illustrates the early waves of transcriptional and alternative splicing responses (right y axis).
Sensitivity of AS to Low Temperatures.
(A) Frequency over time of isoform switches in the RNA-seq time course. Each isoform switch involved “abundant” transcripts (i.e., expression of each transcript makes up at least 20% of the total expression of the gene in at least one time point). Proportion of protein-coding transcripts is also shown and represents either production of different protein-coding isoforms or transcripts encoding the same protein where key AS events are in the UTR region. Data between T17 and T18 represent ISs that occurred between day 1 and day 4.
(B) Proportion of the major types of AS events involved in isoform switches in (A) was measured with SUPPA (Alamancos et al., 2015).
(C) Experimental design for assessing long-term changes in AS induced by small reductions in temperature initiated at dusk. Sampling of 5-week-old Arabidopsis rosettes occurred at dawn, after 12 h of temperature reduction, and is marked by a vertical arrow.
(D) AS of novel cold response gene GEMIN2 (PTC-containing transcript AT1G54380_c3) is sensitive to reductions in temperature of 2°C.
(E) Exon 9 skipping (E9.Skip) of novel cold response gene AT1G15200 (transcripts AT1G15200.1, AT1G15200_JS1, and AT1G15200_JS2) is sensitive to reductions in temperature of 8°C.
(F) Experimental design for assessing immediate changes in AS induced by gradual reductions in temperature initiated at dusk. Sampling of rosettes of 5-week-old Arabidopsis plants occurred at the indicated time points after dusk and is marked by vertical arrows.
(G) AS of novel cold response gene AT5G67540 (transcripts AT5G67540_P1, AT5G67540_P2, and AT5G67540_s1) is affected within 1 h of gradual reduction in temperature.
(H) AS of novel cold response gene LIGHT-RESPONSE BTB2 (LRB2/POB1; transcripts AT3G61600_P1 and AT3G61600_P2) is affected within 1 h of gradual reduction in temperature.
In (D), (E), (G), and (H), Tukey t tests were performed to compare each temperature reduction result against the 20°C control. Significant differences are labeled with asterisks (*P < 0.05, **P < 0.01, and ***P < 0.001). Error bars indicate sd, n = 3 biological replicates.
Cold-Regulated Expression and AS of Transcription and Splicing Factors
The waves of differential expression and AS in response to cold likely reflect regulation by TFs and SFs/RNA binding proteins (SF-RBPs). Differential expression of TFs in response to cold has been well documented (see Introduction). Here, 532 of the 2534 TFs predicted in Arabidopsis (Supplemental Data Set 11) were significantly regulated only at the gene level (DE-only) (Table 1). However, one-third (271) of TFs with cold-induced expression changes involved AS (Table 1), which potentially affects the levels or function of TF proteins. Similarly, of the 798 predicted SF-RBP genes (see Methods; Supplemental Data Set 11), 197 were DE-only, 33 DE+DAS, and 133 DAS-only (Table 1). Thus, many TF and SF-RBP genes were regulated by AS in response to lower temperatures. The majority of them have not previously been associated with the cold response and represent putative novel cold response factors (Supplemental Data Set 12). We next identified the TF and SF-RBP genes with the fastest (0–6 h after onset of cold) and largest changes in expression and AS (log2 fold change ≥1.5 [equivalent to 2.83-fold change] for DE genes and ΔPS >0.25 for at least one transcript in DAS genes) (Supplemental Data Set 12). For DE genes, 117 TFs and 11 SF-RBPs showed rapid and large expression changes while for the DAS genes 59 TF and 47 SF-RBP genes were identified as “early AS” genes. The TF genes included a high proportion of circadian clock genes as well as genes associated with abiotic stress, flowering time, and hormone responses (Supplemental Data Sets 12 and 13). The SF-RBP genes included serine-arginine-rich (SR) and heterogeneous ribonucleoprotein particle (hnRNP) protein genes such as SR30, RBM25, and GEMIN2 known to be involved in stress responses and/or regulation of the clock (Reddy et al., 2013; Schlaen et al., 2015; Cheng et al., 2017; Hartmann et al., 2018) (see Discussion). For many of the early AS genes, the AS changes were either only observed in day 1 at 4°C or persisted through the cold treatment, and many involved isoform switches (Supplemental Data Set 12). Thus, both transcription and AS determine expression levels of TFs and SFs and many of these genes were regulated rapidly in response to reduction in temperature.
Speed and Sensitivity of AS Responses to Small Reductions in Temperature
We previously noted changes in AS with a 4°C reduction in temperature (Streitner et al., 2013). The rapid and large changes in AS seen here suggest that many AS events are sensitive to relatively small changes in temperature. To investigate further, we examined the effect on AS of six early AS genes with large changes in AS (GEMIN2, LRB2/POB1, RDM16, AT1G06960 [U2B”-LIKE], AT1G15200, and AT1G15270). When the temperature was lowered from 20°C to 18°C, 16°C, 12°C, 8°C, and 4°C for 12 h, three genes showed significant changes in the level of at least one AS isoform with only a 2°C reduction in temperature (to 18°C), while the others were affected by 4°C or 8°C reductions (Figures 7C to 7E; Supplemental Figure 10). Thus, relative levels of AS transcript isoforms were dependent on the temperature. We then examined the speed and sensitivity of AS by taking multiple samples between 0 and 3 h after the onset of cold (Figures 7F to 7H; Supplemental Figure 11). Of eight genes examined, one showed significant AS changes within 40 min when the temperature had reached 11°C, five within 60 min (at 8°C), and two within 180 min (at 4°C) (Figures 7G and 7H; Supplemental Figure 11).
U2B”-LIKE Is Regulated at the AS Level and Is Required for Freezing Tolerance
Many of the early AS genes, including TFs and SFs, showed large and rapid changes in AS that alter the levels of protein-coding transcripts (Supplemental Data Set 13). We hypothesized that such significant changes in AS in response to low temperature are important in the overall cold acclimation process of the plant and lead to improved ability to tolerate freezing conditions after acclimation. In support of this, four of the early AS genes have previously been shown to be required for cold acclimation and tolerance to freezing: RCF1 and STA1 (Guan et al., 2013), GEMIN2 (Schlaen et al., 2015), and the LAMMER kinase, AME3 (L. Savitch, unpublished data) (Table 2). To examine whether other early AS genes are involved in cold acclimation, we selected the SF-RBP gene AT1G06960 because it was a novel DAS-only gene with an adaptive expression pattern and appeared to lose rhythmicity of the main transcript in the cold (we refer to AT1G06960 as U2B”-LIKE because of its similarity to U2B” [AT2G30260]; see below). We isolated a knockout mutant of the U2B”-LIKE gene (AT1G06960; Figure 8A; Supplemental Figure 12). U2B”-LIKE has two main AS transcripts, the fully spliced protein-coding mRNA and an isoform with retention of intron 4 (I4R; P1 and P2, respectively; Figure 8A). In wild-type plants, the protein-coding P1 transcript isoform showed rhythmic expression at 20°C and loss of rhythm during day 1 at 4°C, maintaining a high level of expression throughout the remaining cold treatment (Figure 8A). In freezing tolerance tests conducted at −8.0°C and −8.5°C, the u2b”-like mutant plants showed greater sensitivity to freezing; u2b”-like did not survive freezing at −8.5°C after cold acclimation, while wild-type plants recovered (Figure 8B). u2b”-like mutant and wild-type plants both recovered at −8.0°C. Differential sensitivity of the mutant was confirmed in quantitative electrolyte leakage analyses; leaf tissue of u2b”-like suffered significantly increased ion leakage (cellular damage) at −10°C than wild-type plants (Figure 8C), indicating that expression of U2B”-LIKE is required for cold acclimation and freezing tolerance.
AT1G06960 (U2B”-LIKE) Expression Profile and Assays of the Knockout Line.
(A) Structures of highly expressed U2B”-LIKE transcripts (black boxes, UTR; color-coded boxes, exon coding sequence) and gene/transcript expression profile across the time course. I4R, Intron 4 Retention; Alt5′ss, alternative 5′ splice site. Black/white bars below expression plots represent 12-h dark/light cycles.
(B) Freezing sensitivity of cold-acclimated Col-0 and u2b”-like mutants showing recovery of wild-type and nonrecovery of u2b”-like mutant plants at −8.5°C.
(C) Cellular ion leakage in Col-0 (wild-type) and u2b”-like (knock-out mutant) leaf discs subjected to different freezing temperatures before thawing (n = 4). Transformed ion leakage data were used in a one-tailed t test, which confirmed u2b”-like loses more electrolytes than wild-type Col-0 at −10°C (*P = 0.0263). Each bar of the plot represents average ion leakage values. In (B) and (C), plants were grown at 20°C for 4 to 5 weeks and cold acclimated at 5°C for 2 weeks before the freezing assay.
(D) and (E) Relevant regions of gene structures (left) and HR RT-PCR data (right) of plants harvested at dawn of 4°C day 1 are shown (n = 3). AS isoform levels obtained with HR RT-PCR analysis are shown relative to the fully spliced (FS) isoform (AS/FS ratio). The data were compared using a two-way ANOVA.
(D) Relative levels of PTC-containing intron 1 of the PIF7 (AT5G61270) gene is increased in u2b”-like mutant transcripts when compared with Col-0 (**P = 0.002).
(E) Transcription factor HY5 homolog (HYH; AT3G17609). The AS ratio of the PTC-containing I1R transcript shows a highly significant increase in the u2b’’-like mutant compared with Col-0 (***P = 0.001). Significant changes in the AS ratio of the above genes in the u2b’’-like mutant is evidence that the U2B’’-LIKE might be involved in regulating their AS. Black arrows represent the location on the gene of the forward and reverse primers used for HR RT-PCR. IR, intron retention. More data are in Supplemental Data Set 14.
In (A) and (C), error bars are se of the mean (n = 3 biological replicates × 3 sequencing replicates in [A], and n = 4 biological replicates × 6 pseudo-replicates in [C]), whereas in (D) and (E), they represent sd (n = 3 biological replicates).
Arabidopsis contains two U2B”-related genes: U2B” (AT2G30260) (Simpson et al., 1995) and U2B’’-LIKE (AT1G06960). The two proteins are very similar: 80% identical and 90% similar at the amino acid sequence level (Supplemental Figure 13). U2B” is an U2snRNP-specific protein that binds, along with U2A’, to stem-loop IV of U2snRNA in both plants and human (Simpson et al., 1995). In the u2b”-like mutant, there was no expression of U2B”-LIKE but expression of the U2B” paralog (which was neither DE nor DAS in cold) was detected (Supplemental Figure 12), suggesting that U2B” protein could not compensate for the lack of U2B”-LIKE in the u2B”-like mutant and, therefore, that they had functionally diverged. To investigate whether U2B”-LIKE affected AS regulation, we then compared AS patterns of 41 genes (including 34 DAS or DE+DAS genes identified here) in wild-type and u2b”-like mutant plants. Five genes showed significantly different AS (P < 0.05 and >10% difference in splicing ratio between the mutant and the wild type) (Figures 8D and 8E; Supplemental Figure 14 and Supplemental Data Set 14). These included decreased levels of fully spliced, protein-coding transcripts of PIF7 (AT5G61270; Figure 8D), which along with TOC1 and PHYB represses expression of CBFs (Lee and Thomashow, 2012), and HOMOLOG OF HY5 (HYH) (Figure 8E), a clock input gene that is also involved in regulating anthocyanin synthesis at low temperatures. Thus, U2B”-LIKE is one splicing factor that contributes to correct splicing of PIF7, linking U2B”-LIKE-dependent AS to regulation of the major cold response pathway. Therefore, the freezing sensitivity of the u2b”-like mutant may be due to altered AS and expression of specific genes required for cold acclimation.
DISCUSSION
Dynamic changes in expression at both the gene and transcript/AS levels occur in Arabidopsis plants in the process of cold acclimation. The extensive AS information identified here demonstrates a much higher degree of complexity of regulation in response to cold that has been significantly underestimated by analysis of differential gene expression only. In particular, we demonstrate the dynamic contribution of AS by the rapid cold-induced wave of AS activity accompanying the transcriptional response (Figure 6) and the sensitivity of AS of some genes to small reductions in temperature (Figure 7). We also significantly demonstrate the extent of AS by showing that over 2400 genes are regulated by AS in response to cold with over 1600 regulated only at the AS level (Figure 2). Notably, over 1300 of the AS genes were not differentially expressed at the gene level and would not be identified in microarray or gene level RNA-seq analyses and are novel cold responsive genes. The massive changes in expression and AS involved thousands of genes reflecting activation of both transcription and splicing pathways and networks. The speed and extent of the cold-induced AS suggest that AS, along with the transcriptional response, is a major driver of transcriptome reprogramming for cold acclimation and freezing tolerance.
With over 2400 genes regulated by AS, multiple different mechanisms are likely to control the splicing decisions. Reduction in temperature to 4°C is expected to reduce the rate of biochemical reactions and potentially affect transcription and splicing. We observed that the vast majority of introns in the pre-mRNAs of all the cold-expressed genes are efficiently spliced throughout the cold treatment. Therefore, low temperature does not cause a general defect in splicing reflecting the ability of temperate plants to grow in a wide range of fluctuating temperatures. Nevertheless, low temperatures may directly affect AS regulation. For example, in mammals, secondary structures in pre-mRNAs affect splice site selection (Hiller et al., 2007) and cooling could stabilize such structures. Similarly, splicing is largely cotranscriptional and slower rates of RNA polymerase II (Pol II) elongation promote selection of alternative splice sites (Luco et al., 2011). Both of these mechanisms will undoubtedly be involved in the cold-induced AS changes of some of the genes seen here. However, the sensitivity of AS to reductions in temperature of only a few degrees and clear rhythmic expression profiles of AS transcript isoforms in plants exposed to constant 4°C temperature for 4 d (e.g., cluster 7 in Figure 5 and Supplemental Figure 7) argue against such mechanisms being widely responsible for the cold-induced AS changes observed here. Local or global DNA methylation and chromatin modifications can also affect the rate of Pol II elongation or help to recruit SFs to affect splice site selection (Luco et al., 2011). In plants, epigenetic regulation is responsible for suppression of FLC by vernalization in the seasonal response to cold (Berry and Dean, 2015). Furthermore, altered histone 3 lysine 36 trimethylation (H3K36me3) was recently shown to affect some AS events induced by higher ambient temperatures within 24 h (Pajoro et al., 2017). Alongside dynamic changes in histone marks at specific stress-induced genes (Kim et al., 2012; Haak et al., 2017), it is likely that some of the cold-induced AS here reflects local epigenetic changes.
We showed that the levels of hundreds of TF and SF-RBP gene transcripts changed in response to cold at both the transcriptional and AS levels. Therefore, splicing decisions in the physiological response to low temperature are most likely controlled by altered abundance, activity, or localization of SFs or other RNA-interacting proteins (Reddy et al., 2013; Staiger and Brown, 2013; Fu and Ares, 2014; Lee and Rio, 2015; Verhage et al., 2017). In particular, we identified TF and SF-RBP genes with large and rapid changes in AS. Most of the early AS transcription factor genes were regulated only by AS and therefore had not been identified previously as cold response transcription factors. Nevertheless, the rapid cold-induced changes in the AS of some known cold response genes: CAMTA3, which activates the CBFs (Doherty et al., 2009), and VRN2 and SUF4, which are involved in vernalization and silencing of FLC (Ding et al., 2013; Berry and Dean, 2015), have not been described previously, and our results introduce AS as a novel component in their regulation. It will be interesting to address the function of the novel AS-regulated TFs and the function of AS of these and known cold response TFs on cold acclimation and vernalization in future experiments.
By contrast, the early AS SF-RBP genes included SR and hnRNP protein genes known to respond to changes in temperature (e.g., SR30, RS40, GRP8, SR45A, PTB1, and RBM25) (Reddy et al., 2013; Staiger and Brown, 2013; Cheng et al., 2017). Many SF-RBP genes are regulated by AS-NMD, and the rapid induction of AS in these and other early AS genes affects the abundance of protein-coding transcripts and presumably of the splicing factors themselves to alter AS of downstream targets. Various spliceosomal and snRNP protein genes are also among the early AS genes. These include GEMIN2 (snRNP assembly), which is cold-induced, is involved in regulation of the circadian clock, and enhances U1snRNP assembly to compensate for reduced functionality of U1snRNP at low temperatures (Schlaen et al., 2015). Interestingly, a number of U1snRNP core and associated protein genes (U1-70k, LUC7B, LUC7RL, PRP39A, and RBM25) (Barta et al., 2012; Amorim et al., 2017; Kanno et al., 2017) respond rapidly to cold via AS. The early AS genes also include two wound-induced RNA binding proteins, UBA2a and UBA2c (Bove et al., 2008), and may also therefore be involved in the cold response. Three LAMMER kinase genes (AFC1, AFC2, and AME3) that regulate SR proteins via phosphorylation showed changes in their expression due to AS, suggesting that lower temperatures affect activation/deactivation of specific splicing factors that are targets of these kinases (Savaldi-Goldstein et al., 2003) (L. Savitch, unpublished data). In addition, over 20 putative RNA binding proteins, kinases and RNA helicases with little or no known functional information are among the novel early AS genes. Four of the early AS SF-RBP genes (RCF1, STA1, GEMIN2, and AME3; Table 2) have been shown to be involved in freezing tolerance (Schlaen et al., 2015) (L. Savitch, unpublished data). Our results identify over 100 splicing and transcription regulatory genes, whose expression is rapidly and drastically altered by AS in response to cooling. Future work will address the function of these putative regulators and specific transcript isoforms in cold acclimation.
In addition to the early AS genes required for freezing tolerance (above), we provide initial evidence for another early AS gene, U2B”-LIKE, being involved in freezing tolerance and acclimation (Figure 8). The freezing sensitivity and AS effects were surprising as Arabidopsis contains two paralogs: U2B”-LIKE (AT1G06960) and U2B” (AT2G30260). Animals and plants generally have two homologous proteins termed U1A and U2B” that bind to similar stem-loop sequences U1snRNA and U2snRNA, respectively (Simpson et al., 1995; Price et al., 1998; Delaney et al., 2014). The N-terminal RRM1 is responsible for UsnRNA binding, and there is no known function for the conserved C-terminal RRM2. Human contains a single copy of U1A and U2B”, which arose through duplication and specialization (Delaney et al., 2014), while Drosophila melanogaster contains only a single gene (SANS-FILLE [SNF]) whose protein binds to both U1 and U2snRNA. Arabidopsis and some other plants are unique in having and maintaining two copies of U2B” genes and proteins. The Arabidopsis proteins are 80% identical (90% similar) but protein sequence alignments of the U2B” and U2B”-LIKE gene paralogs from Brassicaceae species show clear distinctions in a tyrosine (U2B”) to valine (U2B”-LIKE) substitution in the conserved RNP1 motif of RRM1 and the C-terminal two-thirds of the linker region (Supplemental Figure 13). The Tyr-Val change in RNP1 and the intradomain linker region may affect the RNA binding properties of U2B”-LIKE; the linker region in human U1A and Drosophila SNF affects their RNA binding properties (Williams and Hall, 2010). Mutations in U2B” genes in other organisms have previously been shown to have a general splicing defect. The U2B”-LIKE mutant did not have a general splicing defect but affected the AS of specific genes. Both human U1A and Drosophila SNF are involved in specific AS regulation (Boelens et al., 1993; Cline et al., 1999; Hu et al., 2009); recently, Arabidopsis U1A was also shown to regulate AS (Gu et al., 2018). Thus, U1A/U2B”-related proteins have other functions besides snRNA binding, including binding to pre-mRNAs and the regulation of AS, and Arabidopsis U2B”-LIKE has novel AS regulation functions.
The speed of change of AS may be one of the earliest responses to cooling. We showed significant AS within only 40 to 60 min of cooling and with subtle reductions in temperature of as little as 2°C (Figure 7). Similar responses are seen in mammals where neuronal stimulation and rapid changes of intracellular sterols also activate splicing/AS within minutes without de novo transcription or protein production (Medina et al., 2011; Mauger et al., 2016). In addition, a 1°C change in body temperature activated a program of AS changes within 30 min that involved temperature-sensitive phosphorylation of SR proteins and AS of the U2-associated factor, U2AF26 (Preußner et al., 2017). Therefore, cold-induced AS programs may similarly involve rapid phosphorylation/dephosphorylation of SFs, such as SR proteins, to modulate AS (Stamm, 2008; Preußner et al., 2017). Interestingly, the expression of a third of the early AS SF-RBP genes revealed here is also affected by increased ambient temperatures (Jung et al., 2016), and these genes may represent key temperature-dependent splicing regulators. Master splicing regulators have been postulated to drive splicing networks during cellular differentiation in mammals where regulatory modules of SF-RBPs and/or TFs establish and maintain new expression patterns (Jangi and Sharp, 2014; Fiszbein and Kornblihtt, 2017). Auto- and cross-regulatory modules of some plant SFs are well documented and may be important components of splicing networks (Reddy et al., 2013; Staiger and Brown, 2013).
The cold response pathway in Arabidopsis involves Ca2+-dependent and MAP kinase signaling cascades that affect inducer of CBF expression 1 (ICE1) phosphorylation, lead to activation of CBFs and other TFs, and expression of COR genes (Teige et al., 2004; Knight and Knight, 2012; Li et al., 2017; Zhao et al., 2017). In animals, signaling pathways including calcium-dependent signaling control of phosphorylation and AS of specific SFs regulate AS of downstream targets (Stamm, 2008; Razanau and Xie, 2013). In plants, stress signals affect both the phosphorylation status and subcellular localization of some SR proteins (de la Fuente van Bentem et al., 2006; Rausin et al., 2010), and 16 of the early AS SF-RBPs have been shown to be phosphorylated in Arabidopsis cells (de la Fuente van Bentem et al., 2006). The cold-induced waves of transcription and AS, the rapid AS responses of SF-RBPs, and their potential for phosphorylation suggest a model where cold signaling pathways modulate both transcription and splicing factor levels and activity (Figure 9). These regulatory factors, in turn, drive gene and splicing networks required to determine the overall reprogramming of the transcriptome for cold acclimation and freezing tolerance. These networks are reflected in dynamic changes in the DE and DAS gene sets that we observe across the time series.
Model for the Cold Signaling Pathway and Its Regulation of Genome-Wide Gene Expression.
Cooling activates Ca2+-dependent kinases and MAP kinases that activate or repress transcription factors or splicing factors. These in turn regulate the transcription or AS of downstream genes including other TFs and SFs, thereby driving cascades of cold-induced gene expression.
Plants are exposed to a variety of temperature conditions. They require flexible regulatory systems that modify expression quickly and reversibly upon perception of constantly fluctuating temperatures throughout the day and night and during the regular 24-h cycle of warmer daytime and cooler nighttime temperatures. They must also reprogram the transcriptome in cold conditions to allow the plant to acclimate, avoid freezing, and survive as the intensity and duration of reduced temperatures increase (seasonal changes). The dynamic AS response, sensitivity of AS to small changes in temperature, and the different behavior of AS genes, where changes are transient or persist, demonstrated here, suggest that AS provides a level of flexibility to contribute to different stages of the progression from perception to acclimation. In particular, the speed of AS reactions may contribute to temperature perception by altering the activity of key TFs and SF-RBPs, while the transcriptional response is being activated as seen in animal cells (Mauger et al., 2016). Such control could fine-tune expression of specific genes and pathways throughout the day as temperatures fluctuate. The rapid waves of transcriptional and AS activity within the first few hours of cold exposure, which include transcriptional activation of CBFs, other transcriptional response pathways, and altered expression/AS of clock components may be involved in initial cold responses and in the normal 24-h day/night temperature cycle. Significant changes in expression/AS of many genes occur rapidly in the first day in the cold and persist throughout the cold period and other genes are only activated or repressed by transcription or AS after prolonged cold treatment (day 4); these genes may be important for establishing and stabilizing changes in the transcriptome for acclimation. Thus, temperature-dependent AS is a mechanism to transduce temperature change signals into changes in expression. Dynamic and extensive changes in AS are also likely to drive plant responses to other abiotic stresses, pests, and diseases, and in developmental programs alongside transcriptional responses. Construction of splicing and transcriptional networks from the data here will further define the contribution of AS, as an additional layer of regulation, and the interplay and coordination of the transcriptional and AS responses.
METHODS
Plant Material and Growth Conditions
Arabidopsis thaliana Col-0 seeds were surface sterilized, stratified in the dark at 4°C for 4 d, and grown hydroponically to maturity (5 weeks) in Microclima environment-controlled cabinets (Snijders Scientific), maintaining 20°C, 55% relative humidity, and 12 h light (150 μE m−2 s−1)/12 h dark as described previously (James et al., 2012). Illumination was provided by a mix of the following of cool white tubes: 3 Philips Master TL-D 58W/840, 3 Sylvania BRITEGRO F58W/T8/2023, 16 Osram Lumilux FQ 54W/840, and 16 Sylvania LuxLine Plus FHO 24W/T5/84010-13. Arabidopsis rosettes were harvested and pooled at each sampling time. Harvesting occurred every 3 h beginning with the last 24 h at 20°C, and on days 1 and 4 after transfer to 4°C giving 26 time points in the time series (Figure 1A). Day 1 at 4°C represents the “transition” from 20°C to 4°C when plants first begin to experience the temperature decrease; day 4 at 4°C represents “acclimation” where plants have been exposed to 4°C for 4 d (Figure 1A). Three biological replicates were generated for each time point in separate experiments (78 samples in total). The same growth cabinet was used for all repeats to eliminate the potential effects of minor changes in light intensities and light quality on gene expression. Additionally, to avoid interference in the experiment from possible microclimates within the growth cabinet, 26 trays for each time point were placed in a randomized fashion. The switch to 4°C from 20°C was initiated at dusk. In a temperature reduction, the cabinet used here typically takes 2 h to reach 4°C air temperature (Supplemental Figure 4). Tissue was rapidly frozen in liquid N2 and stored at −80°C until isolation of RNA and preparation of cDNA.
RNA Extraction
Total RNA was extracted from Arabidopsis tissue using an RNeasy Plant Mini Kit (Qiagen), followed by either on-column DNase treatment (for high-resolution RT-PCR; see below), or the Turbo DNA-free kit (Ambion) (for library preparation and RT-qPCR; see below).
Library Preparation and Sequencing
RNA-seq libraries were constructed for 78 RNA samples by following instructions for a TruSeq RNA library preparation (Illumina protocol 15026495 rev. B). In these preparations, poly(A) selection was used to enrich for mRNA, RNA was fragmented for 8 min at 94°C, and random hexamers were used for first-strand cDNA synthesis. The 78 libraries had an average insert size of ∼280 bp, and each library was sequenced on three different sequencing lanes (27 lanes were used in total) of Illumina HiSeq 2500 platform generating 100-bp paired-end reads. The total number of raw reads generated in the RNA-seq data was 9.52 Bn paired-end reads giving ∼360 M paired-end reads per time point (Supplemental Data Set 15).
Residual adaptor sequences at both 5′ and 3′ ends were removed from raw reads using cutadapt version 1.4.2 (https://pypi.org/project/cutadapt/1.4.2) with quality score threshold set at 20 and minimum length of the trimmed read kept at 20. The “–paired-output” option was turned on to keep the two paired read files synchronized and avoid unpaired reads. The sequencing files before and after the trimming were examined using fastQC version 0.10.0.
Quantification of Transcripts and AS
Arabidopsis transcript expression from our RNA-seq experiment was performed using Salmon version 0.82 (Patro et al., 2017) in conjunction with AtRTD2-QUASI augmented by eight genes that were not originally present (Zhang et al., 2017). For indexing, we used the quasi mapping mode to build an auxiliary k-mer hash over k-mers of length 31 (–type quasi –k 31). For quantification, the option to correct for the sequence specific bias (“–seqBias”) was turned on. The number of bootstraps was set to 30 and all other parameters were on default settings. AtRTD2-QUASI is a modification of AtRTD2, a high-quality reference transcript data set for Arabidopsis Col-0 containing >82,000 unique transcripts, designed for the quantification of transcript expression (Zhang et al., 2017). Their use in validation of transcript structures and accurate quantification of individual transcript abundances for alternative splicing analyses was demonstrated previously using the biological repeats of two of the time-points analyzed here (Zhang et al., 2017). Transcript expression results are in Supplemental Data Set 16.
DE Gene and DAS Analysis of the RNA-Seq Data
To carry out differential expression analysis, transcript quantification results generated by Salmon were processed and refined in successive steps. First, transcript and gene read counts were generated from TPM data correcting for possible gene length variations across samples using tximport version 0.99.2 R package with the option “lengthScaledTPM” (Soneson et al., 2015). Second, read count data from sequencing replicates were summed for each biological sample. Third, genes and transcripts that are expressed at very low levels were removed from downstream analysis. The definition of a low expressed gene and transcript was determined by analyzing mean-variance relationships (Law et al., 2016). The expected decreasing trend between the means and variances was observed in our data when removing transcripts that did not have ≥1 counts per million in three or more samples out of 78, which provided an optimal filtering for low expression transcripts. At the gene level, if any transcript passed the expression level filtering step, the gene was included as an expressed gene and then the normalization factor, which accounted for the raw library size, was estimated using the weighted trimmed mean of M values method using edgeR version 3.12.1 (Robinson et al., 2010). Principal component analysis showed significant batch effects within the three biological replicates. Thus, batch effects between biological repeats were estimated using RUVSeq R package version 1.4.0 with the residual RUVr approach (Risso et al., 2014). Normalized read counts in counts per million were then log2 transformed and mean-variance trends were estimated and weights of variance adjustments were generated using the voom function in limma version 3.26.9 (Law et al., 2014, 2016; Ritchie et al., 2015).
General linear models to determine differential expression at both gene and transcript levels were established using time and biological replicates as factors and 18 contrast groups were set up where corresponding time-points in the day 1 and day 4 at 4°C blocks were compared with those of the 20°C block (e.g., block2.T1 versus block1.T1, block2.T2 versus block1.T2, etc.; Supplemental Figures 1 and 2). To detect DE genes, the log2 fold changes of gene abundances were calculated in each contrast group and significance of the changes in expression were determined using t tests. P values were adjusted by the BH method for multiple testing correction (Benjamini and Hochberg, 1995). For DAS analysis, the log2 fold changes (L2FCs) of each individual transcript were compared with the gene level L2FC, which is the weighted average of the L2FC of all the transcripts of the gene based on the size of their standard deviations. Significance of the DAS was obtained by testing the consistency of expression changes among all the transcripts and the changes of expression at the gene level using an overall F-test for the same 18 contrasts using the DiffSplice function (Supplemental Figure 2) (Ritchie et al., 2015). For DTU analysis, L2FCs of each transcript from DAS genes were compared with the L2FC of the weighted average of all the other transcripts from the same gene. Significance of the DTU was obtained by testing the consistency of expression changes between each transcript and the change of expression of the other transcripts was tested using a t test for the same 18 contrasts using the DiffSplice function (Supplemental Figure 2). DTU transcripts which coded for protein isoforms were determined from transcript translations of AtRTD2 (Zhang et al., 2017).
Genes were significantly DE at the gene level if they had at least two contrast groups at consecutive time points with adjusted P < 0.01 and ≥2-fold change in expression in each contrast group (Supplemental Figure 1). Genes/transcripts with significant DAS/DTU had at least two consecutive contrast groups with adjusted P < 0.01 and with these contrast groups having at least one transcript with ≥10% change in expression (Supplemental Figure 2). Gene functional annotation was performed using R package RDAVIDWebService version 1.8.0 (Huang et al., 2009a, 2009b; Fresno and Fernández, 2013). The possibility of a gene and transcript being identified by accident as cold responsive by the statistical method was tested. Time points T1 and T9 (Figure 1A) are virtually identical as they both represent dusk samples at 20°C, the only difference being they are 24 h apart, such that few DE or DAS genes and DTU transcripts were expected when comparing these time points. Indeed, no significant DE or DAS gene, nor DTU transcript, was identified between T1 and T9. This suggests our statistical method to select cold-responsive genes and transcripts is conservative and controls the number of false positives.
Identification of Isoform Switches
For the isoform switch analysis, we used the TSIS R package, which is a tool to detect significant transcript isoform switches in time-series data (Guo et al., 2017). For the analysis, 5317 high abundance transcripts, whose average expression accounts for >20% of total gene expression at least one time point, were selected from the DAS gene transcripts. Switches between any two time points were identified using the default parameters in which (1) the probability of switch (i.e., the frequency of samples reversing their relative abundance at the switches) was set to >0.5; (2) the sum of the average differences of the two isoforms in both intervals before and after the switch point were set at ΔTPM >1; (3) the significance of the differences between the switched isoform abundances before and after the switch was set to P < 0.001; and (4) both intervals before and after switch must consist of at least 2 consecutive time points to detect long lasting switches. SUPPA version 2.1 (Alamancos et al., 2015) was then used to identify the specific AS events (e.g., intron retention, alternative 3′ or 5′ splice site selection, exon skip) that distinguished each pair of switched transcript isoforms.
Reference Gene Lists of RNA Binding, Splicing Factor, and Spliceosomal Protein and TF Genes
To identify putative splicing regulators/RNA binding proteins among the DE, DE+DAS, DAS, and TSIS genes, a reference list of Arabidopsis genes encoding RBPs, spliceosomal proteins, and SFs was assembled. First, published RRM- and KH-domain-containing RBPs from Lorković and Barta (2002) were combined with orthologs of human and yeast splicing proteins (Wang and Brendel, 2004; Barta et al., 2012) to give a nonredundant list of 317 SF-RBP genes. Recently, a protocol where proteins were cross-linked to RNA and then poly(A)+ RNA was isolated and bound proteins were identified by proteomics/mass spectrometry was applied to Arabidopsis (Marondedze et al., 2016; Reichel et al., 2016). Of these, Reichel et al. (2016) provides a comprehensive discussion of the possible RBPs and produced a high confidence list of 300 AtRBPs and a further 446 candidate RBPs. We combined the 317 genes with the AtRBP genes and removed redundancy to generate a high confidence list of 526 genes. In a second approach, all Arabidopsis genes were searched for GO terms “RNA localization,” “RNA catabolic process,” “RNA splicing,” “RNA splicing via transesterification reactions,” “mRNA processing,” “RNA helicase,” “RBP complex,” “SR and RBP proteins,” and “mRNA splicing factors” using PANTHER12.0 software (Mi et al., 2017). In addition, the TAIR Protein Domain database was manually curated for terms “RRM,” “RNA,” “RNA helicase,” “KH,” “ribonucleo,” and “splicing.” The resulting putative RNA binding/interacting proteins and the candidate AtRBPs from Reichel et al. (2016) were combined and duplicates, ribosomal proteins, translation factors, and chloroplast and mitochondrial proteins were removed. Finally, genes that have been shown to affect alternative splicing but that do not encode RNA binding proteins were added manually. This generated a reference list of 798 protein-coding genes with RNA interacting potential referred to as SF-RBP 798.
To generate a reference list of Arabidopsis transcription factors, the TF list from TAIR and Araport 11 (thale mine) of TF families was combined with a published list of TFs (Pruneda-Paz et al., 2014) and genes identified by orthology searches with DNA binding domains. Evidence of DNA binding was taken into account and redundancy removed to generate a list of 2534 TFs (Supplemental Data Set 11).
RT-qPCR
RT-qPCR was performed essentially as described previously (James et al., 2012). cDNA was synthesized from 2 μg of total RNA using oligo(dT) primers and SuperScriptII reverse transcriptase (Thermo Fisher Scientific). Each reaction (1:100 dilution of cDNA) was performed with Brilliant III SYBR Green QPCR Master Mix (Agilent) on a StepOnePlus real-time PCR system (Fisher Scientific UK). The average Ct values for PP2A (AT1G13320) and IPP2 (AT3G02780) were used as internal control expression levels. The delta-delta Ct algorithm (Livak and Schmittgen, 2001) was used to determine relative changes in gene expression. Primer sequences are provided in Supplemental Data Set 17A.
High-Resolution RT-PCR
High-resolution (HR) RT-PCR reactions were conducted as described previously (Simpson et al., 2008). Gene-specific primer pairs were used for analyzing the expression and alternative splicing of different genes (Supplemental Data Set 17). For each primer pair, the forward primer was labeled with 6-carboxyfluorescein. cDNA was synthesized from 4 μg of total RNA using the Sprint RT Complete Double PrePrimed Kit following the manufacturer’s instructions (Clontech Laboratories, Takara Bio). The PCR reaction usually contained 3 μL of diluted cDNA (1:10) as a template, 0.1 μL of each of the forward and reverse primers (100 mM), 2 μL of 10× PCR buffer, 0.2 μL of Taq polymerase (5 units/μL; Roche), 1 μL of 10 mM dNTPs (Invitrogen, Life Technologies), and RNase-free water (Qiagen) up to a final volume of 20 μL. For each reaction, an initial step at 94°C for 2 min was used followed by 24 to 26 cycles of (1) denaturation at 94°C for 15 s, (2) annealing at 50°C for 30 s, and (3) elongation at 70°C for either 1 min (for fragments smaller than 1000 bp) or 1.5 min (for fragments between 1000–1200 bp) and a final extension cycle of 10 min at 70°C.
To separate the RT-PCR products, 1.5 μL of PCR product was mixed with 8.5 μL of Hi-Di formamide (Applied Biosystems) and 0.01 μL of GeneScan 500 LIZ dye or 0.04 μL of GeneScan 1200 LIZ dye size standard and run on a 48-capillary ABI 3730 DNA Analyzer (Applied Biosystems, Life Technologies). PCR products were separated to single base-pair resolution, and the intensity of fluorescence was measured and used for quantification in relative fluorescent units. The different PCR products and their peak levels of expression were calculated using Genemapper software (Applied Biosystems, Life Technologies).
Identification and Characterization of the u2b’’-like Mutant
cDNA of SALK_060577 line was synthesized as described above for HR RT-PCR. PCR was performed using cDNAs and GoTaq Green DNA polymerase (Promega) following the manufacturer’s instructions. Primer sequences are provided in Supplemental Data Set 17.
Freezing and Electrolyte Leakage Assay
Cold-acclimated plants were assessed for damage after freezing conditions. Sterilized seeds were sown on MS-agar plates and after 7 d seedlings were transferred to peat plugs for growth in 12/12 light/dark cycles, 150 to 200 μE/m2/s (illumination was generated from a mix of three cool white F70W/840 CCT-4000K tubes to two Sylvania GRO-LUX F58W/GRO-T8 tubes) at 20°C for 4 weeks. Plants were then transferred to 5°C, 10/14 light/dark cycles (150 μE/m2/s illumination generated from fluorescent lamps FL40SS W/37) for ∼14 d, after which they were used in either a qualitative or quantitative assay. In the qualitative assay, cold-acclimated plants were transferred at dusk to either −8.0°C or −8.5°C for 24 h and then transferred to 5°C, 10/14 light/dark cycles (150 μE/m2/s illumination from fluorescent lamps FL40SS W/37) for 24 h and finally to 12/12 light/dark cycles, 150 to 200 μE/m2/s (illumination generated from a mix of three cool white F70W/840 CCT-4000K tubes to two Sylvania GRO-LUX F58W/GRO-T8 tubes) at 20°C for 1 week, after which they were assessed for signs of regrowth, indicating survival. In the quantitative assay, we performed the electrolyte leakage test (Hemsley et al., 2014). In brief, three leaf discs were collected from each cold-acclimated plant, forming a pseudo-replicate. For each temperature and genotype, six pseudo-replicates were obtained, representing one biological replicate. Ice nucleation was initiated in individual test tubes for each pseudo-replicate by introducing ice chips and tubes were cooled progressively to the sub-zero temperatures indicated. Conductivity measurements were made after thawing and then again after complete loss of all electrolytes, to give a percentage measurement of electrolyte loss in each sample. In total, four biological replicates were analyzed. Percentage ion leakage data were first divided by 100 and then square-root and arc-sine transformed before analysis in one-tailed t tests.
Accession Numbers
Sequence data from this article can be found in the GenBank/EMBL data libraries under accession number PRJEB19974.Major genes mentioned in this study are as follows: PRP39A (AT1G04080), U2B”-like (AT1G06960), SR45A (AT1G07350), SR30 (AT1G09140), AT1G15200, AT1G15270, COR47 (AT1G20440), RCF1 (AT1G20920), ZAT10 (AT1G27730), RDM16 (AT1G28060), ERF11 (AT1G28370), SUF4 (AT1G30970), GEMIN2 (AT1G54380), RBM25 (AT1G60200), PHYB (AT2G18790), CAMTA3 (AT2G22300), U2B” (AT2G30260), WRKY33 (AT2G38470), U1A (AT2G47580), PTB1 (AT3G01150), UBA2c (AT3G15010), HYH (AT3G17609), ICE1 (AT3G26744), U1-70K (AT3G50670), AFC1 (AT3G53570), UBA2a (AT3G56860), LRB2/POB1 (AT3G61600), STA1 (AT4G03430), VRN2 (AT4G16845), AFC2 (AT4G24740), CBF1 (AT4G25470), CBF2 (AT4G25480), CBF3 (AT4G25490), RS40 (AT4G25500), AME3 (AT4G32660), GRP8 (AT4G39260), FLC (AT5G10140), KIN1 (AT5G15960), LUC7B (AT5G17440), LUC7RL (AT5G51410), COR78 (AT5G52310), CRF3 (AT5G53290), ZAT12 (AT5G59820), PIF7 (AT5G61270), and TOC1 (AT5G61380). The Reference Transcript Data Sets AtRTD2 and AtRTD2-QUASI are available in the James Hutton Institute repository (http://ics.hutton.ac.uk/atRTD/). The codes of the DE/DAS/DTU analysis pipeline are available at the Github page https://github.com/wyguo/AtRTD2-DE-DAS-DTU-pipeline. The expression profiles of all Arabidopsis genes in the time-series data can be viewed interactively and downloaded through a web service at https://wyguo.shinyapps.io/atrtd2_profile_app/.
Supplemental Data
Supplemental Figure 1. Differential expression analysis.
Supplemental Figure 2. Differential alternative splicing and transcript usage analysis.
Supplemental Figure 3. Number of genes that are DE and DAS at each time point in day 1 and day 4 at 4°C when compared with the 20°C control.
Supplemental Figure 4. Graph of air temperature reduction over time inside Microclima growth cabinet (Snijders Scientific).
Supplemental Figure 5. Average expression profiles of DE gene clusters from heat map in Figure 4.
Supplemental Figure 6. Expression profiles of cold response DE-only genes.
Supplemental Figure 7. Average expression profiles of DTU transcript clusters from heat map in Figure 5.
Supplemental Figure 8. Expression profiles of cold response DE+DAS genes.
Supplemental Figure 9. Expression profiles of cold response DAS-only genes.
Supplemental Figure 10. Sensitivity of AS to small reductions in temperature.
Supplemental Figure 11. Rapid changes in AS in response to gradual decrease in temperature from 20°C to 4°C in the first 0 to 3 h of the cold treatment.
Supplemental Figure 12. Identification and characterization of u2b’’-like mutant.
Supplemental Figure 13. Alignment of protein sequences of U2B” and U2B”-LIKE homologs in six Brassicaceae species.
Supplemental Figure 14. AS events significantly affected in the u2b”-like knockout plants compared with Col-0 (wild type). The following materials have been deposited in the DRYAD repository under accession number http://dx.doi.org/10.5061/dryad.fk1cj47.
Supplemental Data Set 1. Results of the analysis of differentially expressed, differentially alternatively spliced, and differential transcript usage.
Supplemental Data Set 2. Gene lists of the up- and downregulated DE genes in the different contrast groups.
Supplemental Data Set 3. Summary of previous studies on wild-type Arabidopsis genome-wide cold response.
Supplemental Data Set 4. Gene lists of known Arabidopsis cold response genes in Supplemental Data Set 5.
Supplemental Data Set 5. Novel genes differentially regulated by cold at the expression (DE) and/or alternative splicing level (DAS).
Supplemental Data Set 6. Gene Ontology enrichment analysis of DE genes and DAS genes.
Supplemental Data Set 7. Genes in DE and DTU heat map clusters and Gene Ontology enrichment analysis of individual heat map DE gene clusters.
Supplemental Data Set 8. Gene lists and descriptions of DE and DAS genes organized by the time point at which they first become significantly differentially expressed upon cold.
Supplemental Data Set 9. DAS gene lists and descriptions organized by different ΔPS cutoff values.
Supplemental Data Set 10. Isoform switch and AS events analyses of transcripts from DAS genes.
Supplemental Data Set 11. List of Arabidopsis Transcription Factor and SF-RBP genes.
Supplemental Data Set 12. Summary of TF and SF-RBP DE and DAS genes with the largest and most rapid changes in AS.
Supplemental Data Set 13. Timing of significantly DE and DAS SF-RBP and TF genes in the first 12 h after onset of cold treatment (T9-T13).
Supplemental Data Set 14. AS data for the u2b’’-like HR RT-PCR analysis panel.
Supplemental Data Set 15. Summary of RNA-seq data derived from Arabidopsis rosettes time-series cold response study.
Supplemental Data Set 16. Transcript expression in the Arabidopsis cold response.
Supplemental Data Set 17. Primers used in this study.
Dive Curated Terms
The following phenotypic, genotypic, and functional terms are of significance to the work described in this paper:
Acknowledgments
This work was supported by funding from the Biotechnology and Biological Sciences Research Council (BB/K006568/1, BB/P009751/1, and BB/N022807/1 to J.W.S.B.; BB/K006835/1 to H.G.N.; and BB/M010996/1-EASTBIO Doctoral Training Partnership to J.C.E.) and by the Scottish Government Rural and Environment Science and Analytical Services Division (to J.W.S.B. and R.Z.). We acknowledge the European Alternative Splicing Network of Excellence (EURASNET; LSHG-CT-2005-518238) for catalyzing important collaborations. We thank Janet Laird (University of Glasgow) for technical assistance, and Katherine Denby and Iulia Gherman (University of York) for the list of Arabidopsis transcription factors. RNA-seq was performed at The Genome Analysis Centre (now Earlham Institute), Norwich, UK. We thank Robbie Waugh and Piers Hemsley for critical reading of the manuscript and helpful comments. We apologize to all the authors whose relevant work was not cited in this article due to space limitations.
AUTHOR CONTRIBUTIONS
J.W.S.B., H.G.N., and R.Z. designed the research and are co-corresponding authors. C.P.G.C., W.G., A.B.J., N.A.T., J.C.E., P.E.P., H.K., and R.Z. performed the research. C.P.G.C., W.G., R.Z., and J.W.S.B. analyzed the data. C.P.G.C. and J.W.S.B. wrote the manuscript with contributions from all others.
Footnotes
↵1 These authors contributed equally to this work.
- Received February 26, 2018.
- Revised April 20, 2018.
- Accepted May 10, 2018.
- Published May 15, 2018.
This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
References
