A continuum model of transcriptional bursting

1Transcription of genes is discontinuous, occurring in irregular bursts or pulses of activity, interspersed by irregular intervals of inactivity. Bursting transcription is conserved in all forms of life, from prokaryotes to mammalian cells and tissues. The irregular nature of transcriptional bursting is proposed to be a major driver of spontaneous heterogeneity in gene expression, which in turn drives diversity of cell behaviour in differentiation and disease. Bursting reflects the underlying mechanisms of transcriptional regulation, and measures of bursting can reveal the dynamic processes absent from standard population average measures of RNA expression.

2The standard framework used to describe transcriptional fluctuations compares one state and two state models. In the one state model, transcription occurs with a constant probability, which for moderately and strongly transcribed genes, will generate a low variance in their total transcribed RNA per cell. In some contexts, notably budding yeast, the variance in RNA abundance measured by single molecule RNA fluorescence in situ hybridisation (smFISH) can fit this one state scenario, where the distribution of RNA per cell is well characterised by a Poisson distribution. In many other contexts, the one state model does not fit the smFISH data, with measured RNA abundance showing too much variability between cells than can be produced by a constantly active gene. To explain this increased variance, the more complex random telegraph (or two state) model is often invoked. In this model, the gene switches stochastically between an active state, where mRNA production occurs with constant probability per unit time, and an inactive state, with no mRNA production. The extra state increases the potential variability in output from cells, and can therefore predict the observed extra spread in transcript abundance in the cell population.

3Use of the two-state model in fitting smFISH and protein distributions allows estimates of the parameters of the transcriptional fluctuations, usually the burst size (number of transcripts produced in a burst) and burst frequency (the frequency with which a burst occurs). However, these dynamic properties are usually inferred from a population distribution at a single time point, assuming each cell is part of a homogeneous population with fixed values of the switching rates, transcript production rate and transcript lifetime. In other words, the perception has emerged that transcriptional bursting is a product of molecular noise, rather than a process responsive to the demands of the cell. A rethink is required, not least because of recent work demonstrating burst size and frequency are quantities that can be modulated by extracellular signals and cell properties such as volume and cell cycle stage. These studies challenge the notion, central to the standard two state model, that a population of cells consists of those where the gene of interest is 'off' and those where the gene is 'on' with a constant probability of firing.

4To make accurate models of transcriptional fluctuations and how they are regulated, it is critical to directly observe and quantify how transcription evolves over time. To directly measure features such as burst size and burst frequency requires data capture of complete sequences of bursts, rather than snapshots. Imaging transcriptional output in living cells is possible using RNA detection systems based upon the binding of a bacteriophage coat protein to stem loops of RNA. An array of sequence encoding stem loops, such as MS2 or PP7, is inserted after the promoter of a gene of interest. When the gene is transcribed, the nascent RNA stem loops are detected with a fluorescent MS2 (MCP) or PP7 (PCP) coat protein, which is constitutively co-expressed. The system is visualized using time-lapse fluorescence microscopy, with a nuclear spot indicating active transcription of the gene of interest. The intensity of the spot reflects the instantaneous nascent RNA load at the gene and changes in spot intensity reflect how the RNA load fluctuates over time. These tools allow direct observation of the dynamics of transcription regulation, revealing insights into the mechanics of Poissonian transcription and developmental regulation. Previous work has interpreted the appearance and disappearance of transcription spots in terms of ‘bursty’ transcription. Here, the exponential nature of the 'ON' and 'OFF' time durations of the transcript spot was related to the two state model, with the exponential behaviour proposed to reflect rate-limiting steps in transitions between the active and inactive states. Although this is an appealing inference, how spot fluctuations actually reflect the dynamics of the transcription machinery is unclear.

5In this paper, we test the current models for explaining transcriptional fluctuations, using a combination of live cell imaging, computational modelling and simulation, and targeted mutations of gene and promoter structure. We use a probabilistic approach to infer dynamics at the molecular level from fluctuations in spot intensity. We make quantitative measurements of the transcription site RNA abundance and the retention time of nascent RNA at the gene. We use these measurements to train candidate hidden Markov models to describe the underlying initiation of RNA polymerases, and find that a spectrum or continuum of initiation rates describes experimental data more accurately than a binary off/on model or discrete levels of activity. Finally, we investigate how the processes of transcription elongation and initiation contribute to the transitions of the gene over this spectrum of activity states.

6To monitor transcriptional dynamics in living cells, we integrated an array of MS2 stem loops after the promoter of the endogenous actin5 locus, a strongly expressed actin gene in undifferentiated Dictyostelium cells. Transcription continues after the MS2 loops into the coding sequence, then native terminator, to generate a full length transcript of around 2.5 kb. We visualized the resulting transcription dynamics using time-lapse fluorescence microscopy and extracted time series of spot intensities using custom-built software integrating both cell tracking and spot detection. The movie sequence in Figure 1A illustrates the tracking of a typical cell, showing the fluctuations in spot intensity over time. Figure 1B shows the measured transcription spot intensity for the cell in Figure 1A, with a kymograph of the spot fluctuations. For most genes studied, the durations for which a spot is present or absent are often measured to have approximately exponential distributions, which is the case for act5 (Figure 1—figure supplement 1). Exponential timescales have been inferred to represent modulation of gene activity, between the ON and OFF states of the two-state model, with a rate-limiting step determining switching between states. A more simple possibility is that stochastic fluctuations of a gene with no OFF state - the one state model - could give rise to pulses and intervals between pulses with the experimentally observed lifetimes.

7To test if a one state model, where the probability of a polymerase initiating is constant over time, can give rise to the distributions of spot durations observed experimentally, we constructed a Monte-Carlo simulation framework describing the MS2 system, incorporating polymerase initiation (the rate at which polymerases begin transcription) and elongation (the velocity of the polymerase while transcribing) (Figure 1C,D). Additional features, such as the number of gene states, and termination/release rate (the rate at which the RNA leaves the gene) were omitted at this stage. By defining a signal detection threshold above which a spot is determined to be present, we calculated the pulse and interval length distributions from the simulated data. The average pulse duration depends strongly on the initiation rate, but was also substantially sensitive to the detection threshold and frame interval specified in the simulations. The pulse duration distribution had exponential distributions at long timescales, with deviations at very short timescales (Figure 1E). This was surprising, as the exponential distribution was previously inferred to support the two-state model. Here, our simulations clearly show such a distribution can arise from a more simple one-state scenario.

8An exponential distribution of pulse durations has one parameter, a characteristic timescale, and arises from there being a constant probability of the pulse ending at any given time. This probability is independent of the duration for which the pulse has been sustained. It is tempting to infer from this probability the rate of switching from an active to an inactive state of a two state gene. However, with exponential behaviour observed for our simulation of constant polymerase initiation rates (this is only one state), it is necessary to consider additional contributions to the likelihood of a spot disappearing from one imaging frame to the next.

9To investigate the origins of exponential behaviour, we constructed a theory of spot appearance and disappearance for the live MS2 system. By considering the rate of initiation of polymerases at the start of the gene and the number of RNA stem loops produced by a polymerase as it moves along the gene, we calculate the joint probability distribution linking the spot intensity in the previous frame with the current frame. (Figure 1F, full details in the appendix). The theory successfully predicts the characteristic exponential timescale of pulse duration in terms of the one state model, incorporating the polymerase elongation rate, the frame interval and threshold of spot detection. In simple terms, the pulse duration reflects the balance between polymerase initiation and elongation or termination - in order for a spot to remain present in the next frame, a sufficient number of new polymerases must be initiated to replace the loss of spot intensity from cleaved RNA leaving the transcription site. Intuitively, high initiation rates are likely to give rise to large numbers of initiations per frame, giving long pulses. Increasing the threshold of spot detection decreases the likelihood of sufficient initiations and thus decreases the pulse duration. These predictions are in agreement with the Monte Carlo simulations of the MS2 system (Figure 1G and Figure 1—figure supplement 2). The initiation and elongation rates have approximately inverse effects, with initiation scaling directly with spot intensity (as expected intuitively) and elongation scaling inversely. The effect of slow elongation increasing spot intensity can be thought of in terms of a polymerase 'traffic jam', with several polymerases building up behind a slower one, causing a build up of nascent RNA on the gene. Additionally, the frame interval of observation influences the effective pulse duration, as longer frame intervals increase the likelihood of a gap between pulses being missed between frames. Overall, this theory shows that the length of time for which a spot is visible is not simply related to switches between proposed gene states in two state models and explains the dependence of the pulse duration on the imaging signal-to-noise and frame interval.

10Since the theory shows that the simple presence of a spot does not provide insight into what activity state the gene is currently in, we must make precise measurements of transcription site intensity in order to detect initiation rate changes (changes in gene state). As depicted by the cartoon in Figure 2A, a fluctuation in intensity might be consistent with either the one state (Poisson, top) or multi-state (bottom) models. Without calibration of the intensity fluctuations in terms of numbers of RNAs, it is difficult to discount either model. One strategy that can be used is an autocorrelation analysis, which measures the magnitude and timescale of intensity fluctuations. The autocorrelation measures the ratio of intensity fluctuations relative to the mean intensity. For a Poisson (one state) gene, where initiations behave independently of the previous or subsequent polymerase, this ratio is related to the average number of bound polymerases, allowing the initiation rate to be estimated. This was previously used to infer the initiation and elongation rates from a POL1 promoter inserted upstream of the GLT1 gene in budding yeast, specifying the time for the correlation to decay to zero as the time for which a single polymerase contributes to the spot intensity.

11Applying a similar analysis to transcription of the Dictyostelium act5 gene (Figure 2B), we find that inference of polymerase properties is complicated by an additional long-timescale decay of several minutes to tens of minutes (black dotted line). The additional slow decay means the RNA level at the gene is correlated beyond the dwell time at the gene of a single RNA. This argues against a Poisson model and indicates a model with two or more transcription states would better describe act5 transcription. If the autocorrelation is recalculated using only frames where a spot is present, the slow decay remains, suggesting that the active periods where a spot is present are not defined by a single constant initiation rate; in other words there are multiple 'ON' states. By subtracting a linear estimate of the long timescale decay from the autocorrelation, we calculated estimates of the RNA dwell time and polymerase load on the gene. The point where the autocorrelation deviates from the long decay gives an approximation of the dwell time as 190 ± 20 s. The average number of polymerases contributing to the intensity is estimated from the magnitude of the deviation at zero lag time as 20–27 polymerases. The cells are predominantly in G2, so these autocorrelation-based estimates are for the polymerase load across 2 alleles of the act5 gene comprising the transcription spot.

12As described above, the observed pulsing properties also depend on the RNA dwell time- the time the RNA spends at the transcription site during elongation and termination before it is released from the gene. As an independent measure of dwell time, we used fluorescence recovery after photobleaching (FRAP) measurements on MS2 spots. The MCP-RNA interaction is stable, so recovery of the spot intensity after bleaching is determined by transcription of new MS2 loops. In the absence of photobleaching, spot intensities fluctuate over time, therefore a bleached spot is not expected to recover to its pre-bleach intensity. We therefore normalized the curves based on the recovered intensity. Recovery of a single endogenous gene is highly stochastic, as expected from simulations of photobleaching and recovery (see below). Therefore, we averaged over 20–30 cells in order to extract the typical RNA dwell time. The average recovery curve was fit to a processive polymerase model of the intensity recovery to estimate the elongation rate and termination time (full details in the appendix). We calculated the dwell time for three different gene constructs, with the MS2 loops integrated into the 5’ gene region, the 3’ region, and replacing the entire actin5 coding region (Figure 2C). The dwell time for the 5’ insertion was longer than for the other insertions, as expected as the polymerase will continue to transcribe after the MS2 array, whereas with the other two cell lines, termination will occur immediately upon completion of synthesis of the array. Using the difference in dwell times and assuming a constant rate of elongation across the gene, we calculate the elongation rate as 22 nt.s-1 and termination time as 60-70 s. This gives a total dwell time of around 170 s, for the 5’ MS2 insertion, roughly in agreement with the autocorrelation-derived approximation on the same allele.

13As shown by our pulsing theory, spot intensity variations may reflect either fluctuations in the initiation rate or dwell time. If the dwell time is systematically different for spots of different intensity, we would observe this by studying the recovery times of high and low intensity spots. We divided the FRAP data into three equal groups based on their recovery intensity (inset of Figure 2D) and calculated the average recovery curve for each group. Figure 2D shows that no significant differences in recovery time are observed for the different intensity groups, suggesting dwell time variations are not the primary source of the spot intensity variability.

14We then calibrated live RNA spot intensity in terms of the number of MS2 stem loops by comparison with smFISH using a probe for the MS2 RNA. We used FISH-quant to estimate the intensity of one mature cytoplasmic RNA then used this information to calibrate the number of nascent MS2 loops at the transcription site (Figure 2E). We compared the distribution of these nascent RNA counts with the intensity distribution of transcription site MCP-GFP intensity measurements from live cells immediately prior to fixation for smFISH (Figure 2F). The two distributions are aligned by calculating the percentiles of each distribution and using these values as calibration points (Figure 2G). Due to the differing precision with which the intensity can be calculated in live and fixed samples, there is some deviation from a linear relationship at the extreme ends of the calibration curve, and a kink at low intensities where the detection threshold of live spots is higher than in fixed spots, however the overall trend can be used to calibrate the spot intensities measured in live cells in terms of the number of MS2 stem loops present.

15After calibrating the system, we used our pulsing theory outlined above to estimate the initiation rate. We measured how the mean pulse duration changes as the detection threshold is varied. For a one-state (constant activity) gene, the pulse duration would lie on a contour calculated using our theory. Figure 2—figure supplement 1 shows that the experimental data does not lie on a single contour of constant initiation rate, consistent with the autocorrelation results in Figure 2B. Similarly, contours calculated from a simple two-state implementation, with on- and off-rates chosen to produce the experimentally observed spot frequency (details in the appendix), cannot match the variation of pulse duration and spot frequency as the detection threshold is varied (Figure 2—figure supplement 2). Instead, the behaviour at high thresholds suggests a short-lived higher initiation rate.

16Dictyostelium cells are almost exclusively in G2, owing to a very short S phase and complete absence of G1. Therefore almost all cells will have replicated the MS2-tagged gene and have two genes, held together by sister chromatid cohesion, contributing to a single resolved spot. Transcription of the unreplicated gene will not confuse analysis, as act5 transcription is not robustly detected in the first 15 min after mitosis, when euchromatin is replicated. A higher initiation rate state is possible in the situation where both copies are simultaneously active - giving a potential third state with double the initiation rate of a single copy.

17To test the possibility of this third state, we applied probabilistic modelling to our quantitative data to assess the extent to which three state or higher models can describe the observed transcription dynamics. Since a two-state model of transcriptional activity, with an ‘off’ state and an ‘on’ state of constant initiation rate, does not describe the live data adequately, we asked whether an improvement could be made by including additional active states with different initiation rates. It is not possible to directly observe polymerase initiations using the MS2 system, so we used hidden Markov modelling to infer the likelihoods of different numbers of active states from the measured live cell intensity data.

18In simple terms, a hidden Markov model (HMM) describes a sequence of observations in terms of a sequence of underlying, 'hidden' states. By extension from the two-state model, we define a 'state' as a configuration with a single constant (Poisson) initiation rate. In terms of transcription dynamics, the hidden state corresponds to the underlying initiation rate of the system. This rate cannot be observed directly and instead must be inferred from measurements of the transcription site intensity, which is related probabilistically to the initiation rate. Using standard HMM techniques, the magnitudes and dynamics of the underlying initiation rate can be optimized to maximize the probability of generating the experimental data from the model. For the case of the MS2 system, the spot intensity does not depend only on the instantaneous initiation rate, because polymerases initiated in previous frames are still present on the gene and also contribute to the spot intensity. To take account of these polymerase contributions, we constructed a model (full details in the appendix, architecture depicted in Figure 3A) with hidden states representing the initiation rate (blue circles) and the number of polymerases initiated (orange circles), inferred from the sequence of transcription site intensities (green squares). A benefit of our two-layer HMM framework is that we explicity model the polymerases uncoupled from transitions between states of gene activity. This means that the distribution of polymerases initiated in each state can be estimated (see below).

19To test how well models with discrete states of activity describe the data, we calculate the probability of obtaining the experimental data using the model parameters. This probability is known as the likelihood. We use a standard approach to compare models of different complexities, by calculating the Akaike Information Criterion (AIC) from the likelihood. The AIC is used to find the best approximation in situations where the real scenario is likely to be highly complex. Models are penalized by their number of free parameters, preventing overfitting by excessively complex models. The AIC estimates the prediction error of the candidate model, giving a lower numerical value for models which more accurately describe the experimental data.

20To test the validity of the use of the HMMs and associated algorithms, we tested the framework on simulated data, with known numbers of promoter states, and known timescales of switching between states (Figure 3B). In these cases, the 'hidden' initiation states and numbers of initiated polymerases are known, allowing the accuracy of the probabilistic methods to be measured. Applying the hidden Markov approach to a simulated 3 state model accurately estimated 3 discrete initiation states as the point at which AIC is minimized (Figure 3C). 1 and 2 states give a higher AIC (lower likelihood), whilst 4 shows no improvement over 3 states. In addition, this approach robustly measured the initiation rate of each state, as well as the timescale of switching between the activity states (Figure 3—figure supplement 1).

21We next applied the model fitting to experimental data for transcription of the actin5 gene. We carried out 4 independent time lapse experiments, with multiple fields of view in each experiment, capturing data on 44–145 cells per experiment, giving 1686–6350 individual cell-frames, per experiment. Figure 3D shows the AIC for models with discrete numbers of initiation rate states. The 2-state model shows a strong reduction in AIC compared to a 1-state model, with the 3-state model better still. The curve does not plateau after the 3 state model, instead, the fit to the data is continually and significantly improved by adding further initiation states. The continually improving fit does not depend on the permitted mechanisms of switching allowed from state to state (Figure 3D). This continual improvement in fit occurs if the gene constrained to switching between states one step at a time (Figure 3D- 'ladder'), in steps of up to two states ('Jump 1'), or is unconstrained ('Free'). The most realistic interpretation of the continually improving fit with more states is that the gene switches in activity over a spectrum or continuum of states, rather than a small number of discrete activity states.

22In Figure 3D the largest improvement in fitting occurs up to 3 states; therefore, is a model with three initiation rate states appropriate for the act5 gene? Although the slow rate of improvement in fit for higher numbers of states may suggest three states are adequate, closer analysis of the optimal 3-state fit reveals its inability to capture the full variability of the transcription site intensities. The values of the estimated initiation rates are not uniformly spaced (as would be expected for two alleles undergoing random telegraph transcription), but have a low state and a very high state in addition to an off-state. Furthermore, Figure 3E shows strongly non-Poisson distributions for the predicted number of polymerases initiated per frame in each state of the three state model, inconsistent with a constant initiation rate within each state. In finding the best 3 state fit, the modelling process effectively forced outlying data into non-Poissonian states. The wide polymerase distributions imply either that the initiation rate is varying within a state, or that transitions are occurring between states on a timescale significantly faster than the frame interval.

23To determine the dynamics of the gene state and number of initiated polymerases, we used a Gibbs sampling algorithm (see appendix). The estimated number of polymerases initiated in each frame is approximately independent of the number of states chosen for the model. The initiation rate, the number of events per unit of time, is difficult to define as an instantaneous measurement; instead, by plotting the cumulative number of polymerase initiations as a function of time the initiation rate can be taken from the gradient. For a typical cell in the experimental data, the cumulative polymerase plot is composed of linear segments (Figure 3—figure supplement 2) representing periods of time over which the initiation rate is constant. We measured the initiation rate using an edge-preserving smoothing filter and found a diversity of such gradients, indicating a spectrum of initiation rates, again implying the transcriptional behaviour is not adequately described by a few discrete levels of activity. The most realistic model to account for this additional complexity in initiation fluctuations is a ladder or continuum of activity states. The timescale of initiation rate variation through this spectrum of activity states can be determined using an autocorrelation analysis of non-zero initiation rates. Figure 3F shows such a plot, which reveals the initiation rate fluctuates with an average timescale of around 5–6 min, although it is possible for a roughly constant initiation rate to be sustained for up to 15 min (Figure 3—figure supplement 2).

24To what extent can fluctuations in elongation rate contribute to the complexity of the ladder or continuum of transcriptional states? A transiently slower rate of elongation may cause a build-up in the number of polymerases on the gene and therefore contribute to non-Poisson variations in spot intensity. To test this possibility, we use simulations with three initiation rate states – with average initiation and elongation rates matched to the experimental estimates – to address whether the additional complexity shown experimentally beyond a three state model (Figure 3D) could be accounted for by adding elongation rate fluctuations to the system.

25In addition to increasing the spot intensity, slower elongation rates lead to increased dwell time, which in FRAP measurements, would lead to the brighter spots taking longer to recover. However, the experimental data (Figure 2D) showed no observable difference in recovery between groups of different intensity. It remains possible that fluctuations in elongation rate may occur over a timescale which would not be resolved by FRAP. To investigate the timescale of dwell time fluctuations that might be invisible to FRAP, we developed simulations of the FRAP protocol, matching the experimental procedure as closely as possible. We incorporated temporal fluctuations in the polymerase elongation rate acting either globally (affecting every polymerase on the gene in the same way) or in a polymerase-by-polymerase manner. Since very little is known about the type of elongation fluctuations possible in vivo, we implemented a system which switches randomly between 10 and 30 nucleotides/s, and varied the timescale of fluctuation (how long the system remains in either state).

26The FRAP curves produced by individual simulated cells are treated in the same way as experimental data - dividing the cells into three groups based on their recovery intensity and rejecting cells showing no recovery – to determine the largest fluctuation timescale which shows no difference between intensity groups. As shown in Figure 4A, elongation rates fluctuating independently for each polymerase cannot produce any differences in the recovery curves. For the density of polymerases on the gene for typical spot intensities, fluctuations of individual polymerases are suppressed by catching up and being blocked by a polymerase immediately downstream. This results in a bulk polymerase elongation rate close to the lower speed of 10 nucleotides/s. Thus independent fluctuations in polymerase elongation cannot account for the fluctuations in spot intensity in the continuum model.

27In contrast, global fluctuations in the elongation rate give rise to different recovery curves for different spot intensities (Figure 4B), as the brightest spots tend to be caused by slow-moving polymerases staying at the transcription site for longer, while the intermediate spots have a bias towards faster moving polymerases. As in experiments, the lowest intensity recoveries show a high degree of noise, due to the small normalisation factor. As the timescale of the elongation fluctuations becomes shorter, the difference between the recoveries is reduced, as variability in the elongation rate is blurred out by faster transitions. The point at which the simulated data becomes consistent with the experimental traces (Figure 2D) is around the lower two panels of Figure 4B, that is, an elongation fluctuation timescale with an upper limit of around 10-40 s. What this effectively means is that gene-wide fluctuations in the elongation rate may be present on a timescale of up to around 40 s without being apparent in experimental FRAP recovery curves. The exact timescale depends also on the magnitude of the fluctuations, taken here to be around 50% of the mean elongation rate, nevertheless, some degree of elongation fluctuation may contribute to the temporal variation of the transcription spot intensity.

28We then asked whether adding elongation fluctuations of the magnitude identified above to the three-state model is sufficient to produce the decrease in AIC value beyond 3 states seen in hidden Markov modelling. The difference, ΔAIC, between AIC (3-state) and the minimal AIC value (for a model up to 8-states) represents the additional complexity present in the data – essentially the information missing from a three-state fit. We applied our HMM fitting to simulations of the three-state model with different timescales of elongation fluctuation, and measured ΔAIC to estimate if the addition of elongation fluctuations is sufficient to recapitulate the additional complexity measured experimentally.

29Figure 4C compares the ΔAIC value between a representative experiment and simulations of elongation fluctuations (adjusted to match experimental sample size). For no fluctuations, no improvement is seen as expected, since a three-state model successfully describes the simulated data. For 13 s elongation fluctuations a small increase in complexity is observed. For the 43 s timescale, the magnitude of the improvement is comparable to that observed experimentally. This suggests these elongation fluctuations can potentially recapitulate the variability in spot intensities observed in experiments, with two important caveats. Firstly, in the simulations the 43 s timescale data plateau at 5 states, rather than around 7–8 for the experimental data. Secondly, training a three-state model on experimental data yielded distributions of polymerase initiations in each state with variances greater than expected for Poisson distributions (Figure 3E). This super-Poissonian behaviour of the experimental data is not reproduced by training the same three state model on in silico elongation fluctuation data (Figure 4D, Figure 4—figure supplement 1), indicating further complexity in the experimental data not fully explained by the simple elongation fluctuations.

30To test the effects of perturbing the initiation rate on transcriptional bursting, we generated point mutations in the TATA box of the act5 promoter, T1A and A2C (Figure 5A and Figure 5—figure supplement 1 and 2), which have strong effects on transcriptional output in yeast genes and mammalian expression plasmids. Both TATA mutant lines displayed a slight reduction in overall spot intensity (Figure 5B), with FRAP experiments suggesting a similar dwell time to wild type (Figure 5—figure supplement 3). Analysis of time-lapse experiments using the polymerase HMM framework found that the TATA mutants spend less time in medium and high initiation rate states, and more time at lower initiation rates, compared to control cells (Figure 5C). The overall timescale of initiation rate variability, measured through autocorrelation, was not substantially changed (Figure 5D) with both the wild-type and TATA mutations showing fluctuation timescales of several minutes. In addition, we observed no clear difference in switching rates between the inactive and active states. The rate of switching to the inactive state, k(off), was unchanged between wild-type and TATA mutants (Figure 5E). The tendency of spots to appear, k(on), showed a slight impairment in A2C mutants, although this was not statistically significant. This subtle effect might be also interpreted as the enhanced occupancy of active states of unobservably low intensity in the A2C mutation, rather than simply the absence of transcription. We then addressed the likelihood of the gene switching up or down in initiation rate based upon its current state (Figure 5—figure supplement 4). For all cell lines, both wild-type and TATA mutant, high initiation rates had a tendency to revert to lower initiation rates, and lower initiation rates had a tendency to revert to higher initiation rates. In the TATA mutants, the initiation rate is less likely to switch to a higher activity state, resulting in reduced time spent in high activity states. Together these observations suggest that perturbing the TATA box does not affect the duration or frequency of active states, but rather modulates the initiation rates that are possible.

31The standard analysis of transcriptional bursts using smFISH measurements can extract estimates of the burst size and burst frequency. Implicit in the assumptions of the random telegraph (two state) model is the idea that cells undergoing a burst in a particular gene all have the same underlying initiation rate. However, it seems unlikely, based on our data, that the active initiation rate will be constant over time during a burst, or that it will be the same for all cells in a population. Rather, the rate of initiation will depend on the binding of molecular factors, which turn over on the timescale of seconds to tens of seconds, based upon measured residence times of transcription factors. The cytoplasmic lifetime of mRNA means that these fast fluctuations are blurred out in static measurements, such as smFISH. In most such cases only broad differences in overall transcript content between cells can be determined, which accounts for the success of models with a few numbers of discrete states in fitting the data. Measurements using the MS2 system are integrated over the dwell time of the nascent RNA, a timescale of 2–3 min. The dynamic measurements can begin to resolve systematic variability in the initiation rate over time, which is difficult to assign to a small number of discrete states. Our proposed continuum model is summarized in Figure 6A. This model may also be approximated by a discrete state model with more states than can be effectively detected (Figure 6—figure supplement 1).

32Our continuum model may be consistent with a recent theoretical model which proposed rapid switching between states as an improvement over the standard 2-state model in terms of information transmission (Figure 6—figure supplement 1iii). In this view the initiation rate fluctuations can be thought of as the fraction of time for which the factors required for initiation are bound to the promoter. The rate of binding depends on the spatiotemporal variability in local concentrations, producing a continuum of activities when integrated over the MS2 RNA dwell time. When the TATA box is mutated, the overall decreased initiation rate in the active state may reflect the reduced rate of binding of necessary transcription factors, and therefore reduced occupancy of the initiating state (Figure 6—figure supplement 2). Illustrating the potential equivalence of this rapid switching model with our continuum model, our data show that the optimal discrete 3-state model consigns polymerase initiations to non-Poissonian states (Figure 3E), which implies either multiple initiation rates within a state (in other words more than 3 states) or switching between the discrete states faster than the imaging frame interval.

33Is this continuum view in conflict with the existing smFISH data on 2-state models? To answer this question, we simulated steady state levels of RNA derived from a gene where the binding rate varies continuously over time, overlaid with longer periods of inactivity (Figure 6B, upper). We incorporated autocorrelation-derived timescales and initiation rates extracted from the hidden Markov modelling together with a mRNA lifetime of 40 min, similar to measured lifetime values of act5-MS2 using actinomycinD-treated cells. The distribution of cytoplasmic mRNA molecules was well fit to a standard two-state bursting model, using a negative binomial function (see appendix) capturing primarily the long timescale activity intervals (Figure 6C). The shorter timescale of the MS2 measurements (Figure 6B, lower) is able to discriminate fluctuations within active regions. In other words, when measured in the steady state, using methods such as RNA FISH, a continuum of transcriptional states will appear to consist of far fewer states than can be detected using live cell approaches, highlighting the different timescales measured by different methods. The two-state model approximates the magnitude and broad timescale of transcriptional variability, although it should be noted that the bursting parameters extracted may not directly correspond to the physical quantities of burst size and frequency.

34We have carried out in-depth measurements of the transcription dynamics of an endogenous, highly-expressed gene in Dictyostelium. This analysis has provided an insight into the dynamics of transcriptional regulation at short timescales, allowing models of the underlying mechanics to be discriminated. We interpreted fluctuations in the intensity of the transcription site in terms of the time-resolved rate of polymerase initiation, and found that the behaviour is not consistent with either constitutive, Poissonian activity or switching between a small number of discrete activity states. Instead, we observed a spectrum of activity states, characterised by an initiation rate varying over timescales of several minutes.

35This spectrum of activity may be produced by very short periods of transcription initiation, where the activity lifetime or rate of reactivation is dependent on the time-varying local activity of molecular factors, or the gene locus switching between different topological conformations. When integrated over the retention time of nascent RNA at the transcription site, this ‘microbursting’ behaviour would give rise to a continuum of initiation rates. An alternative model might involve each cell having a distinct initiation rate in the active state, thereby producing a spectrum of initiation rates but with a simple two-state model for each cell. This model is equivalent to a limiting case of our continuum model with infinitesimally slow variation in the activation rate, however, for act5 we observed timescales of variation of minutes by autocorrelation or Gibbs sampling of initiation rates. The differences between our continuum view and models with each cell having it’s own stable transcription rate may reflect differences in experimental system. Alternatively, the differences may arise because live cell analysis of nascent RNA dynamics allows the timescales of fluctuation to be directly extracted from the data, which is not possible from measurements of protein abundance or steady state transcript counts.

36Bursts of transcriptional activity have been described in a wide range of systems. To investigate the molecular origins of the bursting phenomenon, we constructed a simple theory of Poissonian transcription of the MS2 system. Unexpectedly, this simple system with constant promoter activity predicts exponential distributions of transcriptional pulses and intervals, with the pulse duration determined not only by the promoter initiation rate but also by the frame interval of observations and threshold of spot detection. We confirmed the predictions of the theory using Monte Carlo simulations of the MS2-tagged gene. In general, these simulations provide a framework for studying fluctuations in transcription site intensity, incorporating stochastic events for polymerase initiation, transition to elongation, elongation and termination/release. Additional complexity could be added in the form of the probability of promoter-proximal pausing, premature termination or interactions between adjacent polymerases. This approach allows us to predict the change in the spot intensity due to adding or changing a system parameter, and therefore determine whether competing models can be distinguished experimentally. This is illustrated by our analysis of the fluctuations in the elongation rate, where, in comparison with experimental data, we estimated an upper limit for the timescales of fluctuations at around 40 s. Analysis of elongation rate changes suggested such fluctuations could not account entirely for the complexity of the spectrum of transcriptional states, although they could potentially contribute.

37Interspersed with the active regions are periods of relative inactivity, where no transcription spot is visible. The number of steps involved in the reactivation of an inactive gene has been inferred from the distribution of the off-durations. Due to the current detection threshold of our system, a low level of basal transcription cannot be ruled out in this inactive state. Consistent with early views of bursting, the longer periods of inactivity may reflect slower dynamics of the remodelling of chromatin, rather than binding of specific transcriptional regulators, based upon our perturbation of the core promoter. Mutation of the act5 TATA box reduced the overall amount of transcription, primarily by reducing the amount of time spent at high initiation rates rather than changing the switching dynamics between the off and active states. Importantly, neither TATA box mutation we studied abrogated transcription entirely, as might be expected if the pre-initiation complex can no longer bind the promoter. A reduction in the duration of periods of high activity might suggest an impaired duration or frequency of binding, or disruption of normal promoter conformational switches in response to binding. The strong effects of TATA mutations in yeast and plasmid systems may reflect more simple promoter architectures. In more complex systems, the potential for many inputs to transcriptional regulation will buffer the disruption of any single input.

38It must be stressed that our quantitative analysis has been carried out on one gene, in steady-state conditions. The gene is strongly expressed and actin is usually put under the umbrella of 'housekeeping', which is a slight simplification as most of the Dictyostelium actin gene family, including act5, show some developmental regulation, at least at the transcript level. A more strictly induced developmental gene might be expected to show more strict two state ON/OFF behaviour. However, even in cases with more prominent bursting we would argue that the initiation rate in the active state is likely to fluctuate over time and differ from cell to cell. This is supported by recent studies on mammalian transcriptional induction by serum and growth factors, which also suggest more complexity than a standard 2-state model provides, with evidence from luciferase reporter fluctuations and measurements of nascent RNA by smFISH indicating modulation of the transcription rate within the ON state. A more recent smFISH study quantifying nascent RNA revealed modulation of burst size and frequency by cell size and cycle stage, respectively. Whilst we have argued here for a more complex view of transcriptional regulation, it must be considered that the act5 promoter is less than 700 bp long, the gene contains no introns and around 60% of the Dictyostelium genome encodes protein, and so provides little scope for long range regulatory interactions. In the light of these features, we suggest the spectrum of states is likely to have considerably more scope for complexity in a mammalian cell. The continuum of states we infer is likely to be a more realistic view of gene activity fluctuations than the standard views of a small number of fixed discrete states. There are perhaps a hundred different proteins involved in a standard eukaryotic transcription reaction, even ignoring the components of the chromatin template. The likely configurational complexity, in addition to the potential for modulation by protein modification and nuclear context, seems consistent with the continuum view.

39For targeting of MS2 repeats into the actin 5 gene, we utilized a BsrGI-SpeI restriction fragment containing 24 MS2 repeats (1.3 kb) upstream of the blasticidin resistance (bsr) cassette. The resistance cassette is flanked by loxP sites for CRE-mediated removal of the marker, allowing transcription to terminate at the natural 3’ sequence of the gene. For the 5’ tagging of act5, the MS2-bsr was cloned between a promoter fragment of the gene (-680 to +21) and a gene fragment (+108 to +1313), using BsrGI and SpeI sites. Cells derived using this 5’ tagging vector were used for all experiments described in this paper, unless indicated otherwise. For the gene replacement vector, the MS2-bsr was cloned between the same promoter fragment and a fragment from the 3’ coding sequence and terminator of the gene (+1092 to +1899). For the 3’ targeting vector, we used this same 3’ region combined with a 5’ region derived from the act5 coding sequence (+259 to +1113). The ATG corresponds to +1. The translational STOP, TAA is at +1129. The first clear polyadenylation motif (AATAAA) starts at +1193. Sequences were checked at each cloning step to ensure plasmid stocks retained the correct sequences (as specified by dictyBase). For TATA mutations, the T1A (AATAAAT) and A2C (TCTAAAT) mutant promoter fragments were generated by gene synthesis, then spliced into the wild-type promoter sequence using the BstEII site upstream of the TATA box, prior to inserting the MS2-bsr fusion, again using BsrGI and SpeI sites at the same positions.

40Targeting fragments were released from cloning vectors by digestion with polylinker enzymes ClaI and NotI. These targeting fragments were transformed into a Dictyostelium AX3 clone previously engineered to express a red fluorescent nuclear marker, H2Bv3-Cherry, under the control of the endogenous promoter of the rps30 gene. Targeted clones were identified by PCR, then genomic DNA Southern blotted to ensure MS2 repeat integrity and single correct insertions in the targeted clones. Correct targeting of TATA mutations to the act5 promoter was checked by sequencing of PCR products from recombinant clones. Positive clones were then transiently transfected with a plasmid expressing the CRE recombinase, to remove the bsr cassette, allowing the MS2 RNA to fall under the control of the natural act5 terminator. Clones were then transformed with an extrachromosomal vector, based on pDM1096 (from Dr. Douwe Veltman) expressing the MCP-GFP fusion protein, to permit detection of nascent RNA in living cells. Selection used 20 µg/ml G418.

41Cell culture preparation for imaging was carried using cells grown in HL5 medium (FORMEDIUM) supplemented with penicillin+streptomycin. 20 µg/ml G418 selection was added 72 hr after thawing frozen stocks. 18 hr prior to imaging, cells were split into imaging chambers (NUNC LabTek-II) at the appropriate density for imaging the next day at around 20% confluency, drug selection was removed and HL5 was replaced with imaging medium (75% LoFlo medium (Formedium), 10% FBS, 15% HL5). Imaging media was refreshed 1.5 hr before imaging. Imaging was performed using an UltraView Vox spinning disc confocal microscope. Objective, laser lines, camera (Hamamatsu C9100-13 EM-CCD) settings, laser powers and exposure times were optimized to minimize photobleaching and ensure negligible photo-toxicity (measured in terms of the average transcription spot intensity - in trials we found transcription spots were attenuated before reduced cell motility or cell rounding is observed). Data were analysed using custom-built software integrating both cell tracking and spot detection. Code can be accessed at http://www.ucl.ac.uk/lmcb/sites/default/files/Corrigan2016MatlabFiles.zip.

42Cells were imaged live as described above for a single time point and then immediately fixed and prepared for smFISH measurements following the procedure outlined in, using a single probe (CATGGGTGATCCTCATGT; Biosearch, Petaluma, CA) against the repeated spacer between the each MS2 loop sequence, end-labelled with Quasar 570 fluorophore. Nuclei were co-stained with DAPI. Image stacks were acquired using the spinning disc using exposure times of 3 s for single molecule sensitivity. Individual cytoplasmic RNA molecules and nascent transcription sites were detected using FISH-quant software v2, augmented with a custom-written dual-threshold algorithm to segment nuclei and cytoplasm of cells based on the DAPI signal.

43A Monte Carlo simulation framework for the MS2 system was constructed using custom-written tools in MATLAB. For simulation of FRAP experiments, binding and potential unbinding of MCP was simulated with first order kinetics. Tools for analysis of nth-order hidden Markov models, including the forward-backward and Baum-Welch algorithms and Gibbs sampling, were custom-written using MATLAB. Full details of the simulation and probabilistic analysis procedures are presented in the appendix.