# Common bursting relationships underlie eukaryotic transcription dynamics

Po-Ta Chen, Benjamin Zoller, Michal Levo, Thomas Gregor

► To cite this version:

Po-Ta Chen, Benjamin Zoller, Michal Levo, Thomas Gregor. Common bursting relationships underlie eukaryotic transcription dynamics. 2023. pasteur-04116784

HAL Id: pasteur-04116784

https://pasteur.hal.science/pasteur-04116784

Preprint submitted on 5 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Common bursting relationships underlie eukaryotic transcription dynamics

Po-Ta Chen,<sup>1,*</sup> Benjamin Zoller,<sup>1,2,*</sup> Michal Levo,<sup>1,*</sup> and Thomas Gregor<sup>1,2,†</sup>

<sup>1</sup>Joseph Henry Laboratories of Physics & Lewis-Sigler Institute for Integrative Genomics,
Princeton University, Princeton, NJ 08544, USA

<sup>2</sup>Department of Stem Cell and Developmental Biology, CNRS UMR3738 Paris Cité,
Institut Pasteur, 25 rue du Docteur Roux, 75015 Paris, France

(Dated: April 19, 2023)

Transcription commonly occurs in bursts resulting from alternating productive (ON) and quiescent (OFF) periods. Yet how transcriptional bursts are regulated to determine spatiotemporal transcriptional activity remains unclear. Here we perform live transcription imaging of key developmental genes in the fly embryo, with single polymerase sensitivity. Quantification of single allele transcription rates and multi-polymerase bursts reveals shared bursting relationships among all genes, across time and space, as well as *cis-* and *trans*-perturbations. We identify the allele's ON-probability as the main determinant of the transcription rate, while changes in the transcription initiation rate are limited. Any given ON-probability determines a specific combination of mean ON and OFF times, preserving a constant characteristic bursting time scale. Our findings point to a convergence of various regulatory processes that predominantly affect the ON-probability, thereby controlling mRNA production rather than mechanism-specific modulation of ON and OFF times. Our results thus motivate and guide new investigations into the mechanisms implementing these bursting rules and governing transcriptional regulation.

INTRODUCTION

Eukaryotic transcriptional regulation is an inherently
dynamic and stochastic process. Multiple molecular
events orchestrate, in space and time, the initiation of
productive transcription by individual RNA polymerases
(Pol II complexes), leading to the synthesis of nascent
RNA [1, 2]. The amount of transcribed mRNA molecules
in turn shapes protein production and thereby dictates
cellular behavior. Studies across various systems, from
yeast to mammalian cells, have revealed that transcrip-
tion occurs in bursts, namely the release of multiple Pol
IIs in what is often referred to as an ON period, followed
by a quiescent OFF period [3–8]. Yet it remains unclear
how the kinetic parameters of transcriptional bursting
determine mRNA production and govern spatiotemporal
transcription dynamics. Is the transcription rate con-
trolled primarily by tuning the durations of ON or OFF
periods, the initiation rate (the rate of Pol II release
during active periods), or by a combination of these?
Furthermore, are distinct bursting parameters controlled
by specific regulatory processes, and are distinct burst-
ing strategies underlying temporal versus spatial (tissue-
specific) control of transcription?

As a multitude of molecular processes is known to influ-
ence transcriptional activity, several studies aimed to un-
cover links between regulatory determinants and param-
eters of transcriptional bursting [9–11]. Transcriptional
bursting is commonly described by parameters such as
burst size, burst frequency, the kinetic rates governing
ON and OFF times as well as transcription initiation

[3, 12, 13]. Regulatory determinants, such as transcrip-
tion factor (TF) binding, *cis*-regulatory elements, nucle-
osome occupancy, histone modification, and enhancer-
promoter interactions were suggested to affect distinct
bursting parameters [14–22]. Yet, it is difficult to inte-
grate these observations and form a unified understand-
ing of transcriptional control via bursting dynamics.

Much of our quantitative knowledge about transcrip-
tional bursting heavily relies on fixed data [3, 14, 23–26].
Capturing transcriptional bursts in vivo across space and
time remains challenging [20, 27–29]. Adding to this chal-
lenge, live measurements need to be quantifiable in ab-
solute units (i.e., mRNA count) to facilitate comparisons
between different genes and conditions [30–32]. More-
over, to understand the entire spectrum of bursting dy-
namics, there is a need to probe the full dynamic range
of a gene’s activity [25]. Measurements in an endogenous
system where tightly regulated spatiotemporal transcrip-
tional dynamics dictate cell-fate determination can fur-
ther elucidate the functional consequences of bursts and
their relation to different regulatory determinants. The
early *Drosophila* embryo provides a unique system that
meets all these requirements [33].

Here we quantify the endogenous transcription dynam-
ics of key developmental genes in living *Drosophila* em-
bryos. We identify a single control parameter, the in-
stantaneous ON-probability of an allele as the dominant
determinant of transcriptional activity, while the initia-
tion rate is mostly conserved. This finding holds across
spatial domains, developmental times, genes, and pertur-
bations of *cis*-regulatory elements and *trans*-regulators.
Surprisingly, we find a largely constant ON-OFF switch-
ing correlation time of roughly one minute. A corollary
of the latter is that mean ON and OFF times are tightly
coupled, regardless of the *trans*-environment or the *cis*-
regulatory architecture. While perturbations in the up-

<footnote>* These authors contributed equally.</footnote>

<footnote>† Correspondence: tg2@princeton.edu</footnote>
stream regulatory processes lead to dramatic changes in the ON-probability, i.e., spatiotemporal changes in transcription rate, the underlying changes in ON and OFF periods are predicted from wild-type. Instead of a particular perturbation type dictating changes in specific bursting parameters, we observe that more generally lowly transcribing alleles are tuned to higher expression levels by increasing burst frequency, while highly transcribing alleles are mostly tuned by increasing burst size. These results imply that for the examined genes, the bursting phenomenon can be quantitatively understood by a few simple rules: two of the transcription parameters are quasi-constant, and all others are determined by the ON-probability. Hence, future investigations necessitate a re-examination of our mechanistic understanding of transcription, focusing on how regulatory processes influence a unique control parameter.

## RESULTS

**Instantaneous single allele transcription rate measurements.** To study the principles that govern transcription dynamics across space and time, we need access to the endogenous bursting kinetics at a single allele level. We designed an approach to obtain such quantitative measurements in living *Drosophila* embryos [31, 35]. A versatile CRISPR-based scheme is employed to incorporate MS2 cassettes into the gap genes' introns (or 3'UTR) [36]. These form stem-loops in the transcribed nascent RNA that are subsequently bound by fluorescent coat-proteins (Fig. 1A, S1A and Methods) [37, 38]. A custom-built two-photon microscope generates fluorescence images, capturing RNA synthesis at one tagged allele per nucleus with a 68-fold signal improvement over previous studies [31], approaching single-mRNA sensitivity (Fig. S1D-E). An optimized field-of-view yields 10 s interval time-lapses for hundreds of nuclei per embryo, during nuclear cycles (NC) 13 and 14 (Fig. 1A-B; Videos V1-V4), essential for statistical analysis.

To achieve a fully quantitative characterization, we calibrate our measurements to absolute units. We convert the fluorescence signal at the site of transcription into equivalent cytoplasmic mRNA units (C.U.) by matching the mean transcriptional activity to previously calibrated smFISH measurements (Fig. 1B, see Methods) [25, 34]. We find high agreement between smFISH and live measurements, as a single conversion factor adjusts for the difference in fluorescence signal between the two methods, with an average relative error of ~ 5% across all gap genes and MS2 insertion sites (Fig. 1C and S1B). This agreement extends to higher moments (Fig. S1C), despite the fully orthogonal nature of these techniques: one being non-invasive genetically but involving fixation, while the other involves gene editing and stem-loop cassette insertions. This strongly suggests that our live approach captures the endogenous situation and provides

means to express our dynamic transcription measurements in terms of absolute mRNA counts.

Our unique combination of absolute calibration and near single transcript sensitivity (Fig. S1D-E) allows us to reconstruct the underlying single allele transcription initiation events by individual Pol II, namely the events in which Pol II complexes are released onto the gene and engage in productive elongation. To infer these initiation events for each time series, we adopt a Baysian deconvolution approach that accounts for measurement noise (Fig. S1D and Methods). The convolution kernel models the fluorescent signal resulting from the Pol II elongation process through the stem-loop cassette (with constant and deterministic elongation, Fig. 1D) [20, 27]. For each time series, the approach generates multiple configurations of transcription initiation events (Fig. 1E). Averaging these configurations gives us a time-dependent instantaneous single allele transcription rate $r(t)$ per time series (Fig. 1F).

We validate this kernel-based deconvolution approach by performing dual-color tagging of the gene body (a 5' proximal intron and a 3'UTR tag). These measurements support our key assumptions (see Methods) and allow us to extract a Pol II elongation rate of $K_{elo} = 1.8 \pm 0.1$ kb/min, which is in line with previous measurements [31, 39] (Fig. S2). Our inferred transcription rates are thus no longer masked by the Pol II elongation dwell time, unlike the directly measured intensities of transcriptional activity. Transcription rates are thus independent of gene length, enabling the direct comparison between different genes and opening a path to identify common principles underlying transcription dynamics.

**Single allele transcription rates hint at a universal bursting regime.** Before analyzing the gap genes' transcription dynamics at the single allele level, we sought to determine whether the averages of the deconvolved single allele transcription rates $r$ recapitulate the well-documented average protein dynamics [40]. We compute a mean transcription rate $R = \langle r \rangle$ per gene along the anterior-posterior (AP) axis for all time points during NC13 and NC14 (Video V5); averaging occurs over 200-300 nuclei (each contributing one allele) in the same AP and time bin from 10-20 embryos (Fig. 1A and 1B). The extracted mean transcription rate profiles (Fig. 2A and S3A) strongly resemble the ones reconstructed from carefully staged gap gene antibody staining, including the well-documented posterior shifts during NC14. Indeed, with simple assumptions on diffusion and lifetime, the mean transcription rates $R$ predict protein patterns with minimal post-transcriptional regulation (Fig. S4, Video V6). Thus, in this system, the rules governing transcription rate modulation will largely determine function, i.e., protein synthesis.

While the examined genes exhibit a common range of mean transcription rates $R$ (Video V5), they display distinct spatiotemporal profiles, giving rise to gene-specific protein patterns. Yet when we examine the distributions
<figure><img src="image_2.png" /></figure>

FIG. 1. **Live single-cell transcription rate measurements of endogenous gap genes.** (A) Transcriptional activity measured of a single gap gene allele using a custom-built two-photon microscope in a living fly embryo. An 24xMS2 stem-loop cassette is inserted in the first intron of the gap genes. Constitutively expressed MCP-GFP binds stem-loops formed on nascent transcripts, making transcription sites appear bright above background (green hotspots) and enabling quantification of transcriptional activity at single allele resolution along the anterior-posterior (AP) axis of the embryo. (For genomic strategy of all gap gene loci see Fig. S1A). (B) Example single allele transcription time series for the gene *hunchback* across nuclear cycles NC13 and NC14 (sampled every 10 s) from a single AP bin (width ~ 2% embryo egg length *L*) at *x/L* = 0.435 ± 0.010. Intrinsically low embryo-to-embryo variability (compared to the total variance in our data, Fig. S1F) facilitated by the *Drosophila* system allows for pooling alleles from multiple spatially and temporally aligned embryos (*n* = 10–20). Mean transcriptional activity (black line) obtained from pooling 200–350 alleles. (C) Calibration of transcriptional activity in absolute units performed by matching mean activity profiles in a 5-min-interval during NC13 (gray shade in B) from live (color) and previously calibrated fixed smFISH (black) measurements for all examined gap genes [25]. A global conversion factor (one for all genes) leads to a match within 5% error between live and fixed profiles (Fig. S1B), resulting in our activity unit, i.e., cytoplasmic unit (C.U.) equivalent to the intensity of a fully elongated transcript [34]. (D) Reconstruction of transcription initiation events from deconvolution of single allele transcription time series. The signal is modeled as a convolution between transcription initiation events and a kernel accounting for the elongation of a single Pol II through the MS2 cassette and the gene body (using an elongation rate *K*<sub>elo</sub> = 1.8 kb/min, Fig. S2). Bayesian deconvolution is performed by sampling from the posterior distribution of possible configuration of initiation events given the measured activity and measurement noise (Fig. S1D-E). (E) Example deconvolved initiation configuration (gray bars) and corresponding reconstructed signal (red) from a single allele transcription time series (black). (F) single allele transcription rate (gray) from same allele as in E (black). The rate is estimated by counting the number of initiation events within 10 s intervals for a given sampled configuration and averaged over 1'000 of such configurations. The displayed solid line and envelope for transcription rate (gray) and reconstructed signal (red) correspond to the mean and one standard deviation of the posterior distribution.

of single allele transcription rates, *r*, underlying a similar mean transcription rate, *P*(*r*|*R*), we find that these distributions collapse across genes (Fig. 2B). Strikingly, for low- to mid-levels of *R*, the underlying distributions differ starkly from a constitutive regime, which would result in a Poisson distribution. At these levels, the large amounts of non-transcribing or barely transcribing alleles hint at quiescent OFF periods deviating from a constitutive regime. As the mean transcriptional activity increases, the distributions become more Poissonian, sug-

gesting that $P_{\text{ON}}$, i.e., the probability of the genes being ON, increases. These observations are consistent with bursting behavior, where the gene alternates between ON and OFF states.

The collapse of our data and the deviation from a constitutive, Poisson regime, are readily observed also when we compute the relationships between *R* and the higher moments of the $P(r|R)$ distributions (Fig. 2C and S3C). Our data approaches the Poissonian regime only on the extreme ends of the *R* spectrum, implying that the gap
<figure><img src="image_3.png" /></figure>

FIG. 2. Transcription rates display signatures of a universal bursting regime. (A) Snapshot of mean transcription rates *R* as a function of AP position in early NC14 (*t* = 15 min after mitosis) for different gap genes (color). AP profiles are obtained by averaging the deconvolved single allele transcription rates over ~ 200 nuclei within each AP bin. The black dashed lines correspond to the mean activity (Fig. 1C) of each gap gene at the same position and time normalized by the effective elongation time (see Methods, Fig. S3A and Fig. S4A). Tight agreement of colored and dashed profiles supports the deconvolution approach (error bars are one standard deviation across the means of 10–20 embryos). (B) Distribution *P*(*r*|*R*) of single allele transcription rates estimated within 1-min-intervals in both NC13 (color dotted lines) and early NC14 (color solid lines). These distributions are computed over all the nuclei from time points and AP bins whose mean transcription rate *R* is either in a low [2.1, 3.2], mid [7.5, 8.5] or high regime [12.8, 13.9] (as gray shade in A and C). The various gap gene distributions collapse at all regimes and differ from the Poisson distribution (black dashed line), suggestive of a universal bursting regime. (C) Variance of single allele transcription rates as a function of mean transcription rate *R* in both NC13 (square) and early NC14 (circle) (estimated over 1-min-intervals). Note a strong departure of the variance from a constitutive Poisson regime (dashed line, $\sigma^2 \sim R$). All gap genes follow the same trend suggestive of a common bursting regime. For higher moments see Fig. S3C. Vertical gray bars correspond to low, mid, and high *R*, as in A. (D) Auto-correlation (AC) functions of single allele transcription rates estimated within 10 s intervals for *hb* in early NC14 and averaged across time, and within a given AP bin across alleles. Color code stands for position along the AP axis. The AC functions are normalized by the variance and highlight an uncorrelated and a time-correlated component in the single-cell transcription rate fluctuations. The correlated component is characterized by a magnitude $\sigma_{AC}$ and an exponential decay with time scale $\tau_{AC}$. Such correlated fluctuations are expected to arise in a bursting regime due to ON-OFF gene switching (Fig. S3D). (E) Magnitude $\Sigma_{AC}$ of the correlated fluctuations in the single allele transcription rate as a function of mean transcription rate *R*. All gap data (color) collapses showing a universal trend (dashed line, guide to the eye). The fraction of correlated variability decreases as *R* increases, as expected when approaching a constitutive regime of uncorrelated Poisson initiation (Fig. S3E). (F) Correlation time of the correlated fluctuations in the single allele transcription rate as a function of mean transcription rate *R*. The correlation times are obtained by fitting exponentials to the correlated component of the AC function as in D. Strikingly, the correlation time is mostly conserved across genes and transcription levels. Error bars are bootstrapped 68% confidence interval.
genes transition all the way from fully OFF ($P_{ON} = 0$) to fully ON ($P_{ON} = 1$). The gap genes thus provide an opportunity to investigate how an underlying bursting regime can account for a full dynamic range of transcriptional activities.

The dynamic nature of our measurement allows us to
examine single allele transcription rates not only via dis-
tributions pooled across nuclei but also within individual
transcription time series. We find the auto-correlation
function of the single allele transcription rates provides
further evidence for an underlying common bursting
regime (Fig. 2D). An initial sharp drop of magnitude
$1 - \Sigma_{AC}$ at our sampling time scale ($\sim 10$ s) indicates
the presence of uncorrelated noise, consistent with inde-
pendent Pol II initiation events. This drop is followed
by a longer decay of correlated noise at time scale $\tau_{AC}$.
Such correlated noise is expected in a bursting regime,
as the switching between ON and OFF states introduces
temporal correlation in transcriptional activity. A the-
oretical and computational analysis using the two-state
model of transcription [12] supports both the interpre-
tation of the auto-correlation functions (Fig. S3D), and
our ability to estimate $\Sigma_{AC}$ and $\tau_{AC}$ properly from the
deconvolved rates (Fig. S3E-G).

We find both the magnitude of correlated noise Σ<sub>AC</sub>
(Fig. 2E) and the correlation time τ<sub>AC</sub> (Fig. 2F) col-
lapse across AP bins and different genes. Notably, the
magnitude Σ<sub>AC</sub> is highly constrained and drops at high
R, consistent with the behavior of the variance (Fig.
2C). Furthermore, the correlation time τ<sub>AC</sub> is largely con-
served across nuclear cycles, across AP bins, and across
genes, confined within the range of 1–2 min and averag-
ing to a value of 1.37 ± 0.31 min (Fig. 2F). This sur-
prising invariance of τ<sub>AC</sub> with respect to R suggests that
a key temporal characteristic of transcription dynamics
is highly conserved. Thus, both the static moment and
correlation-based analyses point to a common bursting
regime, applicable across genes and space, motivating a
time-dependent analysis of transcriptional bursts at the
level of individual time series.

**Allele ON-probability is the key regulated transcriptional parameter.** To directly quantify individual transcriptional bursts at the single allele level, we take advantage of our deconvolved initiation events and instantaneous transcription rate time series that are unencumbered by the signal-blurring effect of elongation (Fig. 3A). They allow us to identify distinct periods of active transcription, characterized by consecutive initiation events (i.e., multiple Pol IIs released into productive elongation), interpreted as ON periods, followed by quiescent periods, namely OFF periods (Fig. 3A). We define the switch of an allele from an OFF to an ON state when a moving average of the single allele transcription rate exceeds 2 mRNA/min (see Methods). This threshold is consistent with our detection sensitivity of 1–2 mRNAs, and the size of the moving window for averaging is set based on the correlation time scale from the

auto-correlation analysis. The main strength of our burst
calling routine is its sole reliance on a minimal cluster-
ing model, and as such being devoid of any mechanistic
assumptions on the underlying bursts (no explicit mech-
anistic model is needed, see Methods).

Given a computed bursting profile (demarcated ON and OFF periods) for every single allele, we can now ask how bursting dynamics underlie transcription rates. Specifically, the mean transcription rate at time *t*, *R*(*t*), can be decomposed into two parameters: the instantaneous probability of an allele being in the ON state *P*<sub>ON</sub>(*t*) (i.e., the fraction of ON alleles) and the mean initiation rate in the ON state *K*(*t*). Starting with the gene *hb*, we thus estimate, for a given AP bin, the time-dependent parameters *R*(*t*) and *P*<sub>ON</sub>(*t*) by averaging all (~ 250) single allele instantaneous transcription rates and counting the fraction of alleles in the ON state at time *t*, respectively (Fig. 3B-C). To compute *K*(*t*), we average initiation events restricted to the ON state (as opposed to *R*, which is averaging initiation events regardless of allele state). We repeat this procedure for each position of the AP-axis to obtain the full spatiotemporal dependence (Fig. 3D-F). We validate our approach for burst calling and the recovery of bursting parameters from transcription time series on simulated data (based on a 2-state model) with an overall median error of 10% (see Methods, Fig. S5 and Fig. S6).

All three parameters vary significantly across both
space and time (Fig. 3D-F). However, given that these
are related by $R = K \cdot P_{\text{ON}}$ (Fig. S7A), it is possible
that most of the variation in *R* stems from changes to
either *K*, or $P_{\text{ON}}$, or both. When *R* is plotted against
$P_{\text{ON}}$, all data points across time and space collapse in a
tight monotonically increasing function (Fig. 3G), with
$P_{\text{ON}}$ spanning from 0 to 1 (fully OFF to fully ON), echo-
ing back to the noise analysis above (Fig. 2B-C). Simi-
larly, the initiation rates *K* tightly collapse across space
and time when plotted against $P_{\text{ON}}$. However, *K* only
covers a two-fold change in dynamic range (Fig. 3H),
which is marginal compared to *R* spanning from 0 to 15
(mRNA/min) (Fig. 3I and S7B-C). This two-fold change
is largely due to the existence of two optically unresolved
sister chromatids and a modest time dependence of *K*
throughout the nuclear cycle (Fig. S7D-I, see Methods).
With these considerations we estimate the mean Pol II
spacing for a single active chromatid at 303±73 bp, con-
sistent with the classic Miller spreads with average Pol
II spacing of 330 ± 180 bp [41].

Overall, we find that *R* is tightly controlled by *P*<sub>ON</sub>, while *K* is only moderately modulated and has significantly less predictive power over *R*. These results for *hb* suggest that transcriptional activity is mainly controlled through the probability of an allele being in the ON state. Once in the ON state, transcription initiates at a quasi-constant rate.

**Constant switching correlation time restricts ON and OFF periods.** Given the central role of $P_{\text{ON}}$ in con-
<figure><img src="image_4.png" /><figcaption>FIG. 3. Direct estimation of instantaneous mean transcriptional parameters. (A) A simple clustering procedure on single allele transcription rate determines bursts of transcription. The rate (gray curve) is estimated using a moving window of width ~1 min over the deconvolved initiation events (gray vertical bar). A threshold at two mRNA/min (black dashed line) applied on the rate identifies individual bursts (blue curve). (B-C) Heatmaps of deconvolved single allele transcription rates (B) (estimated over 10 s intervals) and of corresponding ON-OFF periods (C) (obtained from burst calling) as a function of time during NC14 for N = 217 nuclei expressing *hb*-MS2 at AP position *x*/L = 0.43. Instantaneous mean transcription parameters such as transcription rate *R* (B, bottom) and ON-probability *P*<sub>ON</sub> (C, bottom) are obtained by vertically averaging the heatmaps (top) over all nuclei, respectively. (D-F) *hb* transcription rate *R* (D), ON-probability *P*<sub>ON</sub> (E), and initiation rate *K* (F) as a function of time in NC14 for all AP positions (color coded). *R* and *P*<sub>ON</sub> are computed as in B and C, respectively, and *K* is obtained by averaging the single allele transcription rate (B) conditioned on the locus being ON (C) over all nuclei in each AP bin. (For tests of burst calling procedure on simulated data see Figure S3.) (G-H) Transcription rate *R* (G) and initiation rate *K* (H) as a function of *P*<sub>ON</sub>, for all time points and positions, demonstrating a massive data collapse, suggesting that *P*<sub>ON</sub> is the central regulatory parameter for transcriptional bursting. (I) Transcription rate *R* as a function of both controlling parameter *P*<sub>ON</sub> and *K* in log-space. Since log (*R*) = log (*K*) + log (*P*<sub>ON</sub>) by construction, changes in *P*<sub>ON</sub> determine changes in *R* below the dashed line (*R* ~ 8.5 mRNA/min, corresponding to *P*<sub>ON</sub> = 0.75).</figcaption></figure>

trolling *hb*'s transcription rate, we aim to examine how
it decomposes into ON and OFF periods in individual
alleles. To this end we compute the mean ON and OFF
times, averaged across alleles of the same AP bin, and at
a given time (T<sub>ON</sub> and T<sub>OFF</sub>, Fig. 4A-C; see Methods).
Near steady state, we expect the ON-probability P<sub>ON</sub> to

be given by the ratio of T<sub>ON</sub> and T<sub>ON</sub> + T<sub>OFF</sub>. We ver-
ify this relationship (Fig. 4D), showing good agreement
(i.e., beyond an initial 7.5 min transient post mitosis, Fig.
S7J-K). This is a strong indication of the self-consistency
of our general approach for extracting these bursting pa-
rameters. Moreover, despite temporal changes in P<sub>ON</sub>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Allele ON-probability controls ON and OFF times. (A) Binarized heatmap from Fig. 3C. Instantaneous mean OFF-time $T_{OFF}$ (bottom, gray) and mean ON-time $T_{ON}$ (bottom, blue) are obtained by the weighted average of the ON and OFF times over all nuclei (see methods). The weights are given by the inverse of the number of time points within each period. (B-C) Mean OFF-time $T_{OFF}$ and mean ON-time $T_{ON}$ as a function of time and position (color coded) for *hb* in NC14. (D) The ratio of $T_{ON}$ over the sum of $T_{ON}$ and $T_{OFF}$ versus ON-probability $P_{ON}$ for all positions and time points beyond the 7.5 min mark in B and C (near-steady state both quantities are expected to be equal after initial transient, see Fig. S7J-K). Thus, temporal changes in transcriptional parameters must be slow enough to allow relaxation. (E-F) Mean OFF-time $T_{OFF}$ (E) and mean ON-time $T_{ON}$ (F) as a function of $P_{ON}$, for all positions and time points beyond the 7.5 min mark in B and C. (G) Effective switching correlation time $T_C$ (defined as: $1/T_C = 1/T_{ON} + 1/T_{OFF}$) as a function of $P_{ON}$, computed using data points in E and F. $T_C$ is mostly conserved across time points and position and is $P_{ON}$ independent.</figcaption></figure>

due to developmental regulation, transcriptional bursting in this system seems to operate in a near-steady-state regime.

While $T_{ON}$ and $T_{OFF}$ change over time and in different AP bins (Fig. 4B-C), when we plot $T_{ON}$ and $T_{OFF}$ against $P_{ON}$ all data points collapse again across time and space onto two tight anti-symmetric relationships (Fig. 4E-F). Various combinations of $T_{ON}$ and $T_{OFF}$ could potentially give rise to any given $P_{ON}$, however, here we observe a highly restricted range for these mean durations. Hence a given $P_{ON}$ is unequivocally linked to a specific pair of $T_{ON}$ and $T_{OFF}$, regardless of space and time.

An allele switching dynamically between ON and OFF states will have a correlation time $T_C$, which determines, on average, the time needed for the single allele transcription rate to become uncorrelated. For such a system, $T_C$ can be computed directly from the mean ON and OFF times and is defined by $1/T_C = 1/T_{ON} + 1/T_{OFF}$ (Fig. 4G, see Methods). Surprisingly, for all AP and time bins, $T_C$ is confined between 1-1.5 min and thus largely constant and independent of $P_{ON}$. Moreover, $T_C$ matches quantitatively the correlation time $\tau_{AC}$ from the autocorrelation analysis (Fig. 2F), which we found to be independent of the transcription rate. Given that $T_C$ characterizes allele switching correlations by construction, this

match suggests that the nature of $\tau_{AC}$ could indeed be related to bursting.

The fact that $T_C$ seems to remain conserved across space and time restricts the mean ON and OFF times. Indeed, $T_{ON}$ and $T_{OFF}$ can be expressed as a function of $P_{ON}$ and $T_C$ (via $T_{ON} = T_C/1 - P_{ON}$ and $T_{OFF} = T_C/P_{ON}$, since $P_{ON} = T_{ON}/(T_{ON}+T_{OFF})$ near the steady state, cf. Fig. 4D). Thus, the constancy of $T_C$ mathematically explains the tight anti-symmetric relationships of $T_{ON}$ and $T_{OFF}$ with respect to $P_{ON}$ (Fig. 4E-F), so that $P_{ON}$ not only governs the mean transcription rate $R$ but also the entire transcriptional bursting dynamics.

**Common bursting relationships underlie the regulation of all gap genes.** Should we expect these bursting parameter relationships found for *hb* to generalize to other gap genes or are they gene-specific? The gap genes differ in their *cis*-regulatory elements, namely different numbers and arrangement of enhancers and promoters (Fig. S1A), and different compositions of transcription factors binding sites within each enhancer. Correspondingly, and as discussed above, the gap genes display distinct transcriptional activities along the body axis (Fig. 1C) and across time (Video V8), which, in the case of *hb* we found to be largely governed by $P_{ON}$.

When we apply our deconvolution and burst calling
<figure><img src="image_6.png" /><figcaption>FIG. 5. Transcriptional parameters collapse for all gap genes. (A) Kymographs of ON-probability $P_{ON}$ for all gap genes as a function of position and time for NC13 and NC14. The spatiotemporal transcriptional pattern of the gap genes arises from a complex regulation of $P_{ON}$ (color map). (B-E) Transcriptional parameters collapse for all gap genes across time and position. Transcription rate $R$ (B), Mean OFF- (C) and ON-time (D) ($T_{OFF}$ and $T_{ON}$, respectively), and switching correlation time $T_C$ (E) as a function of the ON-probability $P_{ON}$. Colored data points represent individual gap gene (same color code as in A, see Fig. S8A-B for *gt* male data); underlying is the remaining data of all other genes (gray). (F) Density of all data points across space and time (NC13 and NC14) of transcriptional parameters for all gap genes. Color code represents the kernel density estimate of the parameters conditioned on $P_{ON}$ and normalized by the maximum density. Putative accessible space (gray shade) for plausible ranges of $K$ (0.1–30 mRNA/min) and $T_C$ (0.1–10 min). Due to the constancy of $T_C$ and the small changes in $K$ (2-fold at most), $P_{ON}$ almost fully determines $R$ and dictates the combinations of $T_{OFF}$ and $T_{ON}$. For $T_{OFF}$ and $T_{ON}$, the dashed lines stand for the 2-state model prediction based on $T_C$, and the solid lines take the finite recording-length into account (see Figure S8D).</figcaption></figure>
procedure (Fig. 3A) to the measured single allele tran-
scription time series of other gap genes (*gt*, *Kr*, and *kni*),
we observe that all genes differ substantially in their spa-
tiotemporal *P*<sub>ON</sub> profiles (Fig. 5A). These differences
are expected and reflective of the above-mentioned dis-
tinct underlying *cis* architectures and *trans* regulators of
these genes. Strikingly, despite these differences, they
show an almost identical mean transcriptional rate *R* to
*P*<sub>ON</sub> relationship (and *K* to *P*<sub>ON</sub> relationship, Fig. S8C),
substantiating *P*<sub>ON</sub> as the governing factor for the tran-
scriptional activity not only across time and space but
also across genes (Fig. 5B).

Computing the various bursting parameters for all gap
genes and plotting these as a function of P<sub>ON</sub>, we find
the genes display the same P<sub>ON</sub>-dependent relationships:
all genes share common T<sub>OFF</sub> to P<sub>ON</sub> (Fig. 5C) and
T<sub>ON</sub> to P<sub>ON</sub> (Fig. 5D) relationships. Thus, when differ-
ent genes display a specific P<sub>ON</sub> value, possibly at differ-
ent spatiotemporal coordinates, the underlying T<sub>ON</sub> and
T<sub>OFF</sub> periods employed are nonetheless largely similar.
This finding can be related to the switching correlation
time T<sub>C</sub>, which we find to be conserved across genes (Fig.
5E). The average T<sub>C</sub> value across all genes, positions, and
times is 1.25±0.37 min, very close to the prediction from
the single rate auto-correlation analysis (Fig. 2F).

Pooling all our data across all genes, times, locations,
and embryos (N > 10⁶ data points) and plotting each
of the computed bursting parameters (R, K, T<sub>C</sub>, T<sub>OFF</sub>
and T<sub>ON</sub>) against P<sub>ON</sub> (Fig. 5F) – including the often-
used burst frequency F = 1/(T<sub>ON</sub>+T<sub>OFF</sub>) and burst size
B = K · T<sub>ON</sub> (Fig. S8E) – reveals highly constrained
relationships. Indeed, all data points occupy only a very
small subset of the parameter space (see methods). Tight
functions can be mapped out (black lines) that confirm
the exceptional predictive power conferred by P<sub>ON</sub>. Sep-
arating all the data into three developmental time win-
dows (NC13 & early NC14, mid NC14, and late NC14),
shows a further tightening of the relationships in these
developmental stages (Fig. S9). This separation thus
confirms that parts of the dispersion results from slow
and moderate changes in K (and to a lesser extent T<sub>C</sub>)
over developmental time (at most 40% reduction of K
and 25% for T<sub>C</sub>).

All measured bursting dynamics seem to adhere to this
simple set of rules, across genes, space, and time. The
mean transcription rate *R* is essentially dictated by *P*<sub>ON</sub>,
with a largely constant *K*. A near-constant switching
correlation time *T*<sub>C</sub> leads to a specific functional relation-
ship for *T*<sub>ON</sub> and *T*<sub>OFF</sub> of inverse proportionality. Thus,
only one of these two parameters is principally modu-
lated and the other remains quasi-constant. While lowly
transcribing alleles are tuned to higher expression lev-
els predominantly by decreasing *T*<sub>OFF</sub>, medium-to-high
transcribing alleles are mostly tuned by increasing *T*<sub>ON</sub>.
These simple rules contain all information necessary to
govern bursting and consequently transcriptional activity
in the system.

**The common bursting relationships predict the effects of cis- and trans-perturbations.** Diverse regulatory mechanisms have been implicated in the control of transcriptional activity, including cis-regulatory elements (e.g., enhancers) and trans-factors (e.g., TF repressors). It is often assumed that distinct regulatory mechanisms directly control distinct bursting parameters. Will the established bursting parameter relationships based on wild-type measurements predict bursting dynamics when we perturb regulatory mechanisms?

To address this question, we devised a strategy to per-
turb the endogenous system in *cis* and in *trans*. For
*cis*-perturbations, we delete a distal *hb* enhancer from
the *hb* locus that has an MS2-stem loop cassette in the
first intron (Fig. 6A). This enhancer removal has a com-
plex effect on *hb* activity, including increased and de-
creased transcriptional activity at different times and AP
locations (Fig. 6B), consistent with previous observa-
tions [42–44]. Despite the stark deviation in the mutant
spatiotemporal transcriptional activity compared to the
wild-type, we find that mean transcription rates across
space and time are governed by *P*<sub>ON</sub>. Thus, the predic-
tive power of this parameter observed for the wild-type
holds for the mutant as well (Fig. 6B). The mutant fur-
ther adheres to the other bursting relationships identified
in the wild-type. In particular, the restrictive *T*<sub>ON</sub> and
*T*<sub>OFF</sub> to *P*<sub>ON</sub> relationships hold, as well as the largely
conserved switching correlation time *T*<sub>C</sub> (Fig. 6C).

A second enhancer deletion, namely, the removal of
the *kni* distal enhancer results in a significant reduction
in *kni* activity. The mutant samples have a smaller dy-
namic range of activity, yet we find a similar data collapse
within that range (Fig. S10A-C). Finally, to examine a
*trans*-perturbation, we measure *kni* activity in embryos
with a *hb* null background. *kni* activity was substan-
tially altered, consistent with earlier studies [45] (Fig.
S10D-E). Yet, again, we observe the collapse of the mu-
tant bursting parameters onto the wild-type busting rules
(Fig. 6D).

The consistency of these bursting rules suggests that
the wild-type derived relationships (Fig. 5F) can predict
how changes in T<sub>ON</sub> and T<sub>OFF</sub> account for the change in
transcriptional activity upon the perturbation. Specifi-
cally, to examine the perturbation's effect on transcrip-
tional activity at a given AP bin and time in devel-
opment, we can consider the pair {P<sub>ON</sub><sup>wt</sup>, P<sub>ON</sub><sup>mut</sup>} at that
spatiotemporal position (two examples of such pairs are
marked in the kymographs in Fig. 6B). Using the depen-
dencies of T<sub>ON</sub> and T<sub>OFF</sub> on P<sub>ON</sub>, we predict whether the
change in transcriptional activity (i.e., the change from
P<sub>ON</sub><sup>wt</sup> to P<sub>ON</sub><sup>mut</sup>) stems predominantly from a change in T<sub>ON</sub>
or in T<sub>OFF</sub> (Fig. 6E). For our two example pairs, one is
predicted to be mostly governed by T<sub>OFF</sub> ("o" mark),
while the other is mostly governed by T<sub>ON</sub> ("★" mark).
This exercise can be generalized to all {P<sub>ON</sub><sup>wt</sup>, P<sub>ON</sub><sup>mut</sup>} pairs
and subsequently verified using the measured T<sub>ON</sub> and
T<sub>OFF</sub> values for wild-type and the perturbations (Fig.
6C-D and Fig. S11A-C). We find that over ~ 90% of the
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

FIG. 6. Modulation of ON and OFF times by *cis*- and *trans*-perturbations. (A) Distal *hb* enhancer removal. The MS2-stem loops are inserted at the same location in the mutant (enhancer deletion) and wild-type fly lines. (B) Quantification of *hb* wild-type and mutant (A) phenotypes. Both transcriptional rate *R* (left) and *P*<sub>ON</sub> level (right, kymograph) display significantly different expression patterns in the mutant. Dotted arrow indicates time point in the kymograph at which rate profiles are depicted. “o” and “*” mark signal two bins with predominant *T*<sub>OFF</sub> modulation and predominant *T*<sub>ON</sub> modulation, respectively. (C-D) Transcriptional parameters for *hb* in *cis*-mutation (C, cyan) and for *kni* in *trans*-mutation (D, light green) collapse on corresponding wild-type parameters (gray). *cis*-mutation is the *hb* distal enhancer removal from A and the *trans*-mutation is a *hb* null background compared to wild-type (Fig. S10D-E). Solid black lines correspond to the endogenous bursting rules from Fig. 5F. (E) Example for *T*<sub>OFF</sub> and *T*<sub>ON</sub> modulation observed upon *hb* enhancer removal in B (circle and star, respectively). Depending on the pair of wild-type and mutant *P*<sub>ON</sub> chosen at the same spatiotemporal location (two examples: circle and star in B), the derived endogenous rules (gray and blue solid lines) predict either larger fold-change in *T*<sub>OFF</sub> compared to *T*<sub>ON</sub> (left, *T*<sub>OFF</sub> dominated modulation) or larger fold-change in *T*<sub>ON</sub> (right, *T*<sub>ON</sub> dominated modulation). Verifying these examples using the estimated *T*<sub>OFF</sub> and *T*<sub>ON</sub> (in C) confirms that both types of modulation can be observed for the same perturbation. (F) Verification of predicted *T*<sub>OFF</sub> and *T*<sub>ON</sub> modulation for all *hb* wild-type and mutant *P*<sub>ON</sub> pairs (defined as in E), for most pairs (> 85%) the prediction is correct (see Fig. S11A-C). (G) Verification of predicted burst size *B* and frequency *F* modulation for all *hb* wt and mutant *P*<sub>ON</sub> pairs. As in E, *B* versus *F* modulation can be predicted for each *P*<sub>ON</sub> pair based on the derived endogenous rules. For most pairs (> 95%) the prediction is correct (see Fig. S11G-I).
pairs correctly verify the type of modulation predicted by our rules (Fig. 6F). We reached the same conclusions predicting changes in burst size and burst frequency (Fig. 6G and S11G-I). Importantly we find that each type of perturbation displays both predominant $T_{ON}$ or $T_{OFF}$ modulation at different times and positions.

The generalization of our bursting rules to cis- and trans-perturbations have strong consequences. From our examination, the type of perturbed regulatory mechanism can hardly be linked to changes in a specific bursting parameter and vice versa. Indeed, as is the case for wild-type, $P_{ON}$ emerges as the main governing parameter of the transcriptional activity also for the mutants, while $K$ and $T_C$ remain largely unchanged. $P_{ON}$ changes upon a perturbation (i.e., wild-type to mutant) are sufficient to determine the corresponding changes in $R$, $T_{ON}$, and $T_{OFF}$. Moreover, the functional form of the relationships implies that the $P_{ON}$ regime (low versus high) is crucial in determining which parameter is predominantly affected. Different regulatory mechanisms have a different propensity for $P_{ON}$ modulation. Yet, the single-parameter regulation set in place by the identified bursting rules points towards general mechanisms conserving $K$ and $T_C$, and linking $T_{ON}$ and $T_{OFF}$.

DISCUSSION

While it is appreciated that mRNA production likely occurs in bursts across various systems, quantitatively measuring endogenous, single allele transcriptional bursting in real-time, across a wide range of gene activities, still poses a challenge. This hinders our understanding of how the kinetic parameters of bursting underlie transcriptional dynamics across genes, space, and developmental time. In this study, we devised an approach to perform such measurements in the context of the developing early *Drosophila* embryo, a system that relies on large changes in transcription rates, as a means to govern protein abundance (Fig. S4A). The spatiotemporal transcriptional activity of the examined genes is regulated by a myriad of processes (e.g. repressor and activator binding, chromatin accessibly, PIC formation and pause-release, histone modifications, etc.) [42, 46–48]. These processes are mediated by gene-specific cis-architectures, with distinct combinations of enhancers and promoters, further differing in their internal motif compositions. Surprisingly, despite the complexity of the regulatory processes involved and the differences between the genes, we find highly restricted, unifying properties of the underlying bursting dynamics.

We observed that the mean transcription rate is governed by a single tunable control parameter, the ON-probability ($P_{ON}$), with a near-constant Pol II initiation rate ($K$). We further found a conserved time scale ($T_C$) of about a minute, over which the allele states, either active (ON) or inactive (OFF), remain correlated. This time constant explicitly links the mean durations of ON

and OFF periods that are thus largely determined by the control parameter ($P_{ON}$). $P_{ON}$, together with the mostly conserved $K$ and $T_C$, fully parameterize the observed bursting dynamics, across genes, space, and developmental time.

While consistent with predictions from our previous measurements on fixed tissues (Fig. S12A) [25], the presented live measurements allow us to go beyond a model-based inference of kinetic parameters from a static snapshot of distributions of nascent transcripts. Our current approach measures the dynamics directly rather than inferring the kinetics indirectly. Therefore, we are relieved from the constraints imposed by a specific mechanistic model and can thus further relax the steady state assumption. Both aspects are often used in the analyses of fixed and live data [3, 6, 25, 49]. Moreover, access to the full dynamics of individual transcription time series (throughout more than 1 hour of development) allows us to estimate time-dependent transcriptional parameters from the underlying bursts, rather than obtaining only a single transcriptional parameter value per condition.

Using synthetic data, we verify the capacity of this procedure to reliably recover a wide range of bursting parameters, including the estimation of $T_C$ (Fig. S5C-D). The conserved nature of $T_C$ and its value (i.e., $1.25 \pm 0.35$ min as estimated from individual bursts) quantitatively agree with the auto-correlation analysis ($\tau_{AC} = 1.37 \pm 0.31$ min), an orthogonal approach not involving burst calling. It is intriguing to consider the functional consequences of a constant correlation time $T_C$ and of the relatively small measured value. $T_C$ not only sets the time scale of the bursting dynamics, linking the mean ON and OFF times, but it has further implications on transcription noise filtering: a small $T_C$ value minimizes noise as bursts are easily buffered by long mRNA lifetimes. A small $T_C$ further allows gene transcription to respond more rapidly to input TF changes, by means of facilitating a fast relaxation to a steady state (Fig. S7J-K).

The identified bursting rules point to a surprising predictive power of $P_{ON}$ on the mean ON and OFF times (Fig. 5F). This observation is only possible because our system allows for the quantification of bursting parameters across a large dynamic range of transcriptional activity, with $P_{ON}$ values ranging from fully off (0) to fully on (1) for most of our examined genes and conditions. This leads us to uncover a strict relationship between $T_{OFF}$ and $T_{ON}$ (Fig. 6E and Fig. 7A). It is commonly thought that distinct regulatory processes will alter transcriptional activity by predominant regulation of specific bursting parameters (e.g., the mean ON durations of a burst versus the mean OFF intervals in between bursts). Yet, when we perform a cis- or trans-perturbation, which substantially alters transcriptional activity, the predominantly modulated bursting parameters ($T_{OFF}$ versus $T_{ON}$ or burst frequency versus burst size) can be predicted by the wild-type and mutant $P_{ON}$ regimes, rather than by the type of the perturbation performed (Fig. 6E-G and
<figure><img src="image_15.png" /><figcaption>FIG. 7. Decoupling between mechanism and response points to unifying rules. (A) Scatter plot of the transcriptional parameters as a function of $P_{ON}$ (color code same as Fig. 5F). Transcriptional parameters from two other *Drosophila* studies are largely consistent with uncovered transcription rules (black lines); one study investigates response to BMP signaling (yellow circles; Hoppe et al.), the other the effect of core promoter manipulations on a transgene (pink circles; Pimmet et al.). Transcriptional parameters resulting from multiple perturbations performed on yeast *GAL10* gene all closely follow the rules (orange circles; Brouwer et al.). Besides the initiation rate *K*, transcriptional parameters inferred from mammalian promoters (cyan circles; Suter et al.) appear mostly compatible with the rules, suggesting that these may apply beyond *Drosophila*. A minimal description of bursts requires three independent parameters, such as initiation rate *K*, switching correlation time $T_C$, and ON-probability $P_{ON}$. $P_{ON}$ is the main regulated parameter controlling the transcription rate *R*, while *K* only accounts for small changes (at most 2-fold), and the correlation time $T_C$ sets a conserved bursting time scale. Conservation of both *K* and $T_C$ implies that $T_{OFF}$ and $T_{ON}$, or, alternatively, the burst size *B* and burst frequency *F*, are fully determined by $P_{ON}$. Thus, whether changes in *R* are mediated by $T_{OFF}$ versus $T_{ON}$, or *B* versus *F* dominated modulation only depends on $P_{ON}$. (B) While *K* and $T_C$ appear mostly conserved, $P_{ON}$ emerges as the key regulatory parameter. It is independent of gene identity (different *cis*-architecture) or the combination of input transcription factor concentrations (that vary along the AP axis), validated by *cis*- and *trans*-perturbations and by data from other studies. Such invariance suggests a mechanistic decoupling, where regulation affects solely $P_{ON}$, which in turn determines unequivocally all of the underlying bursting dynamics.</figcaption></figure>

S11).

To further examine the generality of these results, we investigated previously published transcription measurements in the early fly embryo. These measurements include both endogenous genes and synthetic reporters, where transcription was altered by varying BMP signal (a dorsoventral morphogen) [28] or core promoter composition [20]. Strikingly, we find that these datasets collapse on our established bursting rules (Fig. 7A). As suggested in these studies, the first dataset shows mainly effects on $T_{OFF}$, while the latter principally changes $T_{ON}$. Intriguingly, the two independent datasets cluster in disjoint halves of the full spectrum of $P_{ON}$ values captured by our measurements. Our analysis thus raises the possi-

bility that the predominantly changed parameter ($T_{OFF}$ versus $T_{ON}$) might not be inherent to the examined regulatory manipulation (e.g., input TF concentrations or core promoter elements), but rather a consequence of the limited expression range of these genes.

The striking manifestation of general rules underlying bursting parameters in the *Drosophila* embryo, as well as the conserved nature of regulatory processes and the transcription machinery across eukaryotes, naturally leads to the question of whether our rules apply more broadly [50]. Recent measurements of an extensively perturbed yeast gene [21] provide an equivalent set of bursting parameters that faithfully adhere to our transcription rules (Fig. 7A). Although, the yeast gene initiation rate
$K$ appears mostly constant across conditions (similar to the fly genes), calibration of its value was only possible under specific hypotheses leading to a lower bound that needs to be further tested (Methods). Earlier yeast gene data have shown a significantly smaller initiation rate [5, 38]. A comprehensive study allowing multiple gene classes across multiple organisms will be necessary to verify whether our rules hold more generally in eukaryotes.

Mammalian genes are often lowly transcribed [6, 29, 51], potentially exploring a different parameter regime. Indeed, while bursting parameters inferred from measurements of a luminescent reporter protein [6] are restricted to very small $P_{ON}$ values (0 to 0.05), they appear to be consistent with our established relationships (Fig. 7A). Again, only the initiation rate $K$ seems to deviate. Possible deviations in $K$ across species leave open the possibility that the overall levels of $K$ might be linked to species-specific rates, such as those linked to metabolism [52, 53]. Mammalian time-lapse data with a broader $P_{ON}$ range will be necessary to make a stronger parameter comparison possible.

Revisiting previously performed genome-wide studies in mammalian systems shows trends that are possibly compatible with our established bursting relationships. Similar to our fly genes (Fig. 7A), one study found that while burst frequency is predominantly modulated for low expressing genes, burst size is tuned for high expressing ones, independent on the reporter control sequences [13]. Another study using single-cell RNA-seq found a functional dependence of burst frequency and burst size on mean expression that seems compatible with our established rules (for $P_{ON} < 0.5$) [51]. However, scRNA-seq parameters are significantly influenced by mRNA loss and long mRNA lifetimes, making the mapping of the sets of units between the vastly different approaches challenging.

The potentially wider applicability of our bursting relationships to other species calls for a new framework underlying the regulatory processes governing transcription. Instead of specific regulatory processes being inherently linked to specific bursting parameters, the tuning of transcription seems to be funneled through the sole control

of $P_{ON}$, that all regulatory processes act on (Fig. 7B). Future investigations will have to determine the molecular mechanisms that can implement such a funneling at the level of $P_{ON}$ control. How do diverse processes tune $P_{ON}$? How is the constancy of the switching correlation time $T_C$ implemented molecularly? As for the latter, our work suggests that a highly conserved and general mechanism across eukaryotic genes should be at work. The generality of such a mechanism implies independence from the particularity of a given gene locus and thus could rather be implemented by the environment, including structural consideration of the nuclear architecture [54–56] or the molecular assembly of components of the transcription machinery [57, 58].

Despite the complexity of transcriptional dynamics across species, genes, space, developmental time, and perturbations, our quantitative real-time measurements revealed strict bursting rules, that set strong constraints on mechanistic models of transcriptional regulation. Our work also has some indicators for the generality of these rules across systems, their functional implications, and their molecular underpinning. As is the case with other areas in which organizing principles are increasingly emerging [59–63], these rules offer new ways to think about complex processes and point to conserved mechanisms at their core.

## ACKNOWLEDGMENTS

We thank members of the Gregor laboratory for discussions and comments on the manuscript; E. Wieschaus for suggestions at various stages of the project. This work was supported in part by the U.S. National Science Foundation, through the Center for the Physics of Biological Function (PHY-1734030), and by National Institutes of Health Grants R01GM097275, U01DA047730, and U01DK127429. M.L. is the recipient of a Human Frontier Science Program fellowship (LT000852/2016-L), EMBO long-term postdoctoral fellowship (ALTF 1401-2015), and the Rothschild postdoctoral fellowship.

[1] K. M. Lelli, M. Slattery, and R. S. Mann, Disentangling the many layers of eukaryotic transcriptional regulation., Annual Review of Genetics **46**, 43 (2012).

[2] P. Cramer, Eukaryotic transcription turns 50, Cell **179**, 808 (2019).

[3] A. Raj, C. S. Peskin, D. Tranchina, D. Y. Vargas, and S. Tyagi, Stochastic mRNA synthesis in mammalian cells., PLoS biology **4**, e309 (2006).

[4] J. R. Chubb, T. Trcek, S. M. Shenoy, and R. H. Singer, Transcriptional Pulsing of a Developmental Gene, Current Biology **16**, 1018 (2006).

[5] D. Zenklusen, D. R. Larson, and R. H. Singer, Single-

RNA counting reveals alternative modes of gene expression in yeast., Nature Structural & Molecular Biology **15**, 1263 (2008).

[6] D. M. Suter, N. Molina, D. Gatfield, K. Schneider, U. Schibler, and F. Naef, Mammalian genes are transcribed with widely different bursting kinetics., Science **332**, 472 (2011).

[7] J. P. Bothma, H. G. Garcia, E. Esposito, G. Schlissel, T. Gregor, and M. Levine, Dynamic regulation of eve stripe 2 expression reveals transcriptional bursts in living drosophila embryos., Proc. Natl. Acad. Sci. USA **111**, 10598 (2014).
[8] K. Tantale, F. Mueller, A. Kozulic-Pirher, A. Lesne, J.-M. Victor, M.-C. Robert, S. Capozi, R. Chouaib, V. Bäcker, J. Mateos-Langerak, X. Darzacq, C. Zimmer, E. Basyuk, and E. Bertrand, A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting, Nature Communications 7, 12248 (2016).

[9] D. Nicolas, N. E. Phillips, and F. Naef, What shapes eukaryotic transcriptional bursting?, Molecular BioSystems 13, 1280 (2017).

[10] E. Tunnacliffe and J. R. Chubb, What Is a Transcriptional Burst?, Trends in Genetics 36, 288 (2020).

[11] J. Rodriguez and D. R. Larson, Transcription in Living Cells: Molecular Mechanisms of Bursting, Annual Review of Biochemistry 89, 1 (2020).

[12] J. Peccoud and B. Ycart, Markovian modeling of gene-product synthesis, Theoretical Population Biology 48, 222 (1995).

[13] R. D. Dar, B. S. Razooky, A. Singh, T. V. Trimeloni, J. M. McCollum, C. D. Cox, M. L. Simpson, and L. S. Weinberger, Transcriptional burst frequency and burst size are equally modulated across the human genome., Proc. Natl. Acad. Sci. USA 109, 17454 (2012).

[14] A. Senecal, B. Munsky, F. Proux, N. Ly, F. Braye, C. Zimmer, F. Mueller, and X. Darzacq, Transcription Factors Modulate c-Fos Transcriptional Bursts, Cell Reports 8, 75 (2014).

[15] C. R. Bartman, S. C. Hsu, C. C.-S. Hsiung, A. Raj, and G. A. Blobel, Enhancer Regulation of Transcriptional Bursting Parameters Revealed by Forced Chromatin Looping., Molecular Cell 62, 237 (2016).

[16] T. Fukaya, B. Lim, and M. Levine, Enhancer control of transcriptional bursting., Cell 166, 358 (2016).

[17] C. Li, F. Cesbron, M. Oehler, M. Brunner, and T. Höfer, Frequency Modulation of Transcriptional Bursting Enables Sensitive and Rapid Gene Regulation, Cell Systems 6, 409 (2018).

[18] B. T. Donovan, A. Huynh, D. A. Ball, H. P. Patel, M. G. Poirier, D. R. Larson, M. L. Ferguson, and T. L. Lenstra, Live-cell imaging reveals the interplay between transcription factors, nucleosomes, and bursting., The EMBO Journal 38, 10.15252/embj.2018100809 (2019).

[19] D. A. Stavreva, D. A. Garcia, G. Fettweis, P. R. Gudla, G. F. Zaki, V. Soni, A. McGowan, G. Williams, A. Huynh, M. Palangat, R. L. Schiltz, T. A. Johnson, D. M. Presman, M. L. Ferguson, G. Pegoraro, A. Upadhyaya, and G. L. Hager, Transcriptional Bursting and Co-bursting Regulation by Steroid Hormone Release Pattern and Transcription Factor Mobility, Molecular Cell 75, 1161 (2019).

[20] V. L. Pimmett, M. Dejean, C. Fernandez, A. Trullo, E. Bertrand, O. Radulescu, and M. Lagha, Quantitative imaging of transcription in living Drosophila embryos reveals the impact of core promoter motifs on promoter state dynamics, Nature Communications 12, 4504 (2021).

[21] I. Brouwer, E. Kerklingh, F. van Leeuwen, and T. L. Lenstra, Dynamic epistasis analysis reveals how chromatin remodeling regulates transcriptional bursting, bioRxiv, 2021.12.15.472793 (2021).

[22] V. L. Bass, V. C. Wong, M. E. Bullock, S. Gaudet, and K. Miller-Jensen, TNF stimulation primarily modulates transcriptional burst size of NF-κB-regulated genes, Molecular Systems Biology 17, e10127 (2021).

[23] G. Neuert, B. Munsky, R. Z. Tan, L. Teytelman, M. Khammash, and A. v. Oudenaarden, Systematic Identification of Signal-Activated Stochastic Gene Regulation, Science 339, 584 (2013).

[24] H. Xu, L. A. Sepúlveda, L. Figard, A. M. Sokac, and I. Golding, Combining protein and mRNA quantification to decipher transcriptional regulation, Nature Methods 12, 739 (2015).

[25] B. Zoller, S. C. Little, and T. Gregor, Diverse Spatial Expression Patterns Emerge from Unified Kinetics of Transcriptional Bursting, Cell 175, 835 (2018).

[26] S. Shah, Y. Takei, W. Zhou, E. Lubeck, J. Yun, C.-H. L. Eng, N. Koulena, C. Cronin, C. Karp, E. J. Liaw, M. Amin, and L. Cai, Dynamics and Spatial Genomics of the Nascent Transcriptome by Intron seqFISH, Cell 174, 363 (2018).

[27] K. Tantale, E. Garcia-Oliver, M.-C. Robert, A. L'Hostis, Y. Yang, N. Tsanov, R. Topno, T. Gostan, A. Kozulic-Pirher, M. Basu-Shrivastava, K. Mukherjee, V. Slaninova, J.-C. Andrau, F. Mueller, E. Basyuk, O. Radulescu, and E. Bertrand, Stochastic pausing at latent HIV-1 promoters generates transcriptional bursting, Nature Communications 12, 4503 (2021).

[28] C. Hoppe, J. R. Bowles, T. G. Minchington, C. Sutcliffe, P. Upadhyai, M. Rattray, and H. L. Ashe, Modulation of the Promoter Activation Rate Dictates the Transcriptional Response to Graded BMP Signaling Levels in the Drosophila Embryo, Developmental Cell 54, 727 (2020).

[29] Y. Wan, D. G. Anastasakis, J. Rodriguez, M. Palangat, P. Gudla, G. Zaki, M. Tandon, G. Pegoraro, C. C. Chow, M. Hafner, and D. R. Larson, Dynamic imaging of nascent rna reveals general principles of transcription dynamics and stochastic splice site selection., Cell 184, 2878 (2021).

[30] X. Darzacq, Y. Shav-Tal, V. d. Turris, Y. Brody, S. M. Shenoy, R. D. Phair, and R. H. Singer, In vivo dynamics of RNA polymerase II transcription, Nature Structural & Molecular Biology 14, 796 (2007-09).

[31] H. G. Garcia, M. Tikhonov, A. Lin, and T. Gregor, Quantitative imaging of transcription in living Drosophila embryos links polymerase activity to patterning., Current Biology 23, 2140 (2013).

[32] J. Li, A. Dong, K. Saydaminova, H. Chang, G. Wang, H. Ochiai, T. Yamamoto, and A. Pertsinidis, Single-Molecule Nanoscopy Elucidates RNA Polymerase II Transcription at Single Genes in Live Cells, Cell 178, 491 (2019).

[33] T. Gregor, H. G. Garcia, and S. C. Little, The embryo as a laboratory: quantifying transcription in Drosophila., Trends in Genetics 30, 364 (2014).

[34] S. Little, M. Tikhonov, and T. Gregor, Precise Developmental Gene Expression Arises from Globally Stochastic Transcriptional Activity, Cell 154, 789 (2013).

[35] T. Lucas, T. Ferraro, B. Roelens, J. D. L. H. Chanes, A. M. Walczak, M. Coppey, and N. Dostatni, Live imaging of bicoid-dependent transcription in drosophila embryos., Current Biology 23, 2135 (2013).

[36] M. Levo, J. Raimundo, X. Y. Bing, Z. Sisco, P. J. Batut, S. Ryabichko, T. Gregor, and M. S. Levine, Transcriptional coupling of distant regulatory genes in living embryos., Nature 605, 754 (2022).

[37] E. Bertrand, P. Chartrand, M. Schaefer, S. M. Shenoy, R. H. Singer, and R. M. Long, Localization of ASH1 mRNA particles in living yeast., Molecular Cell 2, 437 (2018).
(1998).

[38] D. R. Larson, D. Zenklusen, B. Wu, J. A. Chao, and R. H. Singer, Real-time observation of transcription initiation and elongation on an endogenous yeast gene., *Science* **332**, 475 (2011).

[39] J. Liu, D. Hansen, E. Eck, Y. J. Kim, M. Turner, S. Alamos, and H. G. Garcia, Real-time single-cell characterization of the eukaryotic transcription cycle reveals correlations between RNA initiation, elongation, and cleavage, *PLoS Computational Biology* **17**, e1008999 (2021).

[40] J. O. Dubuis, R. Samanta, and T. Gregor, Accurate measurements of dynamics and reproducibility in small genetic networks., *Molecular Systems Biology* **9**, 639 (2013).

[41] S. L. McKnight and O. L. Miller, Post-replicative nonribosomal transcription units in D. Melanogaster embryos., *Cell* **17**, 551 (1979).

[42] M. W. Perry, A. N. Boettiger, and M. Levine, Multiple enhancers ensure precision of gap gene-expression patterns in the Drosophila embryo., *Proc. Natl. Acad. Sci. USA* **108**, 13570 (2011).

[43] J. P. Bothma, H. G. Garcia, S. Ng, M. W. Perry, T. Gregor, and M. Levine, Enhancer additivity and non-additivity are determined by enhancer strength in the Drosophila embryo., *eLife* **4**, 10.7554/eLife.07956 (2015).

[44] T. Fukaya, Dynamic regulation of anterior-posterior patterning genes in living Drosophila embryos., *Current Biology* **31**, 2227 (2021).

[45] M. W. Perry, J. P. Bothma, R. D. Luu, and M. Levine, Precision of hunchback expression in the Drosophila embryo., *Current Biology* **22**, 2247 (2012).

[46] G. Farkas, B. A. Leibovitch, and S. C. Elgin, Chromatin organization and transcriptional control of gene expression in Drosophila, *Gene* **253**, 117 (2000).

[47] M. Lagha, J. P. Bothma, and M. Levine, Mechanisms of transcriptional precision in animal development, *Trends in Genetics* **28**, 409 (2012).

[48] J. Park, J. Estrada, G. Johnson, B. J. Vincent, C. Ricci-Tam, M. D. Dragdon, Y. Shulgina, A. Cha, Z. Wunderlich, J. Gunawardena, and A. H. DePace, Dissecting the sharp response of a canonical developmental enhancer reveals multiple sources of cooperativity, *eLife* **8**, e41266 (2019).

[49] J. Liu, D. Hansen, E. Eck, Y. J. Kim, M. Turner, S. Alamos, and H. G. Garcia, Real-time single-cell characterization of the eukaryotic transcription cycle reveals correlations between RNA initiation, elongation, and cleavage, *PLoS Computational Biology* **17**, e1008999 (2021).

[50] A. Sanchez and I. Golding, Genetic Determinants and Cellular Constraints in Noisy Gene Expression, *Science* **342**, 1188 (2013).

[51] A. J. M. Larsson, P. Johnsson, M. Hagemann-Jensen, L. Hartmanis, O. R. Faridani, B. Reinius, r. Segerstolpe, C. M. Rivera, B. Ren, and R. Sandberg, Genomic encoding of transcriptional burst kinetics, *Nature* **565**, 251 (2018).

[52] T. Rayon, D. Stamataki, R. Perez-Carrasco, L. Garcia-Perez, C. Barrington, M. Melchionda, K. Exelby, J. Lazaro, V. L. J. Tybulewicz, E. M. C. Fisher, and J. Briscoe, Species-specific pace of development is associated with differences in protein stability, *Science* **369**, 10.1126/science.aba7667 (2020).

[53] M. Diaz-Cuadros, T. P. Miettinen, O. S. Skinner, D. Sheedy, C. M. Díaz-García, S. Gapon, A. Hubaud, G. Yellen, S. R. Manalis, W. M. Oldham, and O. Pourquié, Metabolic regulation of species-specific developmental rates, *Nature* **613**, 550 (2023).

[54] A. Tsai, A. K. Muthusamy, M. R. Alves, L. D. Lavis, R. H. Singer, D. L. Stern, and J. Crocker, Nuclear microenvironments modulate transcription from low-affinity enhancers, *eLife* **6**, e28975 (2017).

[55] J. E. Henninger, O. Oksuz, K. Shrinivas, I. Sagi, G. LeRoy, M. M. Zheng, J. O. Andrews, A. V. Zamudio, C. Lazaris, N. M. Hannett, T. I. Lee, P. A. Sharp, I. I. Cissé, A. K. Chakraborty, and R. A. Young, RNA-Mediated Feedback Control of Transcriptional Condensates, *Cell* **184**, 207 (2021).

[56] D. B. Brückner, H. Chen, L. Barinov, B. Zoller, and T. Gregor, Stochastic motion and transcriptional dynamics of pairs of distal dna loci on a compacted chromosome, *bioRxiv*, 2023.01.18.524527 (2023).

[57] W.-K. Cho, J.-H. Spille, M. Hecht, C. Lee, C. Li, V. Grube, and I. I. Cisse, Mediator and RNA polymerase II clusters associate in transcription-dependent condensates., *Science* **361**, 412 (2018).

[58] V. Q. Nguyen, A. Ranjan, S. Liu, X. Tang, Y. H. Ling, J. Wisniewski, G. Mizuguchi, K. Y. Li, V. Jou, Q. Zheng, L. D. Lavis, T. Lionnet, and C. Wu, Spatiotemporal coordination of transcription preinitiation complex assembly in live cells, *Molecular Cell* **81**, 3560 (2021).

[59] G. Tkačik, C. G. Callan, and W. Bialek, Information flow and optimization in transcriptional regulation, *Proc. Natl. Acad. Sci. USA* **105**, 12265 (2008), 0705.0313.

[60] D. L. Jones, R. C. Brewster, and R. Phillips, Promoter architecture dictates cell-to-cell variability in gene expression, *Science* **346**, 1533 (2014).

[61] J. Hausser, A. Mayo, L. Keren, and U. Alon, Central dogma rates and the trade-off between precision and economy in gene expression, *Nature Communications* **10**, 68 (2018).

[62] M. D. Petkova, G. Tkačik, W. Bialek, E. F. Wieschaus, and T. Gregor, Optimal Decoding of Cellular Identities in a Genetic Network, *Cell* **176**, 844 (2019).

[63] R. Balakrishnan, M. Mori, I. Segota, Z. Zhang, R. Aebersold, C. Ludwig, and T. Hwa, Principles of gene regulation quantitatively connect DNA to RNA and proteins in bacteria, *Science* **378**, eabk2066 (2022).

[64] K. Becker, E. Balsa-Canto, D. Cicin-Sain, A. Hobermann, H. Janssens, J. R. Banga, and J. Jaeger, Reverse-Engineering Post-Transcriptional Regulation of Gap Genes in Drosophila melanogaster, *PLoS Computational Biology* **9**, e1003281 (2013).
SUPPLEMENTAL FIGURES

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>FIG. S1. Signal calibration, measurement error and embryo-embryo variability. [Caption see next page.]</figcaption></figure>
FIG. S1. **Signal calibration, measurement error and embryo-embryo variability.** (A) The four trunk gap genes, giant (gt), hunchback (hb), Kruppel (Kr) and knirps (kni) were imaged using the MS2/PP7 stem-loop labeling systems. Stem-loop cassettes (vertical black arrow) were inserted either in the first/second intron or in the 3'UTR of each gene. Gap genes harbor different cis-architectures as characterized by the number of promoters, number of enhancers (color boxes), and composition of these cis-regulatory elements (TF binding motifs, core promoter elements, etc.). (B) Relative calibration unit (left) and relative error (right) for each gap gene construct related to Fig. 1C. The conversion of the live signal to absolute units is performed by comparison to a smFISH-based measurement. Calibration was performed by matching the mean full embryo length transcriptional activity profiles (reconstructed by averaging over all nuclei in 2.5% AP bins, within a 5 min time window in NC13) measured by live imaging to those previously measured with smFISH (Zoller et al., 2018). (Left) The procedure was performed using all measured gap profiles at once, leading to a final calibration unit (horizontal back line, dashed lines are plus/minus one standard error). We then repeated the procedure for each individual construct separately (color circles), and the derived units are expressed in percent of the global fit. (Right) Relative error for the calibration unit of each individual gene construct with respect to the global unit and mean relative error (dashed line), which is below 5%. Error bars are 68% confidence intervals. (C) Comparison of higher cumulants versus mean activity relationships obtained by live imaging and smFISH measurements; right column panels are reproduced from Figure 3B-D in Zoller et al. (2018), left column panels are from live data analyzed equivalently. Live cumulants of transcriptional activity (mean, variance, 3rd and 4th cumulant) are estimated over all nuclei in 2.5% AP bins, within a 5 min time window in NC13. Cumulants are converted from equivalent cytoplasmic mRNA units (C.U.) to Pol II counts for a single gene copy of average length (3.3 kb). The cumulants are normalized with respect to $g_0$ defined as the intercept of the Poisson background (dashed line) and the polynomial fit to the data (black solid line for live and doted for smFISH). The number $g_0$ can be interpreted as the mean number of Pol II on a 3.3 kb long gap gene at maximal activity. We get $g_0 = 13.6$ for live and $g_0 = 15.2$ with smFISH measurements, a difference of 12%. Overall, the higher cumulants versus mean relationships obtained from live (left column) and from smFISH (right column) are extremely close (black solid versus dotted line), confirming the quantitative nature and the proper calibration of our live assay. Two independent methods leading to the same quantitative conclusions validate each other reciprocally. It strongly suggests that our synthetic modifications of the endogenous gap gene loci have no currently measurable effect on the transcriptional output of the system. (D) Two independent methods to assess the imaging error. (Left) An interlaced cassette of alternating MS2 and PP7 stem-loops, labeled with two differently colored coat proteins (MCP-GFP and PCP-mCherry), is inserted in the first intron of *Kr*. In absence of imaging error, the transcriptional activity in the green and red channels when calibrated to C.U. should perfectly correlate (on the diagonal). We fitted the spread $\sigma_{img}$ orthogonal to the diagonal (black line, slope one) to characterize the imaging error; assuming $\sigma_{img}^2$ scales as $\sigma_b^2 + \alpha I$ with mean intensity $I$, where $\sigma_b^2$ is the background noise and $\alpha I$ a Poisson shot noise term. The resulting fit for $\sigma_{img}$ is highlighted by the dashed lines (plus minus one std around the diagonal). (Right) Imaging error estimation from the single allele transcriptional time series (with the assumption that the measured transcriptional fluctuations result from the sum of uncorrelated imaging noise and correlated noise due to the elongation of tagged nascent transcripts). We computed the time-dependent mean activity $\mu(t)$, variance $\sigma^2(t)$ and covariance between consecutive time points $Cov(t, t + \Delta t)$ (where $\Delta t$ is 10 s), over all nuclei within 1.5-2.5% AP bins for all measured genes. The uncorrelated imaging variability $\sigma_{img}^2$ is then approximated by $\sigma_{img}^2 - Cov(t, t + \Delta t)$, which is plotted as a function of $\mu(t)$ for all time points. We characterized $\sigma_{img}^2$ by fitting the data with a line $\sigma_b^2 + \alpha\mu$. Fitting results are shown in E. (E) Our two estimates for the imaging error (interlaced dual-color construct (solid line) and correlation-based approach (dashed line)) are consistent. The signal-to-noise-ratio (SNR), defined as $\mu/\sigma_{img}$, is close to one (dotted line) when $\mu \approx 1$, indicating that the sensitivity of our live measurements is close to one mRNA molecule. (F) Fractional embryo variability profiles as a function of AP position and developmental time, for all gap genes. We define embryo variability $\sigma_{emb}^2$ as the variance of the mean activity across embryos, and report the fractional embryo variability as the ratio $\sigma_{emb}^2/\sigma^2$, where $\sigma^2 = \sigma_{emb}^2 + \sigma_{nuc}^2$ is the total variance, and $\sigma_{nuc}^2$ corresponds to the transcriptional allele-to-allele noise across nuclei. Overall, the fractional embryo variability is ~10%, meaning that most of the variability arises from $\sigma_{nuc}^2$. Thus, together D, E, and F show that $\sigma_{tot}^2 = \sigma_{img}^2 + \sigma_{emb}^2 + \sigma_{nuc}^2$ is a good proxy for $\sigma_{nuc}^2$, which is the relevant noise contribution that contains all the bursting phenomenology.
<figure><img src="image_24.png" /><figcaption>FIG. S2. Dual color measurements to validate single-cell deconvolution and measure elongation rate. (A) Validation of the kernel assumption for the deconvolution of initiation events from single allele transcription time series using a dual-color (confocal) imaging approach for *hb* and *Kr*. For *hb* (*Kr*), we generated fly lines with dual-insertions of a MS2 (PP7) stem-loop cassette in the respective first intron and a PP7 (MS2) stem-loop cassette in the 3'UTR. In both cases, the two cassettes were labelled using two different colors (MCP-GFP green and PCP-mCherry red). Since the two signals are correlated through the elongation process, the simultaneously measured pair of time series has a further constrained set of underlying initiation configurations and represents thus a good test for the approach. To deconvolve single allele dual color time series together (i.e., a single train of polymerases needs to match two signals), using two kernels modeling each loop-cassette location and satisfying our key assumptions (i. constant and deterministic elongation rate; ii. no Pol II pausing/dropping in gene body; iii. fast termination). In addition, the dual-color strategy allows estimation of the average elongation rate from the overall delay between the two signals (using the known genomic distance between the MS2 and PP7 insertion sites). (B) Dual-color signal reconstructed from deconvolved single allele transcription time series (black lines for raw measured data). single allele transcription rate (gray line with one std envelope) is deconvolved from the single depicted pair of measured time series (black lines). The signal (red and green lines with one std envelope) is devoid of imaging noise (as it was modeled from Fig. S1D during the deconvolution process) and is reconstructed by convolving back the resulting transcription rate with the kernel of each channel. Qualitatively, the signal (color) matches well (see C) the measured time series (black) in strong support of our kernel assumptions. (C) Distribution of residuals from the dual-color reconstruction. We quantified the mean and standard deviation of the normalized residuals, i.e., of the difference between the measured signal (black in H) and the reconstructed signal (color in B) divided by the standard deviation of the imaging noise, for each recorded individual allele (for *hb* *N* = 2666 (blue) and for *Kr* *N* = 2594 (pink)). Overall, the dispersion of the means and standard deviations of normalized residuals (black line, 95% confidence ellipse) is close to the expected dispersion of a perfect model (dotted line, 95% confidence ellipse). (D-E) Estimated elongation rate *K*<sub>elo</sub> from dual-color measurements. (D) Average elongation rate computed over nuclei across 10 embryos as a function AP position (both *hb* (blue) and *Kr* (pink)) in NC13 (square) and NC14 (circle), with error bars representing one standard deviation across the embryo means. (E) Average elongation rate computed for individual embryos (color code and symbols as D), with error bars representing the standard deviation across the means over positions. The elongation rate is globally conserved across genes and nuclear cycles, with *K*<sub>elo</sub> = 1.8 ± 0.1 kb/min (corresponding to the mean across embryos (black line) plus/minus one standard deviation (dashed line)).</figcaption></figure>
A

<figure><img src="image_25.png" /></figure>

B

<figure><img src="image_26.png" /></figure>

C

<figure><img src="image_27.png" /></figure>

D

<figure><img src="image_28.png" /></figure>

E

<figure><img src="image_29.png" /></figure>

F

<figure><img src="image_30.png" /></figure>

G

<figure><img src="image_31.png" /><figcaption>FIG. S3. Transcription rate fluctuations reveal key bursting characteristics. [Caption see next page.]</figcaption></figure>
FIG. S3. **Transcription rate fluctuations reveal key bursting characteristics.** (A) Snapshots of the gap gene mean transcription rate *R* as a function of AP position in late NC13 and at 40 min in NC14. Gap gene profiles (color) are obtained by averaging the deconvolved single allele transcription rate over all nuclei within each AP bin (width of 2.5% and 1.5% embryo egg length in NC13 and NC14 respectively) and at each time point (10 s temporal resolution). The black dashed lines correspond to the mean activity of each gap gene at the same position and time normalized by the effective elongation time. Both the colored and dashed profiles agree, justifying our deconvolution approach. Error bars are one standard deviation across embryo means. Overall, we have effectively deconvolved $N_g = 7$ “genes” (4+*gt* male and female and anterior and posterior regions), over $N_t = 362$ time points (NC13+NC14), across $N_x = 9–18$ positions, leading to a total of 33′214 bins, each averaging ~200 nuclei (with a single allele per nucleus). (B) Fraction of spatial and temporal bins whose single allele transcription rate distribution $P(r)$ is consistent with the conditional transcription rate distribution $P(r|R)$ determined by pooling nuclei over multiple bins at a given mean transcription rate $R$. We computed the 95% confidence interval on the cumulative distribution of $P(r|R)$ and checked for all the underlying bins at a given $R$ whether their individual cumulative distribution was within the overall confidence interval. We repeated this process for four distinct developmental time windows: NC13 ($6.5 \le t$ min after mitosis) plus early NC14 ($7.5 \le t < 20.5$ min), mid NC14 ($20.5 \le t < 34.5$ min), late NC14 ($34.5 \le t < 48$ min), and a wider NC14 window ($7.5 \le t < 48$ min). Overall, bins that share similar $R$ within the same time window have very similar $P(r)$ distribution (median given by dashed line over 80%), which justifies the pooling of these bins. On the other hand, when pooling bins over the whole NC14 we observe further dissimilarities between bins, suggesting that $P(r|R)$ might moderately change over time. (C) 3rd cumulant and 4th cumulant of single allele transcription rate as a function of mean transcription rate $R$ in NC13 (square) and early NC14 ($7.5 \le t < 20.5$ min; circle). The single allele transcription rates are estimated within 1-min-intervals, highlighting a strong departure of the variance from a constitutive Poisson regime (dashed line, $\kappa_3$ and $\kappa_4 = R$). All gap genes follow the same trend suggesting a common bursting regime. (D) Interpreting auto-correlation function using the 2-state model of transcriptional bursting. In this model for a single gene copy (top), the gene promoter switches stochastically between an OFF and ON state with rates $k_{OFF}$ and $k_{ON}$, where in the latter Pol II can be loaded at rate $K^{(1)}$ and elongate at rate $K_{elo}$. Of note, the superscript (1) specifies the parameter for a single gene copy. The auto-correlation functions displayed here are computed from the model using a switching correlation time $T_C^{(1)} = 1/(k_{OFF} + k_{ON}) = 2$ min, a Pol II elongation time $\tau_{elo} = L_g/K_{elo} = 2$ min (where $L_g$ is the gene length) and an initiation rate $K^{(1)} = 8$ mRNA/min; the steady state ON-probability $P_{ON}^{(1)} = k_{ON}/(k_{ON} + k_{OFF})$ varies from 0 to 1 (i.e., blue to red color code, fraction of nuclei in ON state or fraction of time a nucleus is in ON state). In principle, promoter switching (generating bursts) leads to temporal correlations in the transcriptional activity time series (Activity AC). However, from the raw live measurements, these correlations are hard to distinguish from the ones introduced by elongation (left), specifically when the switching correlation time $T_C^{(1)}$ is close to or smaller than the elongation time $\tau_{elo}$. Instead, performing auto-correlation analysis on deconvolved single allele transcription rates resolves the switching correlations (Transcription AC, right), since correlations due to elongation have been removed. Thus, the switching correlation time $T_C^{(1)}$ can be estimated by fitting an exponential to the decay of the Transcription AC. (E) Expected effect of the ON-probability (left) and Pol II initiation rate (right) on the magnitude of the correlated noise, $\Sigma_{AC}$. Computation of $\Sigma_{AC}$ using same parameters as in D. (Left) As $P_{ON}^{(1)}$ increases, the transcription rate $R = 2K^{(1)}P_{ON}^{(1)}$ increases (here for 2 sister chromatids), and $\Sigma_{AC}$ decreases until it vanishes at $P_{ON}^{(1)} = 1$ (Poisson regime). This behavior is consistent with our data, see Fig. 2E. (Right) At fixed $P_{ON}^{(1)}$ ($P_{ON}^{(1)} = 0.5$) and varying initiation rate $K^{(1)}$, the $\Sigma_{AC}$ increases with growing transcription rate $R$. This behavior is the opposite of what we observed in our data (Fig. 2E). In both cases (Left and Right) the dotted line corresponds to the exact solution for $\Sigma_{AC}$, which is well-approximated by $\Sigma_{AC} = \Delta t K^{(1)}(1 - P_{ON}^{(1)})/(1 + \Delta t K^{(1)}(1 - P_{ON}^{(1)}))$, where $\Delta t = 10$ s corresponds to the data sampling time. (F) Correlation time and correlated noise magnitude are properly retrieved after single allele deconvolution. Using the Gillespie algorithm, we generated simulated data ($N = 200$ and 50 min long cell recordings per condition) according to the two-state model in D. For each input condition ($P_{ON}^{(1)}$ from 0.03 to 0.9 and $T_C^{(1)}$ from 0.5 to 10 min), we performed single allele deconvolution and computed the auto-correlation on the resulting transcription rates. We estimated the correlation time $T_C^{(1)}$ and the magnitude $\Sigma_{AC}$ by fitting exponential. Both parameters are properly retrieved with minimal biases. Color code stand for $P_{ON}^{(1)}$ and dashed line for slope 1. (G) Estimating deconvolution biases due to elongation rate measurement bias. As in F, we generated simulated data ($P_{ON}^{(1)}$ from 0.03 to 0.9 and $T_C^{(1)} = 2$ min) and aimed to deconvolve the data with elongation rate higher (orange dot) or lower (yellow dot) than the correct value (blue dot, used to generate the data as in F). Overall, the estimated parameters are estimated correctly, with larger biases at large $P_{ON}^{(1)}$ and when underestimating the elongation rate (yellow).
A

<figure><img src="image_32.png" /></figure>

B

<figure><img src="image_33.png" /><figcaption>FIG. S4. Transcription rate explains dynamic pattern establishment. (A) A simple modeling attempt for protein accumulation from mean transcription rate measurements. The mean transcription rate (left column) across space and time is estimated by normalizing the measured mean activity by the elongation time and applying a minor correction for the delay (< 1 min) resulting from the loop insertion location. Horizontal white dashed lines correspond to the transition (mitosis) from NC13 to NC14. Protein accumulation (middle column) is computed from the mean transcription rate as the convolution of the latter with a kernel modeling protein decay, diffusion, and delay due to mRNA export, translation, and nuclear import. This simple model introduces three free parameters, a protein lifetime, a diffusion constant, and a time-delay (see B). These three parameters were set by minimizing the mean squared error with previously measured protein patterns (right column; Dubuis et al., 2013). Small residual deviations between middle and right columns might be due to post-transcriptional regulatory processes that our simple model does not account for. (B) Parameters estimated for the modeled accumulation of effective proteins as described in A. The three parameters were either estimated for each gene separately (color bars) or all genes together (dashed lines, used for middle column in A). Overall, the effective parameters are mostly in line with previous estimates [64].</figcaption></figure>
A

<figure><img src="image_34.png" /></figure>

B

<figure><img src="image_35.png" /></figure>

C

<figure><img src="image_36.png" /></figure>

D

<figure><img src="image_37.png" /></figure>

E

<figure><img src="image_38.png" /></figure>

FIG. S5. Validation of transcriptional parameter estimation using stationary simulated data. Burst calling permits effective estimation of transcriptional parameters. For each set of input parameters ($K^{(1)}$, $T_C^{(1)}$, $P_{ON}^{(1)}$), we generated simulated data (200 alleles and 50 min long time series at 10 s intervals and with the same imaging noise as measured in real data). Using the two-state model (Fig. S3D) and the Gillespie algorithm, we generated time series for each individual sister chromatid. We summed the initiation events of both chromatids assuming independence. The resulting initiation events were convolved with the elongation kernel to generate synthetic single allele signal data over which measurement noise was added. We performed single allele deconvolution and burst calling as described (Fig. 3A) to estimate the effective mean transcription parameters ($R$, $K$, $P_{ON}$, $T_{ON}$, $T_{OFF}$ and $T_C$) for each simulated set of alleles. (A-B) Comparison between the theoretically expected (top row) and the estimated from burst calling effective parameters (bottom row) as a function of time, using stationary (constant in time) input parameters ($T_C^{(1)}$ = 2 min, $K^{(1)}$ = 8 mRNA/min and $P_{ON}^{(1)}$ varies from 0.03 to 0.9 that is color coded from blue to orange). Overall, the estimated effective parameters are recovered very well. (A) Despite sampling fluctuations, the constancy of the parameters in time is properly preserved (bottom row), as it should be based on input (top row). (B) Despite the stationary nature of the input parameters $T_C^{(1)}$ and $P_{ON}^{(1)}$, biases for $T_{ON}$, $T_{OFF}$ and $T_C$ (top row) are expected due to the finite length of the simulated time series (censoring), especially noticeable near the beginning and the end, where the estimations are “bent” (the resulting mean time estimate cannot exceed the width of the time window used to perform the estimate). However, our parameter estimation (bottom row) is very much in line with the expected biases (top row). (C) Global comparison between expected and estimated parameters for the stationary case (input parameters are constant in time). The parameter estimation was performed on a large simulated data set that includes data in A and B ($T_C^{(1)}$ = 2 min, $P_{ON}^{(1)}$ from 0.03 to 0.9) and data for other values of $T_C^{(1)}$ comprised between $T_C^{(1)}$ = 0.5 (blue dots) to $T_C^{(1)}$ = 10 (red dots). Each dot results from one combination of input parameters ($K^{(1)}$, $T_C^{(1)}$, $P_{ON}^{(1)}$) and corresponds to the median effective parameter and the error bars to the 68% confidence interval estimated over 50 min. Our deconvolution and burst calling approach leads to an excellent estimation of the effective parameters over a large range of $T_C^{(1)}$ and $P_{ON}^{(1)}$ values, albeit with noticeable biases in $P_{ON}$, $K$ and $T_{ON}$ when $T_C^{(1)}$ approaches 0.5 min (blue dots). Importantly, biases for the effective switching correlation time $T_C$ are small, supporting our ability to detect its constancy in real data. (D) Summary of median relative error for each effective parameter estimated from the data in C as a function of input $T_C^{(1)}$. Parameter estimated from real data (Fig. S9B) suggests that $T_C^{(1)}$ lies within 1 and 3 min (black border). (E) Summary of median relative error for each effective parameter estimated from the whole data set in C as a function of the burst calling parameters. Burst calling depends on two free parameters: the time window $W$ over which the rate is estimated, and the rate threshold $R_d$ applied to call the burst (see Fig. 3A). Our default parameter values are $W$ = 1 and $R_d$ = 2, which should be close to optimal given the estimated correlation time of $\tau_{AC}$ ~ 1 min (Fig. 2F) and our measurement sensitivity of 1-2 mRNA. When testing the effect of different $W$ and $R_d$ values on the median relative error, our default choice leads to the lowest global relative error.
A

<figure><img src="image_39.png" /></figure>

C

<figure><img src="image_40.png" /></figure>

E

<figure><img src="image_41.png" /></figure>

G

<figure><img src="image_42.png" /></figure>

B

<figure><img src="image_43.png" /></figure>

D

<figure><img src="image_44.png" /></figure>

F

<figure><img src="image_45.png" /></figure>

H

<figure><img src="image_46.png" /></figure>

FIG. S6. Validation of transcriptional parameter estimation using *hb*-like non-stationary simulated data. Simulated data is generated as in Fig. S5, except that the input parameters are non-stationary (time-dependent). (A-B) Comparison between expected (A, B, top row) and estimated effective parameters (A, B, bottom row) for time-dependent input parameters (A, B, top row), mimicking the transcriptional output of *hb* in NC14. Color code (blue to orange) stands for virtual AP position. Estimated effective parameters tightly recapitulate the time-dependence of the expected parameters, allowing us to also capture the temporal regulation of the gap genes. (C, D) We show with this realistic test case that our estimated parameters (A, B, bottom row) using burst calling very closely match the expected effective ones (A, B, top row). Indeed, most of the time points (circles) lie on the line of slope one (dashed line), which strongly supports our ability to precisely characterize the transcriptional parameters from real data. (E-F) Comparison between input (E, F, top row) and burst calling estimated single-gene copy parameters (E, F, bottom row) as a function of time. The time-dependent input parameters (E, F, top row) were used to generate the synthetic *hb* data in A and B. Color code (blue to orange) stands for virtual AP position as in A and B. The estimated single-gene copy parameters were computed from the effective ones (A, B, bottom row) assuming the latter originate from two independent sister chromatids. Namely, we get $R^{(1)} = R/2$, $P_{\text{ON}}^{(1)} = 1 - (1 - P_{\text{ON}})^{1/2}$, $K^{(1)} = K(1 + (1 - P_{\text{ON}})^{1/2})/2$, which are exact, and assuming steady state $T_{\text{C}}^{(1)} = 2T_{\text{C}}/(1 + (1 - P_{\text{ON}})^{1/2})$, $k_{\text{ON}} = P_{\text{ON}}^{(1)}/T_{\text{C}}^{(1)}$ and $k_{\text{OFF}} = (1 - P_{\text{ON}}^{(1)})/T_{\text{C}}^{(1)}$. (G-H) Even though the single gene copy parameters are deeply buried in the data, our simple burst calling procedure still manages to recover them correctly. Some discrepancies are observed for $k_{\text{ON}}$, $k_{\text{OFF}}$ and $T_{\text{C}}^{(1)}$, mostly near the beginning when the transient after mitosis violates the near-steady-state assumption. But these are expected as the relationships $k_{\text{ON}} = P_{\text{ON}}^{(1)}/T_{\text{C}}^{(1)}$ and $k_{\text{OFF}} = (1 - P_{\text{ON}}^{(1)})/T_{\text{C}}^{(1)}$ are only valid for near-steady state.
<figure><img src="image_47.png" /><figcaption>FIG. S7. Single gene copy parameter and transient for *hunchback* in NC14. [Caption see next page.]</figcaption></figure>
FIG. S7. Single gene copy parameter and transient for hunchback in NC14. (A) Transcription rate *R* versus the product of the initiation rate *K* and the ON-probability *P*<sub>ON</sub> for *hb* in NC14 at all time points and positions. The color code stands for AP position as in Fig. 3 and 4. As it should be by construction, *R* can be decomposed into the product of *K* and *P*<sub>ON</sub>. (B) Log-derivative with respect to *P*<sub>ON</sub> of the controlling parameter *K* (gray) and *P*<sub>ON</sub> (blue) as a function of *P*<sub>ON</sub>. Since *R* = *K* · *P*<sub>ON</sub>, d log *K*/d *P*<sub>ON</sub> + d log *P*<sub>ON</sub>/d *P*<sub>ON</sub> = d log *R*/d *P*<sub>ON</sub> and changes in *R* are thus dominated by changes in *P*<sub>ON</sub> when d log *P*<sub>ON</sub>/d *P*<sub>ON</sub> = 1/*P*<sub>ON</sub> > d log *K*/d *P*<sub>ON</sub>. The log-derivative is computed from the mean trend in data (Fig. 3I) and the transition is at *P*<sub>ON</sub> = 0.75 (black dashed line). (C) Fraction of changes in *R* explained by *P*<sub>ON</sub> as a function of *P*<sub>ON</sub>. Black line shows d log *P*<sub>ON</sub>/d *P*<sub>ON</sub> normalized by the sum d log *K*/d *P*<sub>ON</sub> + d log *P*<sub>ON</sub>/d *P*<sub>ON</sub> (see B). The solid black line is computed as in B, whereas the dotted black line is the corrected contribution assuming two independent sister chromatids (see G and H). Most of the changes in *R* are thus mostly dictated by *P*<sub>ON</sub>. (D-F) Single gene copy (sgc) parameter computed from the estimated effective parameters for *hb* in NC14, assuming two independent sister chromatids. The color code stands for AP as A. The spatiotemporal regulation of the sgc transcription rate *R*<sup>(1)</sup> (D) and the sgc ON-probability *P*<sub>ON</sub><sup>(1)</sup> (E) are very similar to their corresponding effective parameters, *R* and *P*<sub>ON</sub>. On the other hand, the sgc initiation rate *K*<sup>(1)</sup> (F) no longer varies across position, but only as a function of time. During NC14, the sgc initiation rate *K*<sup>(1)</sup> decreases by 38%, from 8.0 mRNA/min (first dotted line) to 5.0 mRNA/min on average (third dotted line), with most of the decrease happening between the 16<sup>th</sup> and 34<sup>th</sup> minute mark (second dotted line at 6.1 mRNA/min). (G-I) Temporal changes in sgc initiation rate and independent sister chromatid assumption explain the dependence of the effective initiation on the ON-probability. Color code stands for three time-windows in NC14: early (cyan, 2.5-16.7 min), mid (purple, 16.7-34.2 min), and late (magenta, 34.2-50 min). (G) Most of the variation in *K*<sup>(1)</sup> is explained by time, rather than *P*<sub>ON</sub><sup>(1)</sup>. The dotted lines are drawn at the same *K*<sup>(1)</sup> values as in C. (H) The sgc transcription rate *R*<sup>(1)</sup> appears almost linearly dependent on *P*<sub>ON</sub><sup>(1)</sup>. The nonlinearity is mostly explained by temporal changes in *K*<sup>(1)</sup>, as highlighted by dotted lines whose slopes are the *K*<sup>(1)</sup> values from F and G. (I) Under the two independent sister chromatids assumption, the effective initiation rate *K* depends on *P*<sub>ON</sub> and on *K*<sup>(1)</sup>, which varies as a function of time (see F). As *P*<sub>ON</sub> increases, the propensity to observe two gene copies initiating transcription at the same time increases, which explains up to a factor of two in the dependence of *K* on *P*<sub>ON</sub>. Indeed, the dotted lines correspond to the predicted behavior using the same three constant values of *K*<sup>(1)</sup> as in F. In addition, *K*<sup>(1)</sup> varies by up to 38% along time during NC14, as can be seen in F. Together, it explains close to a factor of 3.2 in *K* variation with *P*<sub>ON</sub>. This reasoning explains the observed relationships very well (black line) and argues for a very weak dependence of *K* on *P*<sub>ON</sub>. (J) Close-up of the first 10 min of Fig. 3D shows a rapid transient in the ON-probability for *hb* at the onset of NC14. The vertical dashed line at 7.5 min marks the transition between the transient and the near-steady-state regime as observed in H. (K) Evidence for out-of-steady state transient for *hb* in NC14. Near the steady-state regime, *P*<sub>ON</sub> should be well-approximated by *T*<sub>ON</sub>/*(T*<sub>ON</sub> + *T*<sub>OFF</sub>)*, as it is the case beyond the 7.5 min mark in NC14 (gray circles tracing the diagonal; see also Fig. 4D). However, at the post-mitotic onset of transcription (~3-7.5 min into NC14) we observe strong deviations from the expected near-steady state relationship at all positions (color curves). The system is undergoing a fast transient relaxation that drives it near-steady state within the first 8 min after mitotic exit.
A

gt ♂ anterior

gt ♂ posterior

<figure><img src="image_48.png" /></figure>

B

<figure><img src="image_49.png" /></figure>

C

kr

<figure><img src="image_50.png" /></figure>

gt ♀

<figure><img src="image_51.png" /></figure>

D

<figure><img src="image_52.png" /></figure>

E

<figure><img src="image_53.png" /></figure>

FIG. S8. **Further transcriptional parameters collapse.** (A) Kymograph of ON-probability for *gt* transcription in male embryos, as a function of position and time. As for the other gap genes, the spatiotemporal transcription patterns arise from complex regulation of the ON-probability. (B) Transcriptional parameters for *gt* in male embryos for NC13 and NC14 as a function of PON (orange data points) and all the other gap gene data sets (gray). (C) Initiation rate *K* as a function of PON. *K* collapses for all gap genes across time and position. Colored data points represent individual gap genes (same color code as in Fig. 5A-E and S8A-B); underlying is the remaining data of all other genes (gray). The *K*–PON relationship for *hb* is shown in Fig. 3H. (D) Near-steady state relationship between ON-probability PON and TON/(TON + TOFF) for all gap genes in NC13 (*t* ≥ 6.5min) and NC14 (*t* ≥ 7.5min). Although the data is distributed near the expected relationship (dash line), we observe a slight but clear bias at the extreme ends of the PON spectrum (solid line), namely TON/(TON + TOFF) is slightly above zero at PON = 0 and slightly below one at PON = 1. This is a consequence of the finite nature of our recording (50 min in NC14 and 18.4 min in NC13). Thus, our recording time sets an upper limit on the length of measurable ON and OFF intervals. That limit leads to the observed bias. (E) Global scatter of the burst size *B* and burst frequency *F* (color) as a function of PON for all gap genes in both NC13 and NC14 and putative accessible space (gray region). See definitions of *B* and *F* at top of respective panels. Color code and accessible space defined as in Fig. 5F. Solid line stands for the bursting rules derived with bias (solid line in C) and dashed line without.
<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

FIG. S9. Common bursting rules across developmental time and for single gene copy. (A) Global scatter of effective transcriptional parameters as a function of $P_{\text{ON}}$ for all gap genes estimated within three time windows, corresponding to NC13 ($6.5 \le t \text{ min}$) plus early NC14 ($7.5 \le t < 20.5 \text{ min}$), mid NC14 ($20.5 \le t < 34.5 \text{ min}$), and late NC14 ($34.5 \le t < 48 \text{ min}$). The observed bursting relationships are further refined when accounting for possible temporal changes (color scatter) compared to all time-pooling (gray scatter, Fig. 5F). Indeed, small changes in $K$ and $T_c$ over developmental time ($\sim 40\%$ decrease) explain part of the observed spread in Fig. 5F. (B) As in A, but for single gene copy parameters computed from the effective ones assuming independent sister chromatids. Interestingly, the relationship between sgc transcription rate $R^{(1)}$ and sgc ON-probability $P_{\text{ON}}^{(1)}$ is almost linear, confirming that the sgc initiation rate $K^{(1)}$ does not depend strongly on $P_{\text{ON}}^{(1)}(1)$. Thus, the apparent dependence of $K^{(1)}$ on $P_{\text{ON}}^{(1)}(1)$ is only effective and results from measuring two sister chromatids (two gene copies) together, instead of an isolated single gene copy.
<figure><img src="image_56.png" /><figcaption>FIG. S10. Bursting rules verified by cis- and trans-perturbations. (A) Distal *kni* enhancer removal. The MS2-stem loops are inserted at the same location in the mutant (enhancer deletion) and wild-type fly lines. (B) Quantification of *kni* wild-type and mutant (A) phenotypes. Both transcription rate *R* as a function of *x*/L (left) and the kymograph for *P*<sub>ON</sub> (right) display a significant level decrease underlying the expression patterns of the mutant. Dotted arrow indicates time point in kymograph at which rate profiles (left) are depicted. (C) Transcription parameters for *kni* cis-mutant (olive) collapse on corresponding wild-type parameters (gray), as for *hb* cis-mutant (see Fig. 6C). Solid black lines correspond to the endogenous bursting rules from Fig. 5F. (D) *kni* measurements in a *hb* null background. The absence of *hb* expression alters the network, namely the concentration of input transcription factors sensed by *kni* in the mutant fly line. (E) Quantification of *kni* wild-type and mutant (D) phenotypes. Both transcriptional rate *R* (left) and *P*<sub>ON</sub> kymograph (right) display a significant shift of the anterior boundary in the mutant expression patterns. Dotted arrow indicates time point in kymograph at which rate profiles are depicted.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>FIG. S11. Bursting rules predict ON and OFF modulation by *cis-* and *trans*-perturbations. (A-C) Predicted $T_{OFF}$ versus $T_{ON}$ bursting modulation for mutant based on wild-type-derived rules (Fig. 5F, black lines). (A) The type of modulation is predicted by first approximating $T_{OFF}$ and $T_{ON}$ as a function of $P_{ON}$ using the wild-type rules. The predicted fold change in $T_{OFF}$ ($T_{OFF}^{mut}/T_{OFF}^{wt}$) and $T_{ON}$ ($T_{ON}^{mut}/T_{ON}^{wt}$) are then computed for all possible pairs of $P_{ON}$ (i.e. $P_{ON}^{wt}$ and $P_{ON}^{mut}$). The dotted line delimits the regions where changes in transcription rate are either dominated by changes in $T_{OFF}$ (gray region, $|\log(T_{OFF}^{mut}/T_{OFF}^{wt})| > |\log(T_{ON}^{mut}/T_{ON}^{wt})|$) or $T_{ON}$ (blue region, $|\log(T_{OFF}^{mut}/T_{OFF}^{wt})| < |\log(T_{ON}^{mut}/T_{ON}^{wt})|$). The solid black lines delimit the region, where changes in $T_{OFF}$ and $T_{ON}$ are not significant given the “thickness” of our relationships (95% confidence intervals, see Methods). Thus, this procedure defined a look-up table enabling prediction of the type of modulation using pairs of $P_{ON}$. (B) Scatter plot of all the $P_{ON}$ pairs from *hb* wt and *cis*-mutant (at same spatiotemporal location). Colors correspond to the predicted modulation ($T_{OFF}$ dominated in gray and $T_{ON}$ dominated in blue) using the look-up table in A. (C) Verification of predicted modulation in *B* (color code as in B). For each $P_{ON}$ pair, we computed the $T_{OFF}$ ($T_{OFF}^{mut}/T_{OFF}^{wt}$) and $T_{ON}$ ($T_{ON}^{mut}/T_{ON}^{wt}$) fold change using the estimated $T_{OFF}$ and $T_{ON}$ from data (Fig. 6C). Supporting our ability to predict the modulation, almost all the blue data points (predicted as $T_{ON}$ modulation) are located above the slope 1 diagonal (dashed line), whereas most of the gray ones (predicted as $T_{OFF}$ modulation) are below. Thus, for most data points ($> 85\%$) the prediction is correct (Fig. 6F). (D-F) Initiation rate *K*, burst size *B*, and burst frequency *F* for *hb* in *cis*-mutation (cyan), *kni* in *cis*-mutation (olive) and for *kni* in *trans*-mutation (light green) collapse on corresponding wild-type parameters (gray). Solid black lines correspond to the endogenous bursting rules from Fig. 5F and S8E. (G-I) Predicted *F* versus *B* modulation for mutant based on wild-type-derived rules (Figure S8E, black lines). (G) *F* and *B* are first approximated as function of $P_{ON}$ using the wild-type rules. The predicted fold change in $F$ ($F^{mut}/F^{wt}$) and $B$ ($B^{mut}/B^{wt}$) are then computed for all possible pairs of $P_{ON}$ (i.e. $P_{ON}^{wt}$ and $P_{ON}^{mut}$). The dotted line delimits the regions where changes in transcription rate are either dominated by changes in *F* (gray region, $|\log(F^{mut}/F^{wt})| > |\log(B^{mut}/B^{wt})|$) or *B* (blue region, $|\log(F^{mut}/F^{wt})| < |\log(B^{mut}/B^{wt})|$). As for $T_{OFF}$ and $T_{ON}$ (in A), the solid black lines delimit the region of significance and we have thus defined a look-up table enabling prediction of *F* vs *B* modulation using pairs of $P_{ON}$. (H) Scatter plot of all the $P_{ON}$ pairs from *hb* wt and *cis*-mutant (at same spatiotemporal location). Colors correspond to the predicted modulation (*F* dominated in gray and *B* dominated in blue) using the look-up table in G. (I) Verification of predicted modulation in *H* (color code as in H). For each $P_{ON}$ pair, we computed the $F$ ($F^{mut}/F^{wt}$) and $B$ ($B^{mut}/B^{wt}$) fold change using the estimated *F* and *B* from data (E and F). Supporting our ability to predict the modulation, almost all the blue data points (predicted as *B* modulation) are located above the slope 1 diagonal (dashed line), whereas most of the gray ones (predicted as *F* modulation) are below. Thus, for most data points ($> 95\%$) the prediction is correct (Figure 6G).</figcaption></figure>
A

<figure><img src="image_58.png" /><figcaption>FIG. S12. Validating bursting rules with smFISH data. (A) Bursting relationships in NC13 and early NC14 are consistent with parameters inferred from previous smFISH measurements. We converted the single gene copy parameters of the gap genes in mid-late NC13 (Zoller et al.) into effective parameters for two sister chromatids (color dots). Overall, these parameters closely verify our relationships derived from live measurements.</figcaption></figure>