Abstract
Genes are expressed in stochastic transcriptional bursts linked to alternating active and inactive promoter states. A major challenge in transcription is understanding how promoter composition dictates bursting, particularly in multicellular organisms. We investigate two key Drosophila developmental promoter motifs, the TATA box (TATA) and the Initiator (INR). Using live imaging in Drosophila embryos and new computational methods, we demonstrate that bursting occurs on multiple timescales ranging from seconds to minutes. TATA-containing promoters and INR-containing promoters exhibit distinct dynamics, with one or two separate rate-limiting steps respectively. A TATA box is associated with long active states, high rates of polymerase initiation, and short-lived, infrequent inactive states. In contrast, the INR motif leads to two inactive states, one of which relates to promoter-proximal polymerase pausing. Surprisingly, the model suggests pausing is not obligatory, but occurs stochastically for a subset of polymerases. Overall, our results provide a rationale for promoter switching during zygotic genome activation.
Similar content being viewed by others
Introduction
In all eukaryotes, transcription of active genes into RNA requires the controlled assembly of multiple protein complexes at promoters1. This includes the sequential recruitment of general transcription factors to form the pre-initiation complex (PIC), followed by the recruitment of RNA Polymerase II (Pol II) at the transcriptional start site (TSS)2,3,4 of promoters. Among all these complexes, TFIID stands out as the primary core promoter recognition factor that triggers PIC assembly5. TFIID includes TATA-binding protein (TBP), which binds upstream of the TSS, as well as 13 TBP-associated factors (TAFs), known to bind downstream promoter elements such as the Initiator element (INR) and the downstream promoter element (DPE) motifs3,6.
After initiation, Pol II transcribes a short stretch of 30–80 nucleotides before pausing7. This step is regulated by the TFIID complex, as directly demonstrated in vitro8 and inferred from the over-representation of TFIID-bound core promoter motifs (INR, DPE, pause button) in highly paused genes9,10,11,12,13. Pausing durations are highly variable among genes14,15,16. However, it is thought that exit from the paused state is a kinetic bottleneck in the transcriptional cycle17 and is often used as a checkpoint during development to foster coordination in gene activation, plasticity or priming18,19,20,21. How core-promoter motifs affect this rate-limiting step is still unknown, particularly in the context of a developing embryo.
In parallel to these genomics approaches, single-cell imaging revealed that transcription is not a continuous process over time but occurs through stochastic fluctuations between periods of transcriptional activity and periods of inactivity, called bursting22. Direct labeling of newly synthesized RNA with the MS2/MCP amplification system23,24 has been the method of choice to observe these bursts in real time in living cells25,26,27 or multicellular organisms28,29. In the context of a developing embryo, these studies revealed how enhancers regulate bursting during pattern formation30,31,32. However, relatively less attention has been given to the impact of core promoter motifs on transcriptional bursting.
To build a mechanistic understanding of bursting beyond a qualitative description of burst size and frequency, it is critical to consider the various timescales of the transcription process and to employ mathematical modeling explicitly describing various promoter states22. Indeed, depending on the gene and the cellular context, transcriptional bursting has been shown to occur at multiple timescales33,34, from seconds (e.g., polymerase clusters, Pol II firing35,36,37), to minutes (e.g., TBP binding/unbinding, transcription factor binding38,39) to hours (e.g., nucleosome remodeling, chromatin marks37,40). These different timescales can all occur at a single gene, and the term multiscale bursting was coined to describe such complex bursting kinetics34,37. This multiscale bursting might explain why the simple two-state model (random telegraph), whereby promoters switch from an active to an inactive state, does not always suffice to reliably capture promoter dynamics36,37,41,42.
While it is well understood that cis-regulatory sequences must impact transcriptional bursting, with enhancers primarily affecting burst frequencies43,44, a detailed dissection of the impact of promoter motifs on promoter-state dynamics is lacking. In this study, we sought to examine the role of core promoter sequences on transcriptional bursting in Drosophila embryos. In particular, we focus on two core promoter motifs, the TATA box (TATAWAWR) and the INR (TCAGTY in Drosophila6) which represent pivotal core promoter contact points with the TFIID complex45,46, and are known to regulate both initiation and promoter pausing8,12.
The Drosophila early embryo is an ideal system to decipher transcriptional bursting regulation since spatially distinct patterns of gene expression are deployed within a relatively short time frame in a multinucleated syncytium, rapidly dividing, that is highly amenable to quantitative imaging47. At this early stage of development, the zygotic genome, initially transcriptionally silent, awakens progressively while cell cycle durations lengthen, a process known as the maternal-to-zygotic genome transition (MZT)48. This zygotic genome activation (ZGA) occurs gradually through two waves, a major wave during nuclear cycle 14 (nc14) and a minor wave occurring at earlier stages. Moreover, Drosophila developmental genes show a clear promoter code with well-defined promoter elements6,49, with differential usage during ZGA50.
We employed the MS2/MCP system to monitor nascent transcription. We implemented a machine-learning method that deconvolves single-nuclei mRNA production from live Drosophila embryos to detect all single polymerase initiation events. The waiting times between successive initiation events were further analyzed to infer the number of promoter states and the transition rates among these states. Our results show that TATA box-containing promoters are highly permissive to transcription through long ON durations and high Pol II firing rates. However, the presence of an INR in the core promoter necessitates the use of a three-state model, with a second inactive promoter state that is likely a consequence of Pol II pausing. We propose a renewed view of promoter pausing whereby only subsets of polymerases enter into a paused state.
Results
A synthetic platform to image promoter dynamics in living embryos
To examine how the core promoter sequence influences gene expression variation, we developed a synthetic platform whereby several core promoters were isolated, cloned into a minigene, and inserted at the same site in the Drosophila genome (Fig. 1a). This approach allows for a direct comparison of core promoter activity since differences caused by variations in cis-regulatory context, genomic positioning, and mRNA sequence (particularly the 3’UTR) are eliminated. Core promoters were inserted immediately downstream of the snail (sna) distal minimal enhancer (snaE)51. The core promoters were selected based on the presence of known core promoter motifs49, such as a TATA box in sna or an INR in kruppel (kr) and Insulin-like peptide 4 (Ilp4) (Fig. 1a, Supplementary Fig. 1, and Supplementary Movies 1 and 5). We also selected the brinker (brk) core promoter, which is devoid of any known canonical promoter motifs. None of these four promoters have a DPE, however, Ilp4 possesses a bi-partite bridge element49. These developmental genes are endogenously expressed during nc1452,53 and precise TSS positions were established using embryonic CAGE (Cap Analysis of Gene Expression) datasets54 (Supplementary Fig. 1). For all transgenes, we used 100 bp of a promoter sequence, as previous work established that minimal promoter sequences are sufficient to establish pausing in vivo12,21.
To track transcription, 24 MS2 stem loops were placed in the 5’ UTR downstream of the promoter, followed by the sequence of the Drosophila yellow gene, a gene fragment used as a reporter because of its large size and lack of endogenous expression in the early embryo55. MS2 stem loops in transcribed mRNA were bound by a maternally supplied MCP-GFP fusion protein, which allowed the detection of transcriptional foci as bright GFP spots. We imaged living Drosophila embryos at a high temporal resolution of one frame per 3.86 s and focused on the first 30 min of nc1456. To ensure that all quantified nuclei experience similar peak levels of transcriptional activators, we restricted our analysis to a defined domain within the mesoderm (25 µm on either side of the presumptive ventral furrow; Fig. 1b)57. Signal intensity at the transcription site was retrieved for each nucleus after 3D detection and tracked throughout nc14 (with mitosis between nc13 and nc14 considered as time zero; see “Methods” section) and then associated with their nearest nucleus. This generated individual 4D nuclear trajectories (Supplementary Fig. 2; “Methods”).
Core promoters differentially affect transcriptional synchrony
We hypothesized that differences in core promoter motifs would result in variability in gene expression independently of specific enhancer contacts. To test this, we characterized single-nucleus transcriptional activity (Fig. 1c–e), as well as the initial timing of activation. As previously reported with fixed sample approaches21,58, we observed a spectrum of synchrony profiles (Fig. 1f), defined as the temporal coordination of activation among a spatially defined domain. The TATA-driven sna and INR-driven kr minimal promoters showed a rapid activation with a t50 (the time needed for 50% of nuclei within the region of interest to show GFP puncta) of 11 min and 14 min, respectively while the INR-driven Ilp4 promoter was delayed with t50 of 21 min. The brk core promoter did not surpass 24% of nuclei activated throughout nc14.
In addition to differences in synchrony, the overall mRNA production also varied between transgenes (Fig. 1g). Within our four developmental promoters, promoters showing more rapid reactivation in nc14 also had higher TS intensity in active nuclei, indicating increased instantaneous mRNA production. To examine the total mRNA output, we looked to the integral amplitude, or the area of the curve defined by the average MCP-GFP signal plotted over time (Fig. 1d). The integral amplitude showed that promoters mediating faster activation and higher instantaneous mRNA production had a higher total mRNA output (Fig. 1h).
Thus, we established an experimental setup and quantification pipeline that allows disentangling the contribution of minimal promoter sequences on transcriptional activation of individual, naturally synchronized nuclei within a well-defined, homogeneous spatial domain in live embryos. This initial quantitative comparison of four natural developmental promoters suggests that those with a canonical promoter motif (either TATA or INR) tend to produce higher levels of expression.
A machine-learning method to infer promoter-state transition rates
With current labeling and imaging technology, individual transcriptional initiation events cannot be easily and directly visualized, tracked, and quantified in the early Drosophila embryo. Live imaging of transcription typically shows spots comprised of several newly synthesized mRNAs resulting from the action of multiple polymerases. In order to calibrate fluorescent signals from live imaging, we used single-molecule hybridization experiments28. Using the fluorescence of a single mRNA molecule, we estimated the average number of mRNA molecules present at the TS within the nucleus at a steady state (Supplementary Fig. 3). This calibration step allowed us to express single-nuclei TS intensities as an absolute number of transcribing polymerases (Fig. 2a).
We then used a recently described novel machine-learning method to infer the transcriptional bursting mechanism59. This method involves three major steps: detection of successive initiation events for each nucleus, multiexponential parametric regression of the distribution of waiting times between successive events, and identification of Markovian promoter-state transition models.
In order to detect initiation events, we considered that each trace results from the convolution between (i) the sequence of initiation events marked by a rise in GFP intensity plotted over time, and (ii) the signal produced by a single polymerase (Fig. 2b)33,59. The deconvolution procedure uses a genetic algorithm to determine optimal Pol II positioning within the gene body (Fig. 2c) and thus the time between successive Pol II initiation events (Δt) (Fig. 2d). In this first step, some specific parameters are fixed, such as the temporal resolution (3.86 s), the speed of Pol II elongation (45 bp s−1 in Drosophila embryos60), the size of the transcript (5.6 kb), and the retention time of the mRNA at the transcription site (assumed to be small relative to the time required to produce a transcript).
Having obtained temporal maps of individual Pol II transcription events, the method then statistically analyzes the distribution of waiting times between two successive polymerases (Δt)59. A multiexponential regression fitting was applied to the distribution of Δt, indicating the number of rate-limiting transitions required to best fit the data (“Methods”). We used confidence intervals and a Kolmogorov–Smirnov test to rigorously determine the smallest number of rate-limiting steps fitting our experimental data (Supplementary Table 1). Following the principle of parsimony and to avoid overfitting, models with a larger number of steps and more parameters were not retained even if they also fit well. Importantly, this step uncovered the number of characteristic timescales of transcriptional fluctuations that are tantamount to the number of promoter states.
Finally, a transcriptional model of the promoter was established, with multiple states and timescales inferred from the parameters of the multiexponential Δt distribution (Supplementary Table 2). This step permitted the estimation of the kinetic parameters such as the kON (rate of switching from a transcriptionally nonpermissive state to a permissive one), kOFF (switching rate from a transcriptionally permissive to a nonpermissive state), and kINI (the rate of Pol II initiation events once the promoter is in a transcriptionally permissive state) for the simplest two-state model. However, kinetic estimates of more complex models (three or more promoter states) can also be derived from the Δt distribution. The accuracy and robustness of the deconvolution method were tested by applying it to artificial data with known positions of transcription initiation events and known kinetic parameters (Fig. 2e–h). These artificial data were generated using the Gillespie algorithm and proved the robustness of the approach (see “Methods”).
Importantly, this approach provides a time map of transcription events in a cell population in a model-independent manner. The number of promoter states is evaluated during the multiexponential fit procedure, but with no a priori constraints on the result59. This contrasts with current methods which directly fit a particular transcription model to the data, such as methods based on the autocorrelation functions61,62,63, maximal likelihood estimates36, or Bayesian inference30,64,65.
The sna promoter as a model to investigate the impact of core promoter on transcriptional dynamics
To investigate the quantitative dynamics of the sna promoter, we applied this method to hundreds of nuclei pooled from snaE < sna embryos. A heatmap of positions of single initiation events for each nucleus revealed the sna promoter drove frequent transcriptional initiation events. Silent periods with no initiating polymerase were very rare (Fig. 2i–j, white bars), consistent with the results obtained from fixed embryos21,66. We found that a biexponential (sum of two-exponential functions) fitting correctly fit the survival function of polymerase waiting times (Fig. 2k). The biexponential fitting is only compatible with a two-state model (Fig. 2k). Thus, the transcriptional activity of the sna promoter can be described by the random telegraph model, with a simple random switch between an inactive OFF and a permissive ON state, from which transcription initiates at a given rate (kINI; Fig. 2l). The probability to occupy the permissive ON state was high in the case of sna (0.91; Fig. 2m). The TON was estimated to be 242 s, with a TOFF of 24 s. The kINI of sna promoter was estimated at one initiation event every 9 s (Fig. 2m), consistent with the inferred initiation rate at its endogenous locus66 and in the range of recent estimates of initiation frequencies for other developmental genes in Drosophila32 and in mammalian cells37. Collectively this shows that we can locate individual initiation events, determine the number of promoter states, select a promoter model by the principle of parsimony, and then estimate average promoter switching rates in vivo.
The TATA box regulates the ON and OFF duration but not the number of states
To determine how the TATA box influences transcriptional kinetics in vivo, we developed a series of sna core promoter mutants (Fig. 3a) where the TATA box was replaced with either the TATA-like sequence of kr (snaTATAlight) that is bound by TBP50 (Supplementary Fig. 1) or a non-TATA sequence where the first four bases are mutated (snaTATAmut). Surprisingly, the strong TATA mutation did not completely abolish the activity of the sna promoter, although this mutant promoter was devoid of any known core promoter motifs. In comparison to the sna promoter, the synchrony of the sna TATA mutants was reduced in a graded manner relative to the fidelity of the TATA box (Fig. 3b, Supplementary Movies 1–3). The sna promoter reached synchrony at ~11 min into nc14, while the snaTATAlight promoter reached synchrony at 24 min and the snaTATAmut promoter never reached synchrony. This implies that the TATA box influences synchrony of activation, most likely by promoting and stabilizing TBP binding (Supplementary Fig. 1).
Similar to the activation profiles, instantaneous TS intensities were directly correlated with TATA box sequence fidelity (Fig. 3c). Interestingly, while both the sna promoter and snaTATAlight promoters appeared to reach a steady state, the snaTATAmut was unable to do so. This may have been the result of MCP-GFP signal being too low to consistently surpass the detection threshold, as single-molecule FISH experiments indicated most nuclei in the ventral region were at least weakly active in nc14 snaTATAmut embryos (Supplementary Fig. 4a–c). The total mRNA production of each promoter measured by the integral amplitude was similarly affected (Fig. 3d), with the sna promoter driving the highest total mRNA expression while the TATA box mutants expressed reduced amounts of mRNA correlated to the fidelity of the TATA box.
We next asked whether the snaTATAlight and snaTATAmut transgenes still followed a two-state model. For each nucleus of each of these three transgenic genotypes, we located single Pol II initiation events during the first 30 min of nc14 (Fig. 4a–f and Supplementary Fig. 4d–f). By examining the survival curve of Pol II waiting times, the dynamics of the sna TATA mutants appeared to be accurately recapitulated by the two-state random telegraph model, similarly to the unmutated sna promoter (Supplementary Fig. 4g–i, Supplementary Table 1). We then analyzed which kinetic parameters changed in the mutant promoters. Extracting kinetic parameters for the snaTATAlight transgene revealed a relatively mild effect on initiation rates (1 Pol II per 13 s for snaTATAlight instead of 1 Pol II per 9 s) (Fig. 4g–i and Supplementary Table 2). In contrast, we observed a strong increase in the TOFF (Fig. 4g–k) for the snaTATAlight promoter relative to the sna promoter (Fig. 4g–i, k and Supplementary Fig. 5b). The strong TATA mutation led to a similar trend, although due to the weak activity of this promoter, the number of analyzed nuclei for this mutant is lower than for the other promoters. The kinetic parameters of TATA mutant promoters led to an overall reduction in burst size, defined as the number of transcripts produced during an active period (Fig. 4l). The reduction was primarily due to a decrease in the duration of the ON periods, consistent with a destabilized TBP/TATA box interaction or reduced TBP recruitment to the promoter.
The impact of TATA on the duration of ON/OFF promoter states observed in Drosophila embryos is also in agreement with results from live imaging of HIV-1 transcription in human cells37, and with genomic studies in cultured mammalian cells, where TATA-driven promoters are associated with large burst sizes67. They differ from those obtained with mutations of the native actin promoter in Dictyostelium36, which may be due to differential behavior between induced genes and constitutive genes or to compensation by other genomic elements36.
We conclude that in Drosophila embryos, the TATA box largely controls gene expression through lengthening the ON state duration at the expense of the OFF state duration, but does not alter the number of promoter states.
The INR motif induces a third rate-limiting promoter state
At the biochemical level, transcriptional activation is a multistep process requiring the orchestration of many factors3, and yet TATA-driven promoters can be modeled with a simple two-state model corresponding to one kinetic rate-limiting step in this complex process. We hypothesized that additional rate-limiting steps were likely to exist and could be identified using our imaging-based methodology. As the TATA box did not affect the number of promoter states for the sna promoter, we examined the effect of another core promoter motif, the INR. In Drosophila cells and embryos, the INR is associated with stably paused genes, genome wide7,9,15,68. Moreover, analysis in cell culture indicated that genes with increased pausing stability tend to harbor an INR motif in the core promoter, and in particular a preferential G at the +2 position12. We, therefore, reasoned that manipulating the INR motif in embryos might affect pausing and therefore may induce changes in rate-limiting promoter states.
To examine the role of the INR core promoter motif, we created a series of transgenes in which we manipulated the INR without changing the enhancer or the downstream gene sequence (Fig. 5a). The sna core promoter does not have an INR, while the core promoter of kr has a natural INR sequence with a G at +2 and is stably paused16,18. Moreover, the kr promoter possesses a non-canonical TATA box (TATAlight, Fig. 1a). We exchanged the TSS region of sna with the INR of the kr promoter (sna+INR). We also created a transgene (Fig. 5a) whereby the INR of kr was replaced with the TSS region of sna (kr-INR1). Interestingly, neither loss of the INR in the kr-INR1 transgene nor gain in the sna+INR transgene strongly affected synchrony (Fig. 5b and Supplementary Movies 1, 4, 5). The addition of an INR to sna did not dramatically affect mRNA output (Fig. 5d). Both the sna+INR and kr-INR1 transgenes had higher instantaneous activity than the cognate wild-type promoters (Fig. 5c). Moreover, a second mutation of the INR in kr (kr-INR2) had similar synchrony and mRNA production (average instantaneous intensity) as in the wild-type kr promoter (Supplementary Fig. 6a–c). We then applied the deconvolution procedure to each genotype (Fig. 6a and Supplementary Figs. 6d–k and 7a, d–k). Surprisingly, our analysis of polymerase initiation events indicated that, for sna+INR nuclei, a two-state model was not sufficient to fit the data (Fig. 6a, b, Supplementary Fig. 7i, and Supplementary Table 1). Instead, the sna+INR survival function was well fitted by adding an extra exponential term. Thus, a three-state model appropriately recapitulated sna+INR promoter dynamics (Fig. 6a, b and Supplementary Fig. 7i–iʼ). Similarly, the kr promoter also required a three-exponential fitting (Fig. 6b and Supplementary Fig. 6h–hʼ). In contrast, removing the INR from this promoter (kr-INR1 promoter) led to a two-exponential fitting for the survival function (Fig. 6b and Supplementary Fig. 6i–iʼ). Similar results were obtained with a second INR mutant (kr-INR2) (Supplementary Fig. 6j–jʼ). Furthermore, the natural INR-driven Ilp4 promoter was also associated with a three-exponential fitting, and mutation of its INR sequence (Ilp4-INR) resulted in a reversion to a two-exponential fit (Fig. 6b and Supplementary Fig. 8d, e, h, iʼ). These results, therefore, suggest that the INR motif is associated with a third rate-limiting promoter state.
Pol II pausing is associated with a third promoter state
Given the correlation between the presence of an INR motif and the number of characteristic timescales and thus the number of promoter states (three in the presence of an INR), we reasoned that one of these states could be polymerase pausing. To test this hypothesis, we impaired the establishment of pausing by reducing the expression of the largest subunit of the pause-inducing negative elongation factor complex (NELF), NELF-A. Previous work demonstrated that depletion of NELF by RNAi globally reduced promoter-proximal pausing in Drosophila cells13, and that loss of the human C-terminal NELF-A “tentacle” blocked stabilization of pausing69. We combined the inducible RNAi/GAL4 system with our MS2/MCP-GFP reporter to examine the transcriptional dynamics of the three-state kr promoter. We performed live imaging of the kr-MS2 transgene and compared a control white RNAi knockdown with embryos harboring a maternal depletion of Nelf-A (evaluated to >85% knockdown by RT-qPCR, Supplementary Fig. 9a, Supplementary Movies 7 and 8). Remarkably, reducing Nelf-A transcript levels led to a change in kr promoter dynamics, with a reversion to a two-promoter-state dynamic and much less frequent long waiting times (Fig. 6b and Supplementary Fig. 9e-fʼ, arrowheads, Supplementary Table 1). These data, combined with the effect of INR mutations in cis, led us to propose that one of the two inactive promoter states observed with the three-state promoters could be associated with promoter-proximal polymerase pausing.
To confirm that the INR motif was truly altering Pol II pausing, we performed Pol II ChIP-qPCR on our promoter transgenes from staged embryos. We evaluated the strength of pausing by computing the pausing index. The pausing index (PI) was obtained by quantifying the Pol II binding at the promoter relative to the gene body70 (Fig. 6d). As shown previously, in our transgenic embryonic extracts encompassing various stages from nc13 to late nc14, the wild-type sna promoter is paused18,21,50, but the addition of an INR motif significantly increased its PI (Fig. 6e), consistent with INR-mediated lengthening of the pause duration observed in cell culture12. To examine the inverse scenario, we also quantified the PI of the kr and kr-INR1 transgenic promoters. Similar to what is observed in the endogenous promoter18, the kr transgene was highly paused (Fig. 6e). Mutation of the INR in the kr-INR1 transgene reduced pausing (Fig. 6e). Likewise, the Ilp4 promoter from our transgene showed a high level of pausing, which decreased upon mutation of its INR sequence (Fig. 6d). Both the sna+INR and kr promoters also show increased NELF-E enrichment at these promoters by ChIP-qPCR on embryonic extracts (Supplementary Fig. 9g). Taken together, characterization of the INR core promoter motif suggests the presence of a strong INR motif results in stabilization of pausing in vivo.
Pol II does not undergo systematic pausing
Next, we envisaged how pausing could be modeled by three distinct promoter states. The current view of pausing is that all polymerases systematically enter into pausing, a step followed by either a productive elongation or a termination17,71,72,73. We therefore initially tested if our data could be explained by a model consisting of three promoter states (inactive OFF, permissive ON, and paused) where all productive polymerases undergo pausing prior to transitioning to the permissive ON state (Fig. 6c). This model was clearly incompatible with our data, as the fitting exceeded the bounds of the 95% confidence interval (Supplementary Figs. 6l and 7i”). Recent work on the prototypic model of a highly paused promoter, the HIV-1 promoter, also found that modeling transcriptional dynamics with obligatory pausing was not in agreement with live-imaging data in HeLa cells, and instead proposed an alternative model of pausing where only a fraction of polymerases are subject to pausing in a stochastic manner59. We asked if this alternative non-obligatory pausing model could be applicable to our findings (Fig. 6c). According to the goodness of fit, nonsystematic pausing was compatible (Supplementary Figs. 6hʼ and 7iʼ). Our machine-learning method enabled the estimation of the kinetic transition rates between these states, summarized in Supplementary Table 2 and Supplementary Fig. 5a–c.
We conclude the presence of an INR motif translates into a longer paused state that creates an additional rate-limiting step during transcription in vivo. In our kinetic model, the inactive states could in principle equally correspond to pausing. These two-promoter states are discernible by their distinct timescales, on the order of seconds and minutes, respectively (Supplementary Table 2). Because the NELF knockdown dramatically decreases the frequency of long waiting times (Supplementary Fig. 9e–f), we favor the longer-lived inactive state, which lasts on the order of minutes, as the potential pausing state. The second short-lived inactive state, lasting on the order of seconds and present for both TATA- and INR-containing promoters, would then correspond to a nonpermissive promoter state, possibly TFIID-unbound as its frequency increases in TATA mutants (Fig. 4k and Supplementary Table 2). As the sna and kr endogenous promoters drive a high level of gene expression (Supplementary Fig. 1), it is reasonable to hypothesize that such a nonpermissive state is transient and infrequent (Supplementary Table 2).
INR control of promoter switching kinetics
Like the natural sna promoter, the sna+INR promoter showed a high probability of occupying the transcriptionally permissive ON state, with a long lifetime of 170 s. This ON duration was slightly reduced compared to the sna transgene (Fig. 6f, i), possibly implying competition between the natural canonical TATA box and the added INR in the sna+INR combination (see below). The major difference between these two promoters was the existence of two distinguishable inactive states linked to the presence of the INR sequence. One of these was short-lived (17 s) and one long-lived (302 s) (Fig. 6f, i, Supplementary Fig. 5b, c, and Supplementary Table 2). The stable long-lived paused state was not observed with the sna promoter but did occur in the presence of the INR-containing kr promoter.
In the case of kr and Ilp4 developmental promoters, the paused state lasted ~151 and 105 s, respectively, but was reached relatively rarely (Fig. 6g, j, Supplementary Fig. 5c, and Supplementary Table 2). Pausing was not observed upon two types of mutations of the kr INR motif (Fig. 6j, Supplementary Fig. 6i–jʼ, m, and Supplementary Table 2). Thus our data favor a model where the paused state lasts several minutes (in the case of kr, sna + INR and Ilp4), but occurs infrequently. Interestingly, the sna promoter is also moderately paused in our ChIP assay and at its endogenous locus74, but this promoter is not regulated by a three promoter-state dynamic. We hypothesize that pausing is highly unstable for this promoter at this stage and does not constitute a rate-limiting step that we can capture with our live-imaging assay (see “Discussion”).
In addition to two separable inactive states, the kr transgene exhibited a highly probable permissive state (0.83) with a duration of 94 s (Fig. 6g, Supplementary Fig. 5a, and Supplementary Table 2). When the INR of kr was mutated (kr-INR1), there was a significant increase in ON duration (Fig. 6j and Supplementary Table 2). A similar effect was obtained with an alternative mutation of the INR (kr-INR2) (Supplementary Figs. 5a and 6m and Supplementary Table 2).
Interestingly, in all the transgenes, the differences in the number of states and the active/inactive durations were not accompanied by a change in Pol II firing rates. The sna, sna+INR, kr, and kr-INR1 mutant all exhibited an estimated kINI of 1 Pol II every 8–9 s (Fig. 6f–g, i–j and Supplementary Table 2). These estimates agree with recently documented polymerase initiation rates of Drosophila developmental promoters30,32,75. Our results suggest the rate of initiation is not associated with the presence or absence of an INR motif. Instead, the INR motif leads to a regulation of transcriptional bursting via two inactive promoter states, one of which is associated with stabilized pausing.
Taken collectively, our results demonstrate that transcription dynamics are differentially regulated by the INR and the TATA box. Developmental promoters containing a strong INR motif can travel through multiple inactive states due to an extra regulatory step at early elongation via Pol II pausing. Our live-imaging data favor a model whereby pausing occurs on the order of minutes but is not an obligatory state reached by all engaged polymerases. Instead, we propose that when pausing constitutes a rate-limiting step, it occurs for only a subset of polymerases with a low frequency.
Promoter sequences and timing of initiation
A second property that we found to be common to all promoters examined, was the regime of waiting times between mitosis and first transcriptional activation. The analysis of the lag time between mitosis and onset of activation is “non-stationary” and was previously modeled using a different modeling paradigm, based on mixed gamma distributions of the random time to transcription activation after mitosis in nc1457. This analysis estimates two parameters, the number of promoter states (a) and transition times between them (b). We applied this modeling framework to our data by only focusing on the distribution of waiting times between mitosis and first activation (illustrated in Supplementary Fig. 10a). Interestingly, a multiexponential distribution becomes a mixed gamma distribution when the time parameters of the exponentials are even. Conversely, the equality of these time parameters justifies the use of a mixed gamma regression. Using the more general multiexponential regression, we proved that regardless of the promoter sequence, the transition times between states during this lag regime are homogeneously distributed (Supplementary Fig. 10b, c). This is in contrast with the transition times between states in the stationary bursting regime, which are heterogeneous. This suggests that the states and transitions involved in the two regimes (nonstationary and stationary) are distinct and have separate regulations. This further supports a previous study that demonstrated the delay in postmitotic transcription activation is dependent on the enhancer sequence57.
Interplay between TATA and INR
Recent single-cell RNA-seq quantification suggests that promoters exhibiting both TATA and INR motifs produce higher burst sizes than when only one motif is present67. However, systematic interrogation of human promoters with a synthetic biology approach shows that TATA and INR additively but not synergistically increase gene expression76. A synergistic and/or additive effect does not seem to be evident in Drosophila, as TATA-containing promoters have relatively infrequent INR motifs or exhibit an INR sequence devoid of a G at +2, shown to be associated with long pause durations12,50.
To quantify TATA box and INR interplay on rate-limiting promoter switching states in vivo, we compared the sna + INR and snaTATAlight transgenes to that of a combined mutant snaTATAlight+INR promoter (Supplementary Fig. 7a). Remarkably, after applying the deconvolution procedure to snaTATAlight+INR nuclei (Supplementary Fig. 7g), the distribution of waiting times between initiation events required a fit with two exponentials similarly to snaTATAlight (Supplementary Fig. 7k–kʼ). We note however that the statistical tests allocating the number of promoter states for this particular genotype were not robust (Supplementary Table 1). In terms of promoter kinetics, the snaTATAlight+INR behaved similarly to the snaTATAlight promoter (Supplementary Fig. 5a, b). Interestingly, the snaTATAlight+INR has a p(OFF) nearly twice that of the snaTATAlight promoter (Supplementary Table 2), indicating that the INR may act to maintain a nonpermissive promoter state even in the absence of polymerase pausing.
To gain more insight into TATA/INR interplay, we also generated a kr mutant promoter that does not contain a non-canonical TATA (kr-TATA; Supplementary Fig. 6a). Data from this promoter required a three-state model (Supplementary Fig. 6k–kʼ), similar to the kr promoter. In the kr-TATA promoter, the duration and probability of the paused state increased while the ON state appeared shorter and less probable (Supplementary Fig. 5a, c and Supplementary Table 2). This indicates that the loss of TBP binding may increase pausing induced by the INR motif, potentially providing a mechanistic rationale for the relative scarcity of dual TATA/canonical INR promoters12,50 in the Drosophila genome. Taken together, these results indicate that TATA and INR are not simply synergistic or additive but display more complex interactions.
Discussion
The spatiotemporal organization of gene expression is critical to the development of a functional living organism. While we have accumulated knowledge on how enhancers precisely regulate gene expression, we know relatively little about the impact of core promoters on transcriptional bursting in a developing embryo. Here, we investigated how minimal 100 bp sequences of developmental promoters control transcriptional states and their switching kinetics. We found that a classical two-state model does not suffice to capture promoter dynamics of stably paused promoters containing a strong INR motif but does fit with TATA-containing promoters.
The development of a novel numerical deconvolution method revealed startling insights into the variation of transcriptional initiation between minimal promoter sequences in vivo. Our experimental setup allowed us to unmask gene expression variations, which are only dictated by variations in core promoter sequences. Indeed, cell cycle duration and synchrony, the concentration of input transcription factors, and chromatin states were intentionally kept constant. Our assay revealed that a similar transcriptional activity profile could be obtained from two very distinct promoter sequences via distinct modes of transcriptional initiation. While sna promoter dynamics could be explained by a simple two-state model, the kr promoter required a three-state model, with two distinct inactive promoter states, a short-lived and a longer-lived one. These distinct transcriptional regimes are unlikely due to differences in enhancer–promoter specificity, as evidenced by the synchronous and high transcriptional activation reached with these two promoters. However, knowing the widespread preferential promoter code for specific enhancers77,78, it will be interesting to use our pipeline to examine which rate-limiting step of promoter dynamics is tuned by enhancer–promoter choice.
How could such a small number of states be compatible with the numerous promoter-binding events occurring during transcription initiation? Biochemical studies reveal structural states, while imaging-based approaches unmask key rate-limiting kinetic states. Whole-genome methylation footprinting disentangled five transcription initiation states in Drosophila cultured cells15. Remarkably, the authors demonstrated that TATA-containing promoters were frequently found in PIC-bound configuration (PIC alone or PIC + Pol II). We, therefore, propose that the permissive state in Drosophila embryos corresponds to a state where promoter DNA is bound by the PIC. The inactive state, which we found to be very transient, could then possibly represent a TBP-unbound configuration consistent with TBP dynamics observed in human cells37,39,79.
By using the sna promoter as a model, we found that TATA box directly impacts promoter occupancy of the active state and allows large transcriptional bursts by promoting long ON durations. How could the presence of a TATA box permit these long ON durations? Slow TBP protein turnover at human TATA box-containing genes may foster stable long ON durations37,39,79. However, this hypothesis needs to be confirmed by quantification of TBP kinetics during Drosophila ZGA.
In its endogenous context, the snail promoter is among the first genes to be transcribed during ZGA, in a particularly constrained environment with extremely short cell cycles (<15 min). Remarkably, the majority of genes expressed during this critical period are highly expressed, short, and intron-less with a canonical TATA box and are generally nonpaused18,50. Subsets of these, including snail, are considered “dual promoters”, as they gradually acquire pausing as development proceeds50. Thus, it is possible that the kinetic bottlenecks regulating transcription of developmental promoters evolve as developmental timing proceeds, with paused polymerase gradually emerging as a rate-limiting step. In our study, the promoter model sna shows a moderate level of pausing (as measured by ChIP) but is regulated by a simple two-state model in nc14. This could indicate that pausing occurs on a timescale indistinguishable from the OFF state and is thus embedded within it. Strikingly, early transcription in zebrafish embryos occurs in similarly constrained rapid cell cycles where a subset of zygotically activated genes also display a T/A-rich WW-box motif80. Thus, regulation of transcription via a unique rate-limiting step (OFF to ON transition) that is TBP-dependent might be conserved between fly and vertebrate embryos.
In this work, we provided quantitative evidence that INR-containing promoters are associated with a three-state model (ON, paused, OFF), contrasting with TATA promoters.
Interestingly, genomics studies revealed that transcription initiation code evolves during ZGA in flies and vertebrates50,54,80. Given the results of this study and those obtained in Drosophila cultured cells12,16, it is tempting to link the gradual emergence and stabilization of pausing during ZGA with the presence of an INR, in particular an INR with a G at +2 position12. In light of our results, we propose that the switch in promoter usage from TATA-driven to INR-driven during ZGA may lead to a change in transcriptional dynamics from two to three states, to include an elongation-mediated checkpoint. Acquisition of an extra rate-limiting step might help control cell-to-cell expression variability as well as fine-tune gene expression levels.
Taken collectively, this study establishes our promoter imaging assay and novel mathematical deconvolution and modeling methodology as a valuable tool to probe gene expression dynamics during development. Quantitative analysis of promoter dynamics at high temporal resolution opens the door to a deeper insight into the molecular mechanisms underlying transcriptional regulation in vivo. Future studies involving direct manipulation of pause initiation or duration in living embryos using approaches such as optogenetics81 combined with this framework will help to establish a broader understanding of the nature of promoter states and the role of Pol II pausing in vivo.
Methods
Drosophila stocks and genetics
All crosses were maintained at 25 °C. Transgenic lines were maintained as homozygous stocks. For live imaging, homozygous males carrying the transgene of interest were crossed with homozygous females bearing the MCP-eGFP-His2Av-RFP construct. For smiFISH, homozygous flies were crossed to yw in order to facilitate single-molecule detection. nos:GAL4-VP16 was recombined with MCP-eGFP-His2Av-RFP57. UAS:white RNAi (#35573) and UAS:Nelf-A RNAi (#32897) were obtained from the Bloomington Stock Drosophila Centre (University of Indiana, Bloomington).
Cloning and transgenesis
The sna distal enhancer-24xMS2-y minigene has been previously described51,57. Promoters (Supplementary Fig. 1 and Supplementary Table 3) were amplified from genomic DNA using Q5 polymerase (New England Biolabs) and inserted between the enhancer and 24xMS2 sequences using restriction enzyme-mediated ligation. Mutations were performed with the QuikChange II Site-Directed Mutagenesis kit (Agilent Technologies) or synthesized (Twist Biosciences) and inserted using restriction-mediated cloning. All constructs were sequenced to ensure appropriate insertion. Transgenic flies were generated using PhiC31-mediated recombination (Best Gene, Inc.), and all constructs were inserted into the same genetic background and genomic position (BL 9750). Stocks are available upon request.
Live imaging
Embryos were permitted to lay for 2 h prior to collection for live imaging. Embryos were hand dechorionated and mounted on a hydrophobic membrane prior to immersion in oil to prevent desiccation, followed by the addition of a coverslip.
Live imaging was performed with an LSM 880 with Airyscan module (Zeiss). Z-stacks comprised of 30 planes with a spacing of 0.5 µm were acquired at a time resolution of 3.86 s per stack in fast Airyscan mode with laser power measured using a ThorLabs PM100 optical power meter (ThorLabs Inc.), and maintained across embryos at 3.8 μW for constitutive MCP-GFP/His2A-RFP expression, and 5.0 μW for RNAi analyses. All movies were performed with the following settings: GFP excitation by a 488-nm laser and RFP excitation by a 561 nm were captured on a GaAsP-PMT array with an Airyscan detector using a ×40 Plan Apo oil lens (NA = 1.3) and a 3.0 zoom on the ventral region of the embryo centered on the presumptive ventral midline. Resolution was 512 × 512 pixels with bidirectional scanning. Airyscan processing was performed using 3D Zen Black v3.2 (Zeiss).
Live-imaging analysis
Visualization and analysis of the time series were performed using a custom-made software developed in PythonTM that permits visualization of each analysis step and manual correction if necessary (Supplementary Fig. 2). Activation time traces were collected starting from the Airyscan-processed Z-series described above. Green (MS2) and red (His2Av) channels were clipped to consider only after the start of nc14 as defined by the progression of anaphase across the region of interest. Nuclei were maximum intensity projected and pre-smoothed with a Gaussian filter and then thresholded with an Otsu threshold value. The resulting connected components of the binary images were then labeled and touching nuclei segmented with a watershed algorithm. The software enabled manual correction of the segmentation. Nuclei were finally tracked across the time frames using a minimum distance criterion plus a user-defined distance threshold. Nuclei appearing for only a few frames or those touching the border were removed from the analysis to exclude sources of errors.
GFP puncta representing transcription sites were analyzed in 3D. Because the duplication of DNA occurs relatively early in nc14 and because of the immediate proximity of sister chromatids30,82, it is challenging to independently resolve individual sister chromatid signals using live imaging. This is why the GFP signal at each transcriptional spot can be considered as the sum of both sister chromatids that we treat as a single transcriptional trace. For each time frame, the 3D image was filtered with a 3D Laplacian of Gaussian filter and then thresholded. The threshold value (THR) was expressed as µ + THR * σ, where µ and σ are the average and the standard deviation of the pixels values of the filtered image respectively, while THR is a user-defined value. The threshold value was in this way rescaled with respect to statistical properties of the filtered image, making THR a value independent of the particular data acquisition. All detected spots were filtered to remove (1) all the spots with a volume less than a user-defined volume threshold, and (2) spots present in only one z-frame. For each time frame, detected spots were associated in 2D to the overlapping or closest nucleus, inheriting the tracking from them; a user-defined distance threshold between nucleus and spot was used in order to avoid mis-associations.
Finally, each nuclei-puncta pair passing filtering was described as a time series of intensity, volume, and position. To eliminate intensity variation within the Z-series, spot intensity values were divided by the background fluorescence of the average intensity value of the pixels surrounding the independent spots. Analysis was restricted to the region 25 µm on either side of the center of the gastrulation furrow present at the end of nc14 as positioned using a maximum intensity tile-scan of the entire embryo to determine the coordinate position of the Z-stack. From this data, it is possible to extract: the timing of activation measured as the first timepoint GFP fluorescence crosses the software detection threshold for >1 Z-stack; cumulative activation; intensity profiles for individual nuclei; and individual nuclei burst trajectories. Integral amplitude was calculated using R to determine the surface area under the curve of each individual nuclei and analyzed with Prism (Graphpad 8.0.1) using a Kruskal–Wallis test for significance with multiple comparison adjustments. The calibration method has a minimum detection threshold of >3 transcripts per transcription site. All figures report the number of nuclei and movies used for analysis in figure legends. Movies of genotypes not supplied as Supplementary Movies available upon request.
smiFISH
Embryos heterozygous for the transgene of interest were fixed in 10% formaldehyde/heptane for 25 min with shaking before a methanol quench and stored at −20 °C in methanol before use. smiFISH probes targeting the yellow (y) reporter gene were designed using Oligostan83,84,85 with FLAP-Y for secondary probe recognition (Integrated DNA Technologies, Inc.). Probe sequences are provided in Supplementary Table 4. Secondary probes were conjugated to Cy3 at the 5’ and 3’ ends (Integrated DNA Technologies, Inc.). Probes were resuspended in TE at appropriate equimolar concentrations. Prior to probe addition to embryos, the y targeting primary probes were hybridized to the secondary FLAP-Y probes as described83 and maintained at −20 °C in the dark prior to use.
Embryos were prepared for smiFISH as briefly follows: embryos were dehydrated with 2 × 5 min washes in 100% ethanol, followed by rehydration in PBT for 4 × 15 min and equilibration in 15% formamide/1 × SSC for 15 min. During equilibration, the smiFISH probe mixture was prepared with a final concentration of 1 × SSC, 0.34 µg µL−1 E. coli tRNA (New England Biolabs), 15% formamide (Sigma), 5-µL probe duplex, 0.2 µg µL−1 RNAse-free BSA, 2 mM vanadyl-ribonucleoside complex (New England Biolabs), and 10.6% dextran sulfate (Sigma) in RNAse-free water. The equilibration mixture was removed and replaced with probe mixture, and embryos were incubated overnight in the dark at 37 °C. The following day, embryos were rinsed twice in equilibration mix and twice in PBT, followed by DAPI staining and three PBT washes before mounting in ProLong Gold mounting media (Life Technologies).
Fixed sample imaging and analysis
Fixed sample imaging was performed on an LSM 880 with Airyscan module (Zeiss). Z-planes were acquired with 0.20 µm spacing to a typical depth of 80–100 Z-planes from the apical surface of the embryo, using laser scanning confocal in Airyscan super-resolution mode with a zoom of 3.0. DAPI excitation was performed with a 405-nm laser and Cy3 excitation with a 561-nm laser, with detection on a GaAsP-PMT array coupled to an Airyscan detector. Airyscan processing was performed using 3D Zen Black v3.2 (Zeiss) prior to analysis. Embryos were staged based on membrane invagination. Z-stacks were taken at both the center of the presumptive mesoderm as well as at the border region.
Images were analyzed (Supplementary Fig. 3) with Imaris v9.2.2 by first determining the threshold of detection using the non-mesodermal border region. After applying the threshold to the center pattern, fixed-sized shells (XY radius >0.3 µm, Z radius of >1.0 µm) were created around the centroid of detected objects. The median signal intensity of object shells was used as a proxy for the intensity of single molecules of RNA. The transcription site intensity of each nucleus was summed to account for the presence of sister chromatids and treated as a single transcription site throughout. The mean transcription site intensity was divided by the median single-molecule intensity to determine the average number of mRNA molecules present at the transcription sites.
Mathematical modeling of burst parameters and multiexponential regression fitting
Burst deconvolution and pol II positioning
The Pol II positions were found by combining a genetic algorithm with a local optimization procedure59.
Before the initiation of the analysis algorithm, several key parameters were established. The Pol II elongation speed was fixed at 45 bp s−1(see ref. 60). The reporter construct transcript was divided into three sections consisting of the pre-MS2 fragment (41 bp), 24xMS2 loops (1292 bp), and post-MS2 fragment containing the yellow reporter (4526 bp). The retention time was assumed to be small in relation to the time needed to produce a transcript, and so was fixed at 0 s. The temporal resolution of each movie was 3.86 s per frame. This frame rate is sufficient to detect processes that occur on the order of seconds.
The possible polymerase positions were discretized using a step of 30 bp (or equivalently 2/3 s). This step was chosen, as it is smaller than the minimum polymerase spacing and large enough to have a reasonable computation time. For a movie of 35-min length, this choice corresponds to a maximum number of 3150 positions.
The algorithm was implemented in Matlab R2020a using Global Optimization and Parallel Computing Toolboxes for optimizing Pol II positions in parallel for all nuclei in a collection of movies. The resulting positions are stored for analysis in the further steps of our computational pipeline. At this step, the density of Pol II initiation events can be visualized by binning time and checking the occurrence of Pol II activation in each bin. This was rendered as a heatmap in which rows represent a single-nucleus time series and the number of activation events per 30-s bin (or equivalently 1350 bp) is indicated by the color (Fig. 2e).
The deconvolution step is common to all of the MS2 data analysis pipelines. A detailed description of the algorithm can be found in ref. 59.
Multiexponential regression fitting of the survival function and model reverse engineering using the survival function
Data from several movies corresponding to the same genotype was first pooled together. The entry and exit of each trace corresponding to a unique nucleus were defined using a threshold representing 1/5 of the maximum intensity for the specific trace, in order to restrict the analysis to the stable part of the signal (Fig. 2a, gray box). Waiting times were extracted as differences between successive Pol II positions from all the resulting traces and the corresponding data was used to estimate the nonparametric cumulative distribution function by the Meyer–Kaplan method. This also permits the calculation of a 95% confidence interval for the experimental survival function that is further used to judge the quality of a parametric multiexponential regression fitting.
Then, a multiexponential regression fitting produced a set of 2N − 1 distribution parameters, where N is the number of exponentials in the regression procedure (3 for N = 2 and 5 for N = 3). The regression procedure was initiated with multiple initial guesses and followed by local gradient optimization of the following objective function59:
where \(S\left({t}_{i}\right){,S}_{e}\left({t}_{i}\right)\) are the theoretical (multiexponential) and empirical (estimated by the Meyer–Kaplan method) survival functions, respectively, and \(\alpha\) is a parameter satisfying \(0\le \alpha \le 1\) and representing the weight of linear scale differences in the objective function. We chose an intermediate value \(\alpha =0.6\) for all our parameter estimates (these estimates are nevertheless robust with respect to \(\alpha\)).
The optimization resulted in a best-fit solution with additional suboptimal solutions (local optima with objective function value larger than the best fit). A multiexponential regression is considered acceptable if the predicted survival function is within the confidence bounds of the experimental survival function. This provides a method to select the number N of exponentials in the regression: we progressively increase N starting with N = 2 until an acceptable regression is reached.
The 2N − 1 distribution parameters can be computed from the 2N − 1 kinetic parameters of a N-state transcriptional bursting model. Conversely, a symbolic solution for the inverse problem was obtained, allowing computation of the kinetic parameters from the distribution parameters and reverse engineering of the transcriptional bursting model. In particular, it is possible to know exactly when the inverse problem is well-posed, i.e., there is a unique solution in terms of kinetic parameters for any given distribution parameters in a domain.
The transcriptional bursting models used in this paper are as following (see Tantale et al.59 for a detailed description):
For N = 2, there were three distribution parameters and three kinetic parameters.
The distribution parameters are \({{A}}_{1},{\lambda }_{1}\), \({\lambda }_{2}\), defining the survival function
The solution of the inverse problem for the ON–OFF telegraph model (Fig. 2k and Fig. 6b) is
where \({k}_{2},{k}_{1}^{+},{k}_{1}^{-}\) are the initiation rate, the OFF to ON and ON to OFF transition rates, respectively. Thus, the duration of the OFF and ON states can be calculated as:
For this model, the probability to be in the state ON and OFF is:
For N = 3, there were five distribution parameters and five kinetic parameters59.
The distribution parameters are \({{A}}_{1},{{A}_{2},\lambda }_{1}\), \({\lambda }_{2},{\lambda }_{3}\), defining the survival function
The inverse problem has a unique solution for the three-state model (non-obligatory pause) with two OFF states (OFF and PAUSE) and one ON state (Fig. 6c)
where
and \({k}_{3},{k}_{2}^{+},{k}_{2}^{-},{k}_{1}^{+},{k}_{1}^{-}\)are the transcription initiation, PAUSE to ON, ON to PAUSE, OFF to ON, and ON to OFF rates, respectively.
Thus, duration of the ON, OFF, and PAUSE states can be calculated as:
For this model, the steady-state probability to be in a given promoter state is
The alternative three-state model with systematic pause (Supplementary Fig. 5) satisfies the following relation among distribution parameters59:
This means that only four and not five distribution parameters are free, which further constrains the three-exponential fittings. In order to infer this model, a constrained fitting was performed but the bad quality of fitting recommended rejection of the model (Supplementary Fig. 5).
Testing the method with artificial data
The entire computational pipeline was tested using artificial data (see also Tantale et al.59). Artificial traces were generated by simulating the model using the Gillespie algorithm with parameter sets similar to those identified from data. The simulations generated artificial polymerase positions, from which a first version of the signal was computed by convolution. The results are provided in Fig. 2e–h.
A modified Kolmogorov–Smirnov test for the parametric distribution
A one-sided one-sample Kolmogorov–Smirnov (KS) test was used to check if the waiting times follow the parametric fitted distribution. The output of this test is a P value that is large if the waiting times follow the fitted distribution, and small if not. The KS test is based on differences between estimated and empirical probabilities, being thus sensible to errors in the larger probabilities of shorter waiting times, and less sensitive to rare, very long waiting times. Our fitting procedure combines linear and logarithmic scales to find a balance between short and long time scales (see Tantale et al.59). Although relatively small (in terms of the objective function), the resulting errors on short times are large enough to be considered significant by the KS test for all models. In order to be able to distinguish between models, we have limited the analysis to times larger than 10–20 s.
Error intervals
Distribution parameters result from multiexponential regression fitting using gradient methods with multiple initial data. These optimization methods provide the best fit (global optimum) but also suboptimal parameter values. Using an overflow ratio (a number larger than one, in our case 2) to restrict the number of suboptimal solutions, we define boundaries of the error interval as the minimum and maximum parameter value compatible with an objective function less than the best-fit times the overflow.
Mathematical modeling of postmitotic gaps
The distribution of the postmitotic gap was estimated using the equation below. The fitting suggests that the timescale parameters of the postmitotic gap are even, i.e., \({\lambda }_{1}={\lambda }_{2}={\lambda }_{3}\). In this limit, the five parameters, three-exponential distribution defined by the equation below, becomes the simpler, three parameters mixed gamma distribution described as
where
is the lower incomplete gamma function, and \({p}_{1},{p}_{2},{b},\) are the probabilities of one step, two steps, and the mean step duration, respectively.
Chromatin immunoprecipitation
Homozygous embryos were collected and fixed in 1.4% formaldehyde for 20 min prior to dry storage at −80 °C. Embryos were dissociated on ice in RIPA buffer (150 mM NaCl, 1.0% IGEPAL CA-630, 0.5% sodium deoxycholate, 0.1% SDS, 50 mM Tris-HCl pH 7.8) supplemented with protease and phosphatase inhibitors (Roche). Chromatin shearing was performed in a pre-chilled Bioruptor Pico (Diagenode) for five cycles of 30 s/30 s. Samples were divided into equal volumes and incubated overnight at 4 °C with either rabbit anti-Rbp3 (10 μg) or rabbit anti-Nelf-E (10 μg) (gifts from J. Zeitlinger16) or normal serum. A second incubation overnight at 4 °C with Protein G-magnetic Sepharose beads (GE Healthcare) was used to pull down protein:DNA complexes prior to washing in low salt (120 mM NaCl, 0.1% SDS, 0.5% Triton X-100, 20 mM Tris-HCl pH 8.0) and high salt (500 mM NaCl, 0.1% SDS, 0.8% Triton X-100, 20 mM Tris-HCl pH 8.0) buffers, elution, incubation with protease K and RNAse A, and DNA retrieval using the QiaQuick PCR cleanup kit (Qiagen). qPCR analysis was performed using Light Cycler 480 SYBR Green I Master Mix (Roche) using primers listed in Supplementary Table 5. Analysis was performed using Microsoft Excel and Prism (Graphpad v8.0.1) with a Mann–Whitney test to determine significance.
Quantitative RT-PCR
To test RNAi-mediated knockdown of Nelf-A expression, 0–2 h embryos were homogenized in Trizol (Invitrogen) and RNA was extracted per the manufacturer’s directions. Reverse transcription was performed using the Superscript IV system (Invitrogen) with oligo d(T)20 priming prior to qPCR analysis. nos:GAL4-VP16; UAS:white RNAi embryos were used as the control and all measurements were performed in biological and technical triplicate. qPCR analysis was performed using Light Cycler 480 SYBR Green I Master Mix (Roche) using primers listed in Supplementary Table 5. Analysis was performed using Microsoft Excel and Prism (Graphpad 8.0.1) with a one-tailed Student’s t test to determine significance.
Image analysis software
Live imaging analysis software is available at http://www.igmm.cnrs.fr/segment-track/ along with a video tutorial. Code for the mathematical analysis of burst parameters and multiexponential regression fitting is available59.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Data availability
The data that support this study is available from the corresponding author upon reasonable request. Source data are provided with this paper.
References
Vo Ngoc, L., Kassavetis, G. A. & Kadonaga, J. T. The RNA polymerase II core promoter in Drosophila. Genetics 212, 13–24 (2019).
Bhuiyan, T. & Timmers, H. T. M. Promoter recognition: putting TFIID on the spot. Trends Cell Biol. 29, 752–763 (2019).
Cramer, P. Organization and regulation of gene transcription. Nature 573, 45–54 (2019).
Gupta, K., Sari-ak, D., Haffke, M., Trowitzsch, S. & Berger, I. Zooming in on transcription preinitiation. J. Mol. Biol. 428, 2581–2591 (2016).
Patel, A. B., Greber, B. J. & Nogales, E. Recent insights into the structure of TFIID, its assembly, and its binding to core promoter. Curr. Opin. Struct. Biol. 61, 17–24 (2020).
Haberle, V. & Stark, A. Eukaryotic core promoters and the functional basis of transcription initiation. Nat. Rev. Mol. Cell Biol. 19, 621–637 (2018).
Core, L. J. & Adelman, K. Promoter-proximal pausing of RNA polymerase II: a nexus of gene regulation. Genes Dev. 33, 960–982 (2019).
Fant, C. B. et al. TFIID enables RNA polymerase II promoter-proximal pausing. Mol. Cell 78, 785–793.e8 (2020).
Hendrix, D. A., Hong, J.-W., Zeitlinger, J., Rokhsar, D. S. & Levine, M. S. Promoter elements associated with RNA Pol II stalling in the Drosophila embryo. Proc. Natl Acad. Sci. USA 105, 7762–7767 (2008).
Lee, C. et al. NELF and GAGA factor are linked to promoter-proximal pausing at many genes in Drosophila. Mol. Cell. Biol. 28, 3290–3300 (2008).
Li, J. et al. Kinetic competition between elongation rate and binding of NELF controls promoter-proximal pausing. Mol. Cell 50, 711–722 (2013).
Shao, W., Alcantara, S. G.-M. & Zeitlinger, J. Reporter-ChIP-nexus reveals strong contribution of the Drosophila initiator sequence to RNA polymerase pausing. eLife 8, 1–25 (2019).
Gilchrist, D. A. et al. Pausing of RNA polymerase II disrupts DNA-specified nucleosome organization to enable precise gene regulation. Cell 143, 540–551 (2010).
Henriques, T. et al. Stable pausing by RNA polymerase II provides an opportunity to target and integrate regulatory signals. Mol. Cell 52, 517–528 (2013).
Krebs, A. R. et al. Genome-wide single-molecule footprinting reveals high RNA polymerase II turnover at paused promoters. Mol. Cell 67, 411–422.e4 (2017).
Shao, W. & Zeitlinger, J. Paused RNA polymerase II inhibits new transcriptional initiation. Nat. Genet. 49, 1045–1051 (2017).
Jonkers, I. & Lis, J. T. Getting up to speed with transcription elongation by RNA polymerase II. Nat. Rev. Mol. Cell Biol. 16, 167–177 (2015).
Saunders, A., Core, L. J., Sutcliffe, C., Lis, J. T. & Ashe, H. L. Extensive polymerase pausing during Drosophila axis patterning enables high-level and pliable transcription. Genes Dev. 27, 1146–1158 (2013).
Ghavi-Helm, Y. et al. Enhancer loops appear stable during development and are associated with paused polymerase. Nature 512, 96–100 (2014).
Lagha, M., Bothma, J. P. & Levine, M. Mechanisms of transcriptional precision in animal development. Trends Genet. 28, 409–416 (2012).
Lagha, M. et al. Paused Pol II coordinates tissue morphogenesis in the Drosophila embryo. Cell 153, 976–987 (2013).
Tunnacliffe, E. & Chubb, J. R. What is a transcriptional burst? Trends Genet. 36, 288–297 (2020).
Bertrand, E. et al. Localization of ASH1 mRNA particles in living yeast. Mol. Cell 2, 437–445 (1998).
Fusco, D. et al. Single mRNA molecules demonstrate probabilistic movement in living mammalian cells. Curr. Biol. 13, 161–167 (2003).
Chubb, J. R., Trcek, T., Shenoy, S. M. & Singer, R. H. Transcriptional pulsing of a developmental gene. Curr. Biol. 16, 1018–1025 (2006).
Golding, I., Paulsson, J., Zawilski, S. M. & Cox, E. C. Real-time kinetics of gene activity in individual bacteria. Cell 123, 1025–1036 (2005).
Larson, D. R., Zenklusen, D., Wu, B., Chao, J. A. & Singer, R. H. Real-time observation of transcription initiation and elongation on an endogenous yeast gene. Science 332, 475–478 (2011).
Garcia, H. G., Tikhonov, M., Lin, A. & Gregor, T. Quantitative imaging of transcription in living Drosophila embryos links polymerase activity to patterning. Curr. Biol. 23, 2140–2145 (2013).
Lucas, T. et al. Live imaging of bicoid-dependent transcription in Drosophila embryos. Curr. Biol. 23, 2135–2139 (2013).
Lammers, N. C. et al. Multimodal transcriptional control of pattern formation in embryonic development. Proc. Natl Acad. Sci. USA 117, 836–847 (2020).
Falo-Sanjuan, J., Lammers, N. C., Garcia, H. G. & Bray, S. J. Enhancer priming enables fast and sustained transcriptional responses to notch signaling. Dev. Cell 50, 411–425.e8 (2019).
Hoppe, C. et al. Modulation of the promoter activation rate dictates the transcriptional response to graded BMP signaling levels in the Drosophila embryo. Dev. Cell 54, 727–741.e7 (2020).
Pichon, X., Lagha, M., Mueller, F. & Bertrand, E. A growing toolbox to image gene expression in single cells: sensitive approaches for demanding challenges. Mol. Cell 71, 468–480 (2018).
Lammers, N. C., Kim, Y. J., Zhao, J. & Garcia, H. G. A matter of time: using dynamics and theory to uncover mechanisms of transcriptional bursting. Curr. Opin. Cell Biol. 67, 147–157 (2020).
Cisse, I. I. et al. Real-time dynamics of RNA polymerase II clustering in live human cells. Science 341, 664–667 (2013).
Corrigan, A. M., Tunnacliffe, E., Cannon, D. & Chubb, J. R. A continuum model of transcriptional bursting. eLife 5, 1–38 (2016).
Tantale, K. et al. A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting. Nat. Commun. 7, 12248 (2016).
Donovan, B. T. et al. Live-cell imaging reveals the interplay between transcription factors, nucleosomes, and bursting. EMBO J. 38, 1–18 (2019).
Teves, S. S. et al. A stable mode of bookmarking by TBP recruits RNA polymerase II to mitotic chromosomes. eLife 7, 1–22 (2018).
Nicolas, D., Zoller, B., Suter, D. M. & Naef, F. Modulation of transcriptional burst frequency by histone acetylation. Proc. Natl Acad. Sci. USA 115, 7153–7158 (2018).
Zoller, B., Nicolas, D., Molina, N. & Naef, F. Structure of silent transcription intervals and noise characteristics of mammalian genes. Mol. Syst. Biol. 11, 823 (2015).
Bartman, C. R. et al. Transcriptional burst initiation and polymerase pause release are key control points of transcriptional regulation. Mol. Cell 73, 519–532.e4 (2019).
Fukaya, T., Lim, B. & Levine, M. Enhancer control of transcriptional bursting. Cell 166, 358–368 (2016).
Bartman, C. R., Hsu, S. C., Hsiung, C. C. S., Raj, A. & Blobel, G. A. Enhancer regulation of transcriptional bursting parameters revealed by forced chromatin looping. Mol. Cell 62, 237–247 (2016).
Louder, R. K. et al. Structure of promoter-bound TFIID and model of human pre-initiation complex assembly. Nature 531, 604–609 (2016).
Patel, A. B. et al. Structure of human TFIID and mechanism of TBP loading onto promoter DNA. Science 362, eaau8872 (2018).
Gregor, T., Garcia, H. G. & Little, S. C. The embryo as a laboratory: quantifying transcription in Drosophila. Trends Genet. 30, 364–375 (2014).
Schulz, K. N. & Harrison, M. M. Mechanisms regulating zygotic genome activation. Nat. Rev. Genet. 20, 221–234 (2019).
Sloutskin, A. et al. ElemeNT: a computational tool for detecting core promoter elements. Transcription 6, 41–50 (2015).
Chen, K. et al. A global change in RNA polymerase II pausing during the Drosophila midblastula transition. eLife 2, e00861 (2013).
Ferraro, T. et al. Transcriptional memory in the Drosophila embryo. Curr. Biol. 26, 212–218 (2016).
Lécuyer, E. et al. Global analysis of mRNA localization reveals a prominent role in organizing cellular architecture and function. Cell 131, 174–187 (2007).
Wilk, R., Hu, J., Blotsky, D. & Krause, H. M. Diverse and pervasive subcellular distributions for both coding and long noncoding RNAs. Genes Dev. 30, 594–609 (2016).
Hoskins, R. A. et al. Genome-wide analysis of promoter architecture in Drosophila melanogaster. Genome Res. 21, 182–192 (2011).
Perry, M. W., Boettiger, A. N., Bothma, J. P. & Levine, M. Shadow enhancers foster robustness of Drosophila gastrulation. Curr. Biol. 20, 1562–1567 (2010).
Fernandez, C. & Lagha, M. Lighting up gene activation in living Drosophila embryos. Methods Mol. Biol. 2038, 63–74 (2019).
Dufourt, J. et al. Temporal control of gene expression by the pioneer factor Zelda through transient interactions in hubs. Nat. Commun. 9, 5194 (2018).
Boettiger, A. N. & Levine, M. Synchronous and stochastic patterns of gene activation in the Drosophila embryo. Science 325, 471–473 (2009).
Tantale, K. et al. Stochastic pausing at latent HIV-1 promoters generates transcriptional bursting. https://doi.org/10.1038/s41467-021-24462-5 (2021).
Fukaya, T., Lim, B. & Levine, M. Rapid rates of Pol II elongation in the Drosophila embryo. Curr. Biol. 27, 1387–1391 (2017).
Coulon, A. & Larson, D. R. Fluctuation analysis. Methods Enzymol. 159–191. https://doi.org/10.1016/bs.mie.2016.03.017 (2016).
Desponds, J. et al. Precision of readout at the hunchback gene: analyzing short transcription time traces in living fly embryos. PLoS Comput. Biol. 12, 1–31 (2016).
Ferguson, M. L. & Larson, D. R. Measuring transcription dynamics in living cells using fluctuation analysis. Methods Mol. Biol. 47–60. https://doi.org/10.1007/978-1-62703-526-2_4 (2013).
Rodriguez, J. et al. Intrinsic dynamics of a human gene reveal the basis of expression heterogeneity. Cell 176, 213–226.e18 (2019).
Larson, D. R., Singer, R. H. & Zenklusen, D. A single molecule view of gene expression. Trends Cell Biol. 19, 630–637 (2009).
Boettiger, A. N. & Levine, M. S. Rapid transcription fosters coordinate snail expression in the Drosophila embryo. Cell Rep. 3, 8–15 (2013).
Larsson, A. J. M. et al. Genomic encoding of transcriptional burst kinetics. Nature 565, 251–254 (2019).
Bressmann, T. Self-inflicted cosmetic tongue split: a case report. J. Can. Dent. Assoc. 70, 156–157 (2004).
Vos, S. M., Farnung, L., Urlaub, H. & Cramer, P. Structure of paused transcription complex Pol II-DSIF-NELF. Nature 560, 601–606 (2018).
Adelman, K. & Lis, J. T. Promoter-proximal pausing of RNA polymerase II: emerging roles in metazoans. Nat. Rev. Genet. 13, 720–731 (2012).
Nechaev, S. et al. Global analysis of short RNAs reveals widespread promoter-proximal stalling and arrest of Pol II in Drosophila. Science 327, 335–338 (2010).
Gressel, S. et al. CDK9-dependent RNA polymerase II pausing controls transcription initiation. eLife 6, 1–24 (2017).
Graveley, B. R. et al. The developmental transcriptome of Drosophila melanogaster. Nature 471, 473–479 (2011).
Zeitlinger, J. et al. RNA polymerase stalling at developmental control genes in the Drosophila melanogaster embryo. Nat. Genet. 39, 1512–1516 (2007).
Zoller, B., Little, S. C. & Gregor, T. Diverse spatial expression patterns emerge from unified kinetics of transcriptional bursting. Cell 175, 835–847.e25 (2018).
Weingarten-Gabbay, S. et al. Systematic interrogation of human promoters. Genome Res. 29, 171–183 (2019).
Haberle, V. et al. Transcriptional cofactors display specificity for distinct types of core promoters. Nature 570, 122–126 (2019).
Zabidi, M. A. et al. Enhancer–core-promoter specificity separates developmental and housekeeping gene regulation. Nature 518, 556–559 (2015).
Hasegawa, Y. & Struhl, K. Promoter-specific dynamics of TATA-binding protein association with the human genome. Genome Res. 29, 1939–1950 (2019).
Haberle, V. et al. Two independent transcription initiation codes overlap on vertebrate core promoters. Nature 507, 381–385 (2014).
Krueger, D. et al. Principles and applications of optogenetics in developmental biology. Development 146, dev175067 (2019).
Blumenthal, A. B., Kriegstein, H. J. & Hogness, D. S. The units of DNA replication in Drosophila melanogaster chromosomes. Cold Spring Harb. Symp. Quant. Biol. 38, 205–223 (1974).
Tsanov, N. et al. smiFISH and FISH-quant—a flexible single RNA detection approach with super-resolution capability. Nucleic Acids Res. 44, e165–e165 (2016).
Wang, Y.-L. et al. TRF2, but not TBP, mediates the transcription of ribosomal protein genes. Genes Dev. 28, 1550–1555 (2014).
Little, S. C., Tikhonov, M. & Gregor, T. Precise developmental gene expression arises from globally stochastic transcriptional activity. Cell 154, 789–800 (2013).
Acknowledgements
We thank Julia Zeitlinger, Jean-Christophe Andrau, Jeremy Dufourt, and all members of the Lagha lab for their critical reading of the manuscript and constructive discussions. We are grateful to F. Juge for the insightful discussion and H. Faure-Gautron for technical assistance. We acknowledge the MRI imaging facility, a member of the national infrastructure France-BioImaging supported by the French National Research Agency (ANR-10-INBS-04, «Investments for the future»). This work was supported by the ERC SyncDev starting grant to M.L and CNRS; by CNRS grant PEPS MIGHTY to O.R. and E.B., and by ANRS grants ECTZ62561 to E.B.. O.R. acknowledges support from the French National Research Agency (ANR-17-CE40-0036, project SYMBIONT). Machine-learning calculations were performed on the HPC facility MESO@LR run by the University of Montpellier.
Author information
Authors and Affiliations
Contributions
Conception and design: M.L.; acquisition of the data: M.D., V.P., and C.F.; analysis: M.L., M.D., V.P., and O.R.; software: A.T. and O.R.; modeling: O.R.; interpretation of the data: M.L., M.D., V.P., and O.R.; writing: V.P. and M.L. with input from O.R. and E.B.; visualization: M.D., V.P., and O.R.; supervision: M.L.; fund acquisition: M.L. O.R. and E.B. conceived and developed the machine-learning deconvolution procedure and the facultative model of paused polymerase. All authors read the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Peer review information Nature Communications thanks David Suter and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Source data
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Pimmett, V.L., Dejean, M., Fernandez, C. et al. Quantitative imaging of transcription in living Drosophila embryos reveals the impact of core promoter motifs on promoter state dynamics. Nat Commun 12, 4504 (2021). https://doi.org/10.1038/s41467-021-24461-6
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467-021-24461-6
This article is cited by
-
Tissue-specific RNA Polymerase II promoter-proximal pause release and burst kinetics in a Drosophila embryonic patterning network
Genome Biology (2024)
-
The Polycomb system sustains promoters in a deep OFF state by limiting pre-initiation complex formation to counteract transcription
Nature Cell Biology (2024)
-
Single-cell nascent RNA sequencing unveils coordinated global transcription
Nature (2024)
-
Lola-I is a promoter pioneer factor that establishes de novo Pol II pausing during development
Nature Communications (2023)
-
Dynamic epistasis analysis reveals how chromatin remodeling regulates transcriptional bursting
Nature Structural & Molecular Biology (2023)