INTRODUCTION
============

Although the abundance of messenger RNAs (mRNAs) in the cell is the net result of mRNA synthesis and degradation, expression profiling experiments usually take a snapshot of the mRNA levels in the cell and do not capture the dynamics of mRNA synthesis and breakdown. This lack of information on the dynamic component of the regulation of mRNA levels limits biological investigations in general and proper modeling of transcriptional networks in systems biology approaches in particular. When mRNA synthesis is shut down completely, the effect may be immediate for very short-lived mRNAs but only manifest after hours or days for mRNAs with long half lives. The current study fills in part of this gap and surveys mRNA degradation rates during muscle differentiation on a genome-wide scale. With the availability of these data, the genetic network orchestrating muscle differentiation can be revised and improved.

From the results of previous genome-wide analyses of mRNA degradation rates, we know that there are large differences between the stabilities of individual mRNAs, and that mRNA stability is strongly related to the molecular function of the encoded protein ([@B1; @B2; @B3; @B4; @B5; @B6; @B7; @B8]). However, the mechanisms underlying differential stability of transcripts remain poorly understood. MicroRNA (miRNA)-mediated mRNA degradation emerged as an important controlling factor for mRNA half life. Earlier research also characterized a plethora of RNA binding proteins that may affect mRNA stability ([@B9]). Recent studies in differentiating embryonic stem cells ([@B7]) and T-cells ([@B6]) demonstrated that mRNA stability may be affected by the cell's biological state. Furthermore, an important subgroup of transcripts use alternative polyadenylation sites depending on the differentiation status ([@B10]). This results in 3′-untranslated regions (3′-UTRs) with variable lengths and probably contributes to the differential stability of these transcripts ([@B10]). The deciding factors determining the choice for shorter or longer 3′-UTRs are currently unknown.

We were interested in the control of decay rates during myogenic differentiation and studied proliferating and differentiated myogenic C2C12 cells, a well-established system to study differentiation of muscle cells. Cells were treated with actinomycin D to arrest *de novo* RNA synthesis, with sampling at short time intervals. We used Affymetrix whole-genome exon arrays covering all known and predicted exons. Where previous genome-wide studies determined expression levels and stability of the 3′-ends of transcripts, these arrays allowed us to assess expression level per exon and thereby, for the first time, the effect of alternative splicing on mRNA decay rates. Our results demonstrated that changes in mRNA abundance during differentiation of muscle cells are correlated with changes in decay rates and that the ratios of specific splice variants are controlled at the level of mRNA stability.

MATERIALS AND METHODS
=====================

Cell culture, RNA isolation, microarray hybridization
-----------------------------------------------------

C2C12 cells were grown in proliferation medium \[Dulbecco's modiﬁed Eagle's medium (DMEM) with 10% fetal bovine serum (FBS), 2% Glutamax, 1% glucose and penicillin/streptomycin\]. Myogenic differentiation was induced by shifting to differentiation medium (DMEM with 2% FBS, 2% Glutamax, 1% glucose and penicillin/streptomycin). C2C12 cells were maintained in differentiation medium for 8 days, when long and multinucleated myotubes had formed. Cultures of proliferating and differentiated C2C12 muscle cells in 6-well dishes containing 3-ml medium were treated with 5 µg/ml actinomycin D to arrest the transcription machinery. At seven different time points after addition of actinomycin D (0, 10, 20, 30, 60, 150, 480 min), cells were lysed and RNA was isolated with the Machery Nagel Nucleospin RNAII Kit. There was no notable decrease in cell viability within this time period. The experiment was performed in duplicate. After depletion of ribosomal RNAs, mRNA was reverse transcribed with random primers containing a T7 tag, and labeled with biotin during the T7 RNA amplification (which retains strand specificity), according to standard Affymetrix protocols. Equal amounts of amplified RNA were hybridized to 28 Affymetrix Mouse Exon v1 arrays \[7 time points × 2 replicates × 2 conditions (proliferating/differentiated)\].

Affymetrix exon arrays
----------------------

Affymetrix Mouse Exon v1 arrays contain 1 236 087 non-control probe sets. Most probe sets consist of four probes and the intensity value of each probe set is the robust mean of the individual probe sets. Most exons are covered by at least one probe set. There is a subset of exons, called core exons, which represent all the exons for which there is substantial experimental proof for their existence. Next to the core set, the array covers a huge set of computationally predicted exons. The intensities of probe sets in most of these latter exons (in our experiment 91 versus 51% for core exons) are not above background, and these have not been considered in the analysis described here. Data from these arrays can be analyzed at the exon or the gene level. The presented gene level summarization, giving only one expression value per gene, considers the intensities of probes in core exons only.

Data analysis---normalization
-----------------------------

In our experimental set-up, normalization is not straightforward, since an important assumption underlying standard normalization procedures, equality of the total amount of mRNA under all conditions studied, is violated. Despite the decrease in the total amount of mRNA in the cell at later time points, equal amounts of RNA have been hybridized to the chips. This means that the samples from the later time points will be enriched for stable transcripts and depleted for fast decaying transcripts. We tried to use spike-in RNAs to correct for the differences in input RNA but this was not successful due to the low number of spike-in RNAs used. Thus the only alternative was to use a standard normalization method, with the intensities of more stable transcripts increasing at later time points.

The median polish summarization method \[as implemented in RMA ([@B11])\] was applied to summarize the intensities of probes in a probe set and probes in a gene. VSN \[variance stabilization and normalization ([@B12])\] was used for normalization. VSN and RMA normalization generated similar values in the high intensity range. In the lower intensity range, the variance stabilizing properties of the VSN algorithm corrects for noise-induced artificially inflated differences between arrays in near-background reads. As an additional advantage, the VSN data are on a natural logarithm scale compatible with our exponential decay model, whereas the RMA data are on a log~2~ scale. Since the arrays contain many negative and positive control spots, the effect of removing those before the normalization was investigated. This had no effect on the normalization of the non-control spots. As a threshold for background intensity levels, we used the mean + 2 SD of around 4000 negative control probes on the array.

All normalization and subsequent data analysis was performed in the statistical programming language R ([@B13]).

Raw and normalized data have been submitted to GEO and are available under series GSE14387.

Data analysis---model fit
-------------------------

To model decay rates, a standard linear model ([@B1]) was fitted with normalized intensity *Y* of each probe set *j* as dependent variable and time as independent, continuous variable. where time is a vector (0, [@B10], [@B20], [@B30], 60, 150, 480); α*~j~* represents the intercept (estimated expression level at *t* = 0); β*~j~* represents the slope (negative of the decay constant).

The appropriateness of this model is discussed in the [Supplementary Methods](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1) and in the caption of [Supplementary Figure S1](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1). Replicate time points were independent data points in the model. Two separate models were fit for the proliferating and differentiated cells. Decay constants were said to be significantly different between proliferating and differentiated cells when there was no overlap between their 95% confidence intervals (as determined from the standard error on *β*, and given the *t*-distribution and 12 degrees of freedom). Decay constants for entire transcripts were calculated using the exact same model, but after summarizing core probe sets in a transcript into a single expression value by RMA. The half life is given by --ln([@B2])/*β*.

For identification of probe sets with a behavior different from the other probe sets in the gene (possible alternative splicing events), we considered the mixed effect model ([@B2]) (nlme package; lme function in *R*), where time is a continuous variable and treated as a fixed effect, where the probe sets were considered as random effects and where we included an interaction term between the probe sets and the differentiation status. In formula form ([@B14]): where *b~j~* ∼ *N* (0, σ~1~^2^), *b~jk~* ∼ *N* (0, σ~2~^2^), ε*~ijk~* ∼ *N* (0, σ^2^) and where *Y~ijk~* represents normalized intensity of probe set *j* (*j* = 3, ... , *n* probe sets in the gene) with differentiation status *k* (*k* = 1, 2) in gene *i*; time represents time vector (0, [@B10], [@B20], [@B30], 60, 150, 480); *β* represents fixed effect; *b~j~* represents random intercept; *b~jk~* represents the interaction between the probe set and differentiation status and ε*~ijk~* represents error term.

Data analysis---annotation and visualization
--------------------------------------------

The R package Xmap ([@B15]), interfacing to a local Ensembl install (mouse release 47; NCBI genome build 37) was used for annotation and visualization of the exon-level data.

Quantitative PCR analysis
-------------------------

Quantitative polymerase chain reaction (qPCR) analysis with SYBR-Green detection was performed as described earlier ([@B16]). *Bgn* was used in the RNA decay studies as a reference gene to correct for differences in input cDNA. *Bgn* was found to be extremely stable and did not change significantly in Ct value over the time period analyzed. *Gapdh* was used as a reference gene for the nuclear versus cytoplasmic comparison.

3′-RACE
-------

To determine the polyadenylation sites in the *Itga7* transcript, we performed 3′-RACE (rapid amplification of cDNA ends). We performed a reverse transcription (RevertAid RNaseH-M-MuLV, MBI-Fermentas) on 500 ng of total RNA starting at the polyA site with an anchored primer oligo-dT primer (sequence: GCTCGCGAGCGCGTTTAAACGCGCACGCGTTTTTTTTTTTTTTTTTTVN). Subsequently, we performed two rounds of PCR \[30 cycles, Faststart Taq (Roche)\]. The forward primer was designed in exon 22 (sequence: CTGTTAGTCCTGGCCTTGCT); the sequence of the reverse primer in the first round was GCTCGCGAGCGCGTTTAAAC, whereas the sequence of the reverse primer in the second round of nested PCR was GCGTTTAAACGCGCACGCGT. Bands were excised from agarose gel and sequenced on an ABI-3700 sequencer.

Isolation of nuclei
-------------------

Cells were scraped in phosphate buffered saline and pelleted by centrifugation for 5 min at 1500*g*. The cells were resuspended in cold lysis buffer (10 mM Tris--HCl pH 7.4, 3 mM MgCl~2~, 10 mM NaCl, 1mM EDTA) and pelleted by centrifugation for 5 min at 1500*g*. After resuspending in lysis buffer including 0.5% (v/v) Triton X-100, cells were disrupted in a Dounce homogenizer with 4 × 10 strokes of the pestle. Nuclei were pelleted by centrifugation for 2 × 5 min at 1500*g* and washed with lysis buffer (without Triton X-100). Nuclei were purified by sucrose gradient centrifugation. To this end, a centrifuge tube was filled with 4 ml of high sucrose buffer (2 M sucrose, 5 mM MgAc, 0.1 mM EDTA, 10 mM Tris--HCl pH 8.0, 1 mM DTT) and nuclei were added after mixing with 4.4 ml of low sucrose buffer \[0.32 M sucrose, 3 mM CaCl~2~, 2 mM MgAc, 0.1 mM EDTA, 10 mM Tris--HCl pH 8.0, 1 mM DTT, 0.5% (v/v) Triton X-100\], with subsequent centrifugation for 30 min at 30 000*g* in an SW41 rotor. Pellet was resuspended in 300 µl of storage buffer \[50 mM Tris--HCl pH 8.3, 40% (v/v) glycerol, 2 mM MgCl~2~, 0.1 mM EDTA\], aliquotted and snap frozen in liquid nitrogen.

RESULTS
=======

Experimental setting
--------------------

Cultures of proliferating C2C12 muscle precursor cells and cultures of C2C12 cells differentiated to elongated and fused myotubes were treated with actinomycin D to arrest the transcription machinery. At seven different time points after addition of actinomycin D (0, 10, 20, 30, 60, 150, 480 min), RNA was isolated and hybridized to Affymetrix Mouse Exon v1 arrays (in duplicate, 28 arrays in total). With the selected time intervals, it is possible to accurately estimate the degradation rates of relatively short-lived transcripts. For very stable transcripts (*t*~1/2~ \> 480 min) the estimate becomes inaccurate, but longer incubation times resulted in a significant decrease in cell viability severely complicating the analysis.

To model mRNA degradation, we consider an exponential decay model. This model becomes linear on the logarithmic scale used to report the normalized exon-array signals (model 1; see [Supplementary Methods and Supplementary Figure S1](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1) for details on the model and the goodness-of-fit). The slope in this linear model is the (negative of the) decay constant, whereas the intercept reflects the expression level before addition of actinomycin D.

Estimation of decay rates
-------------------------

All mRNAs decay over time. The normalization method used (see 'Materials and Methods' section), centers the average slopes around zero ([Figure 1](#F1){ref-type="fig"}A). Thus, stable transcripts appear to have positive slopes relative to the average, whereas transcripts with fast degradation will have negative slopes. The relative differences in decay rates (slow decaying versus fast decaying) can be well estimated using this model ([Figure 1](#F1){ref-type="fig"}B, [Supplementary Figure S1 and Supplementary Methods](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1)). The maximum slope for expressed genes is 0.0028 intensity units/minute, the mean is −0.00042 (median −0.00020), and the minimum −0.0073. If we assume that the slowest decaying transcript is stable over the period analyzed, all slopes would need to be corrected with −0.0028. The half life is calculated as −ln([@B2])/slope and the shortest half life we observed is then 69 min, the median half life of all expressed transcripts 231 min. This is in a similar range as reported before for proliferating C2C12 cells ([@B8]). The decay rates determined on the microarrays correlate well with those measured by quantitative PCR \[[Figure 2](#F2){ref-type="fig"}, Pearson's correlation coefficient: 0.94 (*P* = 3.2 × 10^−6^)\], in spite of the array data appearing to overestimate the half life of short-lived transcripts and to underestimate those of stable transcripts. Relative decay rates for all genes with detectable expression (6852 in total) are reported in [Supplementary Table S1](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1). Figure 1.Estimated relative decay rates measured on microarrays. Distribution of the slopes (negative of degradation rate constants) (**A**) and *R*^2^ values of the linear model (**B**) for all probe sets (dashed line) and for probe sets in genes expressed above background (solid line). The majority of probe sets do not give detectable signals and have very flat (slope near zero) and noisy (low *R*^2^ values) decay curves. In contrast, expressed genes frequently demonstrate more positive or negative slopes, and their decay rate constants can be reliably measured (high *R*^2^ values). This is expected since background signals remain constant over time, thus generating a slope close to zero, whereas true signals tend to increase or decrease, generating a slope away from zero. Figure 2.Correlation of microarray- and quantitative PCR-based degradation rates. For both microarray and qPCR assays, the natural logarithm of the expression level was plotted against the time after start of actinomycin D treatment and a linear curve was fitted. The slope of this curve is the negative of the decay constant. The slopes for the *Bmp4, Ccng2, Cebpg, Pde4dip-short, Pde4dip-long* and *Itga7-exon27* genes in proliferating and differentiated myoblasts measured with qPCR with primers designed in the 3′-UTR are plotted on the *x*-axis, and the slopes measured with microarray probe sets in the same regions are plotted on the *y*-axis. The correlation coefficient is the Pearson correlation coefficient.

Factors determining mRNA half life
----------------------------------

It has been noted that mRNA stability appears to correlate with the cellular function of the encoded proteins ([@B1; @B2; @B3; @B4; @B5; @B6; @B7; @B8]). To study this correlation in our data we functionally classified the most stable and unstable mRNA transcripts in proliferating and differentiated cells ([Table 1](#T1){ref-type="table"}). We see a similar correlation. Stable RNAs code for proteins involved in cell maintenance and structure, short-lived RNAs code for proteins involved in regulation of gene expression, including histones. RNAs coding for structural and contractile muscle proteins, which are only expressed in differentiated cells, appear to be very stable. Transcripts in the category 'muscle development' are either very short- or very long-lived, because this category contains both regulatory and structural proteins. Table 1.Functional classification of fast and slow decaying transcripts in proliferating and differentiated cellsFast*P*-valueSlow*P*-valueProliferatingRegulation of cellular processes3.7 × 10^−8^DNA-dependent DNA replication3.6 × 10^−4^Establishment and/or maintenance of chromatin architecture7.1 × 10^−4^DifferentiatedRegulation of cellular procesess4.6 × 10^−10^Muscle development1.1 × 10^−4^Regulation of transcription8.3 × 10^−5^Muscle contraction7.6 × 10^−28^Negative regulation of protein kinase activity3.1 × 10^−7^Cell cycle5.0 × 10^−5^Development2.1 × 10^−7^Muscle development1.3 × 10^−4^[^1]

Short mRNAs tend to be less stable than average ([Supplementary Figure S2](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1)). The median half lives of mRNAs consisting of one, two, three exons in proliferating cells are 172, 189 and 199 min, respectively. This is significantly shorter (*P* ≤ 1 × 10^−5^ Kolmogorov--Smirnov test) than the average half life in proliferating cells (237 min). Similar results were obtained for differentiated cells. [Figure 3](#F3){ref-type="fig"} visualizes the statistically significant shift to shorter half lives in mRNAs with three, two or one exons. Figure 3.One, two and three exon mRNAs have shorter half lives. Distribution of slopes (negative of degradation rate constants) for all mRNAs (dashed line), and mRNAs consisting of one (solid line), two (dash-dotted line) or three (dotted line) exons. The latter three distributions demonstrate statistically significant shifts to the left (shorter half lives), with Kolmogorov--Smirnov *P*-values of 4.3 × 10^−9^, 1.3 × 10^−7^ and 1.2 × 10^−5^, for one, two and three exons, respectively. mRNAs with 4 exons were only marginally different from the overall distribution (*P* = 0.015). The plot is for proliferating cells only, but the plot and results for differentiated cells are highly similar.

There is no correlation between 3′-UTR length and transcript stability ([Supplementary Figure S3](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1)).

mRNA stability is dependent on the cellular state
-------------------------------------------------

We calculated the differences in decay rates between differentiated and proliferating cells. There were 170 mRNAs displaying significant differences in stability between the two states; 99 of those were more stable in proliferating cells and 71 were more stable in differentiated cells ([Supplementary Table S1](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1)). Interestingly, many of the mRNAs with differences in stability were also differentially expressed between the two conditions. [Table 2](#T2){ref-type="table"} summarizes some of the results for transcripts playing an established role in muscle development. Slower degradation in differentiated than in proliferating cells goes with higher expression levels in differentiated cells, and vice versa. *Myog* and *Cdkn1a* (p21), two important regulators of differentiation previously shown to have increased expression and stability in differentiated cells ([@B17; @B18; @B19]) are amongst the genes with clearly different mRNA degradation rates: The half life of the *Myog* mRNA increased from 165 min in proliferating cells to 268 min in differentiated cells; the half life of the *Cdkn1a* mRNA increased from 138 to 371 min. The generality of these observations is illustrated by the strong positive correlation obtained after plotting the differences in gene expression against the differences in decay constants ([Figure 4](#F4){ref-type="fig"}, Pearson's correlation coefficient: 0.57; *P* \< 2.2 × 10^−16^). Thus, control of mRNA stability and degradation emerges as an important determinant for the regulation of steady-state mRNA levels. Figure 4.Correlation between differential expression and differences in decay rates. Correlation of the differences in slopes (per minute) between differentiated and proliferating cells (*y*-axis) and differences in normalized expression levels between differentiating and proliferating cells (*x*-axis) for gene level data (Pearson's correlation: *R* = 0.57). Transcripts that are higher expressed under differentiating conditions (right side of the plot) seem to have less negative slopes (upper part of the plot) and thus to be more stable under differentiating conditions. Probe set level data give similar results but are more difficult to visualize since the plot contains many more data points. Table 2.Half lives of myogenesis-related genes change during differentiationGeneFold-change differentiated versus proliferating cellsHalf life proliferating cells (min)Half life differentiated cells (min)*Tnc*−3.86444229*Mbnl3*−2.72342264*Musk*2.82246438*Lrp4*3.29236416*Cdkn1a*3.73138371*Tnnt2*4.75242500*Tmod1*7.23201455*Acta1*7.75322555*Myog*9.09165268*Cacna1s*17.01229798[^2]

Relative abundance of splice isoforms is controlled by mRNA stability
---------------------------------------------------------------------

The observation that the cellular state determines the overall decay rates of transcripts triggered us to investigate the role of mRNA degradation in controlling the ratio between different splice isoforms. To this end, we used the mixed effect model ([@B2]) to identify genes where the decay constant of a particular probe set differed more than 2 SDs from the decay constant of all other probe sets in the gene. Sixty-nine probe sets met these criteria ([Supplementary Table S2](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1)). After further filtering for probes that show a low intensity across all conditions, we were left with 39 probe sets in 36 genes meeting our very stringent criteria. The positions of these probe sets in the mRNAs were evaluated ([Table 3](#T3){ref-type="table"}). As expected, we found 10 probe sets in 3′-UTRs, the length of which is known to affect transcript stability ([@B10]). The mRNAs include *Ankrd2, Myom2, Itgb1, Cacna1s, Pde4dip* and *Xirp1*, with a known function in muscle development, highlighting the importance of the shift in 3′-UTR usage during myogenic differentiation. Another 13 probe sets were located in the 5′-UTR and one in the coding part of a first exon. Given the known variability in the use of transcription start sites ([@B20]), this implicates that the choice for a transcription start site and thus the length of the 5′-UTR may influence the stability of the transcript. A further 15 probe sets were located in internal exons, many of them known to be subjected to alternative splicing. This demonstrates for the first time that transcript isoforms that differ in just one or two internal exons may differ in stability. This category includes mRNAs coding for structural muscle proteins, like *Myom1*, *Myom2*, *Itga7*, *Tnnt3* and *Atp1a2*. Table 3.Location of probe sets with aberrant slopes as compared with other probe sets in the geneNumber of probe setsProbe set locationGenes145′-UTR, alternative 5′-UTR, or coding part of first exon*Tbrg4, Foxm1, Casq1, Trim37, Ly75, Casq2, Ezh2, Atp2a1, Ncapd3, Serpine1, Osbp2, Mamdc2, Bub1b, Mybph*103′-UTR or alternative 3′-UTR*Kif20a, Ankrd2, Myom2, Foxo6, Cp, Mcm5, Itgb1, Cacna1s, Pde4dip, Xirp1*9Internal exon subject to alternative splicing*Adssl1, Myom1*[^a^](#TF3){ref-type="table-fn"}*, Itga7*[^a^](#TF3){ref-type="table-fn"}*, Pygm, Rbm9, Sec24c, Tnnt3*6Internal exon not known[^b^](#TF4){ref-type="table-fn"} to be subject to alternative splicing*Tcfap4, Atp1a2*[^a^](#TF3){ref-type="table-fn"}*, Cabc1, Myom2, Mmp15*[^3][^4]

The decay of the transcript isoforms of two mRNAs was investigated in more detail: *Pde4dip* representing one of the mRNAs with differential use of 3′-UTRs, and *Itga7*, representing one of the mRNAs subject to alternative splicing of internal exons.

### Pde4dip

*Pde4dip* codes for the protein myomegalin, a protein associated with the Z-discs of skeletal muscles ([@B21]). In the ENSEMBL database, there are three annotated transcripts for this gene: a long transcript (indicated as *Pde4dip-long*) and a shorter transcript with an alternative 3′-UTR (indicated as *Pde4dip-short*) ([Figure 5](#F5){ref-type="fig"}A). The expression levels of all exons are higher in differentiated than in proliferating cells ([Figure 5](#F5){ref-type="fig"}A, note that exon expression levels reflect the cumulative expression of all transcripts containing that exon). Quantitative PCR analysis with primers specific for the 3′-UTRs of the long and short isoforms demonstrated that the short isoforms are expressed at similar levels in differentiated and proliferating cells, but that the long isoform is 25-fold higher expressed in differentiated compared with proliferating cells ([Figure 5](#F5){ref-type="fig"}C). [Figure 5](#F5){ref-type="fig"}B shows that the decay of all exons except for those unique to the short isoforms is slower in differentiated than in proliferating cells. In addition, the exons unique to the long isoform are significantly more stable in differentiating than in proliferating cells. qPCR experiments ([Figure 5](#F5){ref-type="fig"}D) validate the longer half life of the long isoform in differentiated cells. Thus, the ratio between the long and the short isoform is regulated, at least in part, at the level of mRNA degradation, and depends on the differentiation status of the cells. Figure 5.Association between RNA degradation and alternative splicing in *Pde4dip*. (**A**) Exon structure of the *Pde4dip-short* and the *Pde4dip-long* mRNA and the natural logarithm of the ratio (*y*-axis) of the expression levels of the different exons (*x*-axis) between differentiated and proliferating cells, as assayed before actinomycin D treatment. Note: using CAGE-seq (cap analysis gene of expression analyzed by high-throughput sequencing, investigating transcription start sites) data, we found that the annotation of the first exon from ENSMUST00000090750 in ENSEMBL is incorrect ([@B20]), explaining its aberrant behavior. (**B**) Plot of the decay rates in differentiated (closed symbols) and proliferating (open symbols) cells for the different probe sets in the *Pde4dip* gene, ordered from 5′ to the 3′ end of the gene. Less stable parts of the RNA have more negative slopes. (**C**) Relative abundance (arbitrary units; a.b.u) of the short and long transcripts of the *Pde4dip* gene, as assayed by qPCR with primers in the two different 3′-UTRs. (**D**) Estimated half lives of the short and long transcripts of the *Pde4dip* gene, as assayed by qPCR with primers in the two different 3′-UTRs. Error bars reflect standard error of the mean.

### Itga7

*Itga7* is an example of a gene that codes for different transcript isoforms that share the same 3′-UTR but differ in stability. *Itga7* codes for integrin α7, a transmembrane receptor important for maintaining the link between the extracellular matrix and the cytoskeleton of the muscle fiber ([@B22],[@B23]). *Itga7* is known to undergo alternative splicing, resulting in differences in both the extracellular (*N*-terminal) and cytoplasmic (*C*-terminal) domain. The splice isoform containing an additional exon (here indicated as 23, referred to in other publications as transcript α7A) is specific for skeletal muscle and expressed mainly during myogenic differentiation ([@B24]). Since this exon is out-of-frame, inclusion of this exon gives rise to a protein with an alternative *C*-terminal amino acid sequence, which may activate alternative signaling cascades. By 3′-RACE, with a primer in exon 22, we demonstrated that the transcripts with or without exon 23 indeed share the same 3′-UTR and use the same polyadenylation site ([Figure 6](#F6){ref-type="fig"}A). As expected, transcripts containing exon 23 are almost exclusively present in differentiated cells ([Figure 6](#F6){ref-type="fig"}B and D). Both exon-array ([Figure 6](#F6){ref-type="fig"}C) and qPCR experiments ([Figure 6](#F6){ref-type="fig"}E) with primers in exon 23 indicate that transcripts containing exon 23 are very unstable in proliferating cells. Moreover, analysis of the nuclear fraction of the mRNA shows that the splice variant that includes exon 23 is made in proliferating cells and is more abundant in the nucleus than in the cytoplasm ([Figure 6](#F6){ref-type="fig"}F and G), and therefore preferentially degraded in the cytoplasm of the proliferating cells. This example shows that the stability of splice variants is actively controlled and dependent on the differentiation status. Figure 6.Association between RNA degradation and alternative splicing in *Itga7*. (**A**) Agarose gel displaying the products of 3′-RACE for *Itga7* with forward primer in exon 22 in proliferating (P) and differentiated (D) cells and a 200-bp size marker (M). The longer fragment (more abundant in differentiated cells) contains exon 23, whereas the shorter fragment lacks exon 23. The size difference equals the size of exon 23 (113 bp). Subsequent sequencing confirmed this and mapped the polyadenylation site for both fragments at the previously identified and annotated polyA site (position 4018 in GenBank entry NM_008398.2). (**B**) Exon structure of three different transcript isoforms (indicated by their Ensembl transcript identifier) and the natural logarithm of the ratio (*y*-axis) of the expression levels of the different exons (*x*-axis) in the *Itga7* gene, as assayed before actinomycin D treatment. The arrow indicates exon 23, which is preferentially incorporated in transcripts expressed in differentiated myoblasts. (**C**) Plot of the decay rates in differentiated (closed symbols) and proliferating (open symbols) cells for the different probe sets in the *Itga7*gene, ordered from 5′ to the 3′ end of the gene. Less stable parts of the RNA have more negative slopes. The arrow indicates the probe in exon 23. The transcript containing exon 23 is much more stable in differentiated than in proliferating myoblasts. (**D**) Relative abundance (arbitrary units; a.b.u) of the different transcripts of the *Itga7* gene, as assayed by qPCR with primers in exon 23 (assaying only transcripts that include exon 23), exons 22 and 24 (a combination of transcripts that do and do not include exon 23), and exon 24 (assaying all transcripts), in proliferating (open bars) and differentiated (closed bars) cells. (**E**) Estimated half lives of the different transcripts of the *Itga7* gene, as assayed by qPCR with the primers described under (D) in proliferating (open bars) and differentiated (closed bars) cells. Error bars reflect standard error of the mean. (**F**) Melting curves obtained after qPCR with with forward primer in exon 22 and reverse primer in exon 24 of the *Itga7* gene in total (left) and nuclear (right) RNA of proliferating (top) and differentiated (bottom) cells. Two different PCR products were obtained. The product containing exon 23 melts at a higher temperature than the product without exon 23. The transcript containing exon 23 is present in the nuclei of proliferating myoblasts but barely detectable in total RNA (mainly cytoplasmatic). (**G**) Abundance of the *Itga7* transcript containing exon 23 in total and nuclear RNA of proliferating (left) and differentiated (right cells), as determined by qPCR with forward primer in exon 22 and reverse primer in exon 23. Note the different scales for the left and right graphs.

DISCUSSION
==========

We reported around 7000 mRNA half lives in proliferating and differentiated myoblast cells. Given the 480 min time interval analyzed, the estimates for transcripts with half lives shorter than 480 min are more precise than those for very stable transcripts. The reported half lives can guide the choice of evaluation time points in future studies where transcript levels are induced or repressed \[with the obvious exception of small interfering RNA (siRNA) experiments influencing mRNA degradation itself\]. Knowledge on mRNA half life can be used in the analysis of gene expression profiling studies to identify genes that are in fact completely shut down, whereas the data as such only indicate that expression is slightly reduced. In addition, the data can be used to flag transcripts or transcript isoforms that can not be valued since their mRNA half lives exceed the time points analyzed, or demand special attention due to significantly altered half lives.

Furthermore, we show that transcripts, including specific splice variants, with a role in myogenesis and differential expression between proliferating and differentiated myoblasts frequently demonstrate different stabilities in the two states. Transcripts with higher expression in differentiated than in proliferating cells are more rapidly degraded in proliferating cells and vice versa. A cell is in some cases reminiscent of a Dutch auction. When the clock is down to zero and there appears to be no demand, the commodities are to be destroyed.

Despite the characterization of a multitude of proteins and pathways involved in mRNA degradation ([@B9]), we are only beginning to understand the mechanisms that control the stability and degradation of individual transcripts. The 3′-UTR is very important for transcript stability given the presence of miRNA target sites, AU-rich elements (AREs) and binding sites for other RNA-binding proteins. Our observation that decay constants of first or last probe sets in a gene frequently differ (are higher or lower) from the other probe sets is consistent with the reported extensive use of differential 3′-UTRs (and alternative promoters) depending on the tissue or cellular state ([@B10],[@B25; @B26; @B27]). In differentiating T-cells ([@B10]) and in cancer cells ([@B27]), differential expression of transcripts has been related to alternative 3′-UTR (polyadenylation site) use, resulting in differences in mRNA stability.

It is not yet clear which RNA binding proteins are responsible for the control of mRNA stability during myogenic differentiation. The HuR and KSRP can transiently affect the stability of transcripts containing AREs during myogenic differentiation ([@B17; @B18; @B19],[@B28]). However, only a minority of the transcripts with variable decay rates contain AU-rich elements \[ARE motifs, predicted by the ARED algorithm ([@B29])\] in their 3′-UTR. From our digital gene expression studies ([@B20]), two members of the tristetraprolin family of RNA binding protein emerge as other important candidate factors for differentiation state-dependent control of mRNA degradation: ZFP36L1, which has been shown to accelerate degradation of tumor necrosis factor-alpha (TNF-α), vascular endothelial growth factor (VEGF) and granulocyte macrophage-colony stimulating factor (GM-CSF) ([@B30; @B31; @B32]) and ZFP36L2. The former is 10-fold increased and 3-fold decreased in expression during myogenic differentiation. The muscle-specific RNA binding protein CUGBP has been shown to stabilize mRNAs by binding to GU-rich elements in the 3′-UTR ([@B8],[@B33]). However, CUGBP shows little variation in expression between proliferating and differentiated cells and is therefore less likely to be responsible for the differences in decay rates observed in proliferating and differentiated myoblasts. In addition to RNA stabilizing and destabilizing proteins, also miRNAs and components of the nonsense-mediated decay machinery have been shown to have tissue-dependent expression levels and activities and may contribute to the observed differences in mRNA stability ([@B34; @B35; @B36; @B37; @B38]).

One of the new aspects of the current study is the notion that also splice variants that share the same 3′-UTR but include alternative internal exons may differ in stability. It was previously thought that the splicing machinery determines the abundance of the different splice isoforms aided by tissue-specific splicing factors like FOX, CELF, PTB, MBNL, hnRNP A/B ([@B26],[@B39],[@B40]). From our study, it appears that the presence of stabilizing or destabilizing elements in the alternative exons or differences in structural characteristics of the synthesized transcripts may also contribute to their differential expression. We show this in detail for the *Itga7* transcripts variants with or without exon 23. The transcript including exon 23 appears to be me far more stable in differentiated cells, whereas in proliferating cells it is subject to fast degradation. A comparison of the nuclear fraction of the RNA with total RNA (nuclear and cytoplasmic) shows that the degradation of the transcript containing exon 23 mainly occurs in the cytoplasm. This observation is in line with a recent publication showing major differences between transcript variants in their nucleocytoplasmic distribution, at least partly attributable to mechanistic differences between nuclear and cytoplasmic mRNA degradation ([@B41]). We investigated whether exon 23 contains any known destabilizing elements such as AREs or PUFs ([@B7]), but could not find any. It is still possible that exon 23 contains an miRNA binding site but prediction programs failed to predict a good candidate. Thus, the mechanism for its preferential cytoplasmic degradation in proliferating myoblast remains elusive.

We conclude that mRNA degradation is a key controlling factor in gene expression regulation, possibly as important as transcription factor-mediated induction of mRNA synthesis. The mechanisms of mRNA degradation are currently less well understood than the mechanisms behind transcription factor activity regulation. We have gained several mechanistic insights in the process of mRNA degradation, but many mechanistic aspects deserve more attention in future studies. Knowledge on mRNA degradation rates, as presented in this study, is important for the proper interpretation of many biological studies. Decay rates need to be taken into account when choosing a time point for evaluation of the perturbation of a biological system, aid the modeling of transcriptional networks, and should be incorporated as important parameters in systems biology approaches.

ACCESSION NUMBER
================

Gene Expression Omnibus GSE14387.

SUPPLEMENTARY DATA
==================

[Supplementary Data](http://nar.oxfordjournals.org/cgi/content/full/gkq790/DC1) are available at NAR Online.

###### Supplementary Data

FUNDING
=======

Funding for open access charge: Centre for Medical Systems Biology within the framework of The Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO).

*Conflict of interest statement*. None declared.

[^1]: Gene annotation was derived from Gene Ontology. The biological processes significantly overrepresented in the 100 fastest and slowest decaying transcripts in proliferating and differentiated cells are displayed, together with the *P*-value from the hypergeometric test used for testing overrepresentation. Given multiple testing issues, *P* \< 1 × 10^−4^ was considered to be significant.

[^2]: Only genes annotated with 'muscle development' in Gene Ontology and a minimum absolute fold change of two are shown.

[^3]: ^a^Two probe sets in two different exons affected.

[^4]: ^b^According to Ensembl annotation.
