|Home | About | Journals | Submit | Contact Us | Français|
mRNA levels in cells are determined by the relative rates of RNA production and degradation. Yet, to date, most analyses of gene expression profiles were focused on mechanisms which regulate transcription, while the role of mRNA stability in modulating transcriptional networks was to a large extent overlooked. In particular, kinetic waves in transcriptional responses are usually interpreted as resulting from sequential activation of transcription factors.
In this study, we examined on a global scale the role of mRNA stability in shaping the kinetics of gene response. Analyzing numerous expression datasets we revealed a striking global anti-correlation between rapidity of induction and mRNA stability, fitting the prediction of a kinetic mathematical model. In contrast, the relationship between kinetics and stability was less significant when gene suppression was analyzed. Frequently, mRNAs that are stable under standard conditions were very rapidly down-regulated following stimulation. Such effect cannot be explained even by a complete shut-off of transcription, and therefore indicates intense modulation of RNA stability.
Taken together, our results demonstrate the key role of mRNA stability in determining induction kinetics in mammalian transcriptional networks.
mRNA levels in cells are determined by the relative rates of RNA production and degradation. Transcript levels at steady state therefore reflect equilibrium of RNA synthesis and decay. Gene expression microarrays, and more recently RNA deep-sequencing, are valuable means for genome-wide profiling of the cellular transcriptome and its modulation during normal development and in response to external perturbations. Standard microarray analyses on total cellular RNA provide a measure of mRNA abundance but cannot discriminate whether changes are due to alterations in RNA transcription or decay. New techniques have been developed to allow the measurement of actual rates of transcript production. When these are carried out in parallel to recording overall RNA abundance, they can reveal the relative contribution of alterations in gene transcription and mRNA stability to the observed net change in RNA abundance [1,2].
Despite the above, to date, the contribution of RNA degradation to global changes in cellular transcriptome was largely overlooked. This was mostly due to the fact that most microarray studies, explicitly or implicitly, ascribed alterations in RNA levels to correlated changes in gene transcription. However, several pioneering studies have shed light on the critical role that modulation of mRNA stability plays in the regulation of cellular transcriptome in response to various stresses [3-7].
Time-course analysis is a very common design for microarray analysis, which allows researchers to follow the dynamics of the cellular response to perturbations. Clustering analysis applied to time-course data partitions the response into distinct kinetic waves, distinguishing between early-, intermediate- and late- responding genes. Such kinetic waves can be the result of sequential activation of primary and secondary transcription factors (TFs) [8,9]. Nevertheless, mathematical modeling of changes in RNA levels predicts that mRNA degradation rate plays a pivotal role in shaping the kinetics of genes' response [10,11]. A standard measure for the speed of a transition between two steady states is T1/2, which is the time at which half of the change between the new and the former steady-state levels is achieved. (In the phase of mRNA decay T1/2 is usually referred to as 'T half-life' which measures the period of time it takes for a transcript undergoing decay to decrease its level by half). Importantly, a simple kinetic model predicts that T1/2 is determined by RNA degradation rate not only in decay of expression but also in the phase of induction. The model assumes that mRNA is produced at a constant rate (β) while the rate of its degradation is proportional to its concentration. Accordingly, the rate of change in mRNA concentration (X) is given by the equation: dX/dt = β-αX (α denotes the degradation rate constant). At steady state mRNA concentration reaches equilibrium (that is, there is no change: dX/dt = 0), and therefore, the steady state level (Xss) is determined by the ratio between the synthesis and degradation constants: Xss = β/α. Solving the above equation gives the change in mRNA concentration over time: ΔX(t) = [β/α-X0]*(1-e-αt) (X0 represents the mRNA concentration at t0 where the perturbation was applied to the system, assuming that change in transcription rate occurs instantaneously). Of note, according to this solution, the rapidity of a transition between former and new steady states (T1/2) is determined only by α, and is inversely proportional to it: T1/2 = ln2/α; Therefore, genes with unstable RNA are predicted to respond in fast kinetics, whereas genes with stable RNA should respond more slowly (Figure (Figure1;1; For a thorough discussion of the kinetic model see ). Furthermore, this model predicts that the time it takes a transcript whose transcription rate is increased by a factor L to achieve a k-fold induction in expression is a linear function of T1/2: Tk = -log2(1-f)*T1/2; where f = (k-1)/(L-1) ; That is, for a similar increase in transcription rate in response to a stimulus, genes encoding unstable mRNAs are predicted to achieve a certain fold of induction faster than genes which encode stable transcripts.
Recently, the critical role for mRNA stability in the induction kinetics of genes encoding inflammatory proteins was demonstrated at single-genes level . Here, we set out to examine this role on a large-scale utilizing a global atlas of mRNA stability recently generated in mammalian cells and analyzing numerous gene expression time-course datasets that collectively cover many different aspects of the cellular physiology.
To examine the role of mRNA stability in determining the dynamics of transcriptional networks in mammalian cells, we scanned expression microarray repositories for time-course datasets that profiled cellular transcriptome in human and murine cells at dense kinetics. As a test case, we first analyzed a dataset that recorded expression profiles at multiple time points (0, 0.5, 1, 2, 4, 6, 8, and 10 hrs) after stimulation of murine T cells with IL2 . We partitioned the genes that were induced by this treatment into distinct kinetic clusters according to their response time (i.e., according to the first time point in which a gene was induced above a certain fold-change threshold) (Figure (Figure2A).2A). A clear distinction between early-, intermediate- and late- responding genes emerged.
Next, we tested if there is a relationship between the observed kinetics of gene induction and mRNA stability. Namely, we tested whether fast responding genes were characterized by high instability (low T half-life), while late responders were more stable (high T half-life). In these tests we utilized a global atlas of mRNA stability in mammalian cells that was recently generated in murine fibroblast and human B cells [2,12]. Importantly, these reports noted that mRNA half-life times were generally well conserved between the examined cell types and species (the median T1/2 in human and mouse cells was 315 min and 274 min, respectively). Using this data source, we found a striking correlation between response time and transcript stability: T half-lives of early-induced genes were significantly shorter than those of late-induced genes (Figure 2B, C), as predicted by the kinetic model.
Next, we examined the generality of the association between kinetics of induction and RNA stability. We analyzed a variety of time-course expression datasets recorded in human and murine cells, collectively covering many different aspects of cellular physiology. In the vast majority of the datasets analyzed, we found a highly significant anti-correlation between rapidity of induction and mRNA stability (Figure (Figure3,3, Additional file 1, Additional file 2). This widespread relationship points to the critical role played by mRNA stability in shaping the dynamics of gene induction in complex transcriptional networks. It also indicates a broad conservation of RNA stability under different conditions, which therefore reflects, to a large extent, an intrinsic property of the mRNA molecules. Inspection of the early-induced genes revealed a core set of genes whose induction-response was very rapid in many different datasets and which encode highly unstable transcripts (e.g., Fos, Jun, Ier3, Dusp1, Atf3, Btg2 and Zfp36; all have T half-life lower than 1 hr). (Additional file 3 lists the core set genes, defined as the set of genes that were induced before or at 2 hrs after stimulation in at least three of the six datasets recorded in murine cells.) However, the broad relationship between mRNA stability and induction kinetics is not merely explained by this common core set of rapidly induced genes, as the relationship remained highly significant also after the removal of this core set from the analysis (Additional file 4, data not shown). This indicates that many other genes with an unstable mRNA were rapidly induced in a stimulus-specific manner.
mRNA stability is regulated mainly by cis-regulatory elements embedded in the transcript 3'-UTR [5,15]. Seeking for major mechanisms that control stability in the datasets that we analyzed, we searched for enriched sequence patterns in the 3'-UTRs of the core set of early-induced genes. In agreement with previous reports [15,16], we found that these highly unstable mRNAs are significantly enriched for the AU-rich element (ARE); the most highly enriched 7-mer in the 3'-UTRs of these genes was UAUUUAU, which appeared in 52% of the 3'-UTRs in this set compared to background frequency of 18% in all 3'-UTRs of mouse genes (p-value = 7.6*10-9, hypergeometric tail).
Another physical factor which limits the rapidity at which genes are induced is the genomic transcribed length. The earliest time in which a transcript can be induced is bound by the length of its encoding gene and the velocity at which the RNA polymerase elongates along it. Therefore, we expected that genes which are induced at a very rapid kinetics would be characterized by short genomic transcribed length. Indeed, in the vast majority of the datasets that we analyzed we detected a significant correlation between rapidity of induction and genomic transcribed length (Figure (Figure4,4, Additional file 1 and Additional file 5). However, in contrast to mRNA stability that affected the induction kinetics over a wide range of time points, the effect of genomic transcribed length was noticeable only in the very early time points after stimulation (compare Figure Figure2B2B and Figure Figure4).4). Furthermore, no correlation was observed between response time and length of mature transcripts (that is, the length of mRNA transcripts after introns are spliced out; data not shown), and there was no overall correlation between mRNA stability and genomic transcribed length (data not shown). Yet, the core set of early induced genes was characterized by both very short T half-lives and very short genomic transcribed length (Additional file 6). Therefore, to achieve very rapid gene induction, a combined strategy that couples substantial increase in transcriptional rate with rapid degradation rate and short genomic transcribed length is undertaken. Of note, this core set is enriched for TFs (16 out of 52 genes, p = 4.7*10-5 after FDR correction), as agile regulation of the regulators is a critical property of networks.
The mathematical kinetic model also predicts a similar relationship between RNA stability and kinetics of gene suppression. Genes whose expression is down-regulated rapidly are expected to encode RNAs with lower stability than genes whose expression is down-regulated at slower rate (Figure (Figure1).1). Interestingly, while in several datasets we observed a very good agreement with this expectation, deviations from it were frequent. In five out of the ten datasets that we analyzed, we observed no relationship between mRNA stability and kinetics of gene suppression, while this relationship was highly significant for the induced genes (Additional file 1, Additional file 7). In those cases, the expression of many mRNAs which have high T half-life under standard conditions (above 5 hrs) was already down-regulated by at least a factor of 2.0 after only 1-2 hrs post treatment. Such a rapid down-regulation of stable RNAs cannot be explained even by a complete shut-off of transcription, and therefore suggests intense modulation of mRNA stability of the rapidly suppressed genes. We could not detect statistically enriched sequence motifs in the 3'-UTR of these genes that might point to the mechanism (e.g., microRNA, RNA-binding protein) which controls this stability modulation.
Two mechanisms which underlie the temporal order of gene induction are sequential activation of primary and secondary TFs and mRNA stability. These mechanisms operate in parallel, and therefore the observed induction kinetics is the result of their superposition. The kinetic model we used in our study is oversimplified as it assumes constant rates of transcription and degradation over time. In practice, however, most genes are regulated by modules of transcriptional activators and repressors rather than by single TF, the activity of TFs themselves is modulated over time and they often form interlocked feedback loops. Therefore, the kinetic patterns exhibited by responding genes are much more complicated than the simple exponential pattern predicted by the model. Yet, much insight can be gained on the dynamics of transcriptional networks using the simplified description. While many studies delineated the way various transcriptional networks are propagated by sequential cascades of TFs ([9,17,18], much less attention was given to the role of mRNA stability in shaping the induction dynamics. Recently, Hao and Baltimore  showed that mRNA stability significantly influences the induction kinetics of genes encoding inflammatory proteins. Here, by analyzing numerous gene expression datasets that collectively cover many different aspects of cellular physiology, we demonstrated on a global scale, the critical role that mRNA stability plays in shaping the dynamics of transcriptional networks. This global relationship agrees with the simple kinetic model for mRNA concentration which predicts that unstable transcripts respond faster than stable ones.
Of note, the striking anti-correlation that we detected in a diverse panel of cell types between mRNA stability and induction kinetics was derived using an atlas of mRNA half-lives measured in murine fibroblasts and human B-cells. This result indicates a broad conservation of RNA stability under different cell types and conditions, which therefore reflects, to a large extent, an intrinsic property of the mRNA molecules. This observation supports and strengthens Friedel et al. conclusion that mRNA stability is largely conserved between cell types and species .
Another factor which limits gene induction time is the genomic transcribed length. In accord, we observed that genes that were induced very rapidly were significantly short. Therefore, to achieve very rapid gene induction, substantial increase in transcriptional rate is coupled with rapid degradation rate and short genomic transcribed length, as exhibited by the core set of early-induced genes.
Interestingly, while we observed very good agreement with the kinetic model when analyzed gene induction, we found major deviations from the model predictions when analyzed gene suppression. In response to many stimuli, we observed very fast down-regulation of mRNAs whose half-life time is very high under normal conditions (as measured in fibroblasts (mouse) and B-cells (human)). Even a complete turn-off of transcription is not enough to achieve such a rapid reduction in the concentration of stable mRNAs. Therefore, those mRNAs that need to be down-regulated at a very high speed require regulatory mechanisms that decrease their stability. Key regulators of mRNA stability are RNA binding proteins (RBPs) and microRNAs (miRs), and ample information on the activation of these regulators in response to various stresses has been already accumulated. We speculate that activation of stimulus-specific RBPs and miRs is a major factor that underlies the deviation from model's prediction in the case of gene suppression.
Comprehensive understanding of gene expression networks can be gained only once we obtain a global delineation of the orchestrated modulation of transcription and degradation rates carried out by cells in normal development and in response to perturbations. In this study, we elucidated the key role of mRNA stability in shaping the kinetics of gene induction in intricate gene networks in mammalian cells.
Expression data were downloaded from public repositories (GEO and ArrayExpress). All datasets used Affymetrix arrays. Expression levels were calculated using the rma method  (implemented in Affymetirx Expression Console tool). For each dataset, presence flags were calculated using MAS5, and only probe-sets that were flagged as 'Present' in at least two chips were retained for subsequent analysis. For genes represented by multiple probe-sets, we chose the one probe-set with the highest median intensity in the dataset.
Kinetic clustering and all statistical analyses of T half-life and genomic transcribed length distributions were done in R.
Transcript sequences and lengths for all human and mouse genes were obtained using Biomart (Ensembl v54) . For genes encoding multiple transcripts, the transcript with the longest genomic transcribed length was used in length distribution tests.
Enrichment test for k-mers in 3'-UTR of murine genes was carried out using the AMADEUS tool . For genes with multiple forms of 3'-UTR, the longest one was used. Enrichment test for GO functional categories was carried out using DAVID web-service . In both, the core set of early-induced genes was compared to a background set of all murine genes.
RE and RA conceived and designed the study. RE carried out the statistical analyses. EZ and KZ performed gene expression experiments. RE and RA wrote the manuscript. All authors read and approved the final manuscript.
Relationship between mRNA stability, response kinetics and genomic transcribed length in ten expression datasets.
Core set of early-induced genes.
Examination of the relationship between mRNA stability and kinetics of induction in the IL2 dataset after removing from the analysis the core set of early induced genes. p-value was calculated for the comparison between the distribution of T half-life of early and late induced genes, as done in Additional file 1, but after the removal of the core genes.
Relationship between kinetics of induction and genomic transcribed length. The effect of genomic transcribed length on the response time is evident only at the very early time points (up to 1-2 hrs after stimulation; see legend of Additional file 1).
The core set of early induced genes is characterized by both (a) very short T half-life (mean T1/2 of 0.90 h vs. 6.45 h, in the core and background sets, respectively) and (b) very short genomic transcribed length (mean genomic transcribedlength of 9,474 bp vs. 39,789 bp, in the core and background sets, respectively). P-values (Wilcoxon test) were calculated for the comparison between the core set and a background set which contained all the rest of genes for which T half-life and genomic transcribed length (i.e., CDS and UTRs annotations) data are available.
We thank Chaim Linhart for critical reading of the manuscript. RE is supported by an EMBO long-term fellowship.