Synthetic deconstruction of hunchback regulation by Bicoid

During development, cell identity is established reproducibly among individuals through the expression of specific genes at the correct time and correct location in space. How genes extract and combine both positional and temporal information from different transcription factor (TF) profiles along polarity axes remain largely unexplored. Here, we showcase the classic hunchback gene in fruit fly embryos, with focus on 3 of its main TFs: Bicoid, Zelda and Hunchback proteins. We constructed a series of synthetic MS2 reporters, where the numbers and combination of binding sites for each TF are varied. Using live imaging of transcription dynamics by these synthetic reporters and modeling tools, we show that i) a Bicoid-only synthetic reporter needs 3 more Bicoid binding sites than found in the hunchback promoter to recapitulates almost all spatial features of early hb expression but takes more time to reach steady state; ii) Hunchback and Zelda binding sites combined with Bicoid sites both reduce the time to reach steady state and increase expression at a different step in the activation process: Zld sites lower the Bicoid threshold required for activation while Hb sites increase the polymerase firing rate and reduce bursting; iii) the shift of the Bicoid-only reporter induced by a reduction by half of Bicoid concentrations indicates that the decay length of the Bicoid activity gradient is lower than the decay length of the Bicoid protein gradient. Altogether, this work indicates that Bicoid is the main source of positional information for hunchback expression and places back the Bicoid system within the physical limits of an equilibrium model.


Introduction
Morphogen gradients are used by various organisms to establish polarity along embryonic axes or within organs. In these systems, positional information stems from the morphogen concentration detected by each cell in the target tissue and mediates the determination of cell identity through the expression of specific sets of target genes. While these processes ensure the reproducibility of developmental patterns and the emergence of properly proportioned individuals, the question of whether the morphogen itself directly contributes to this robustness or whether it requires the involvement of downstream cross-regulatory networks or cell-communication remains largely debated. This question becomes even more pressing with the recent discovery that when studied at the single cell level, transcription is frequently observed to be an extremely noisy process, hardly suggestive of such precise control.
Here, we chose a synthetic approach to decipher the functioning of Bcd in the transcription process at the mechanistic level. We built Bcd-only reporters with specific numbers of Bcd BS as well as reporters with 6 Bcd BS in combination with BS for the two known maternal Bcd co-factors binding to the hb-P2 promoter, namely the Hb protein itself Simpson-Brose et al., 1994) and the Zld pioneer transcription factor (Hannon et al., 2017;Xu et al., 2014). We show that 6 Bcd BS are not sufficient to recapitulate the hb-P2 expression dynamics while a reporter with only 9 Bcd BS recapitulates most of its spatial features, except a slightly lower steepness of its expression boundary and a longer period to reach steady state. To account for the bursty behavior of Bcd-only reporters in excess of Bcd, we fitted our data to a model involving a first step of Bcd binding/unbinding to the BS array and a second step where the bound Bcd molecules activate transcription. Synthetic reporters combining Bcd BS with either Hb or Zld BS indicated that both Hb and Zld sites reduce the time to reach steady state and increase expression by different means: Zld sites contribute to the first step of the model by drastically lowering the Bcd concentration thresholds required for activation while Hb sites act in the second step by reducing Bcd-induced burstiness and increasing the polymerase firing rates. Lastly, the boundary shift of the Bcd-only synthetic reporter with 9 Bcd BS in embryos maternally expressing one (1X) vs two (2X) bcd functional copies was smaller than the expected shift given the decay length of the Bcd concentration gradient. Thus, we propose the existence of an effective Bcdactivity gradient underlying the Bcd-protein gradient. Altogether, this work reinforces the notion that the Bcd gradient is the main source of positional information for the early expression of hb while both Zld and Hb accelerate the process.

Results
Nine Bicoid binding sites alone recapitulate most features of the hb-P2 pattern We first investigated the transcription dynamics of Bcd-only MS2 reporters carrying exclusively 6, 9 or 12 strong Bcd binding sites (BS) (Hanes and Brent, 1989;Treisman et al., 1989) upstream of an hsp70 minimal promoter (Table S1), all inserted at the same genomic location. Movies were recorded and analyzed from nuclear cycle 11 (nc11) to 13 (nc13) but we focused on nc13 data, which are statistically stronger given the higher number of nuclei analyzed. Unless otherwise specified, most conclusions were also valid for nc11 and nc12.
The expression of the B6 (6 Bcd BS), B9 (9 Bcd BS) and B12 (12 Bcd BS) reporters ( Fig. 1A) harbored similar features as expression of the hb-P2 reporter (Lucas et al., 2018), which carries the ~300 bp of the hb-P2 promoter and the hb intron (Fig. 1B, Table S1): during the cycle, transcription was first initiated in the anterior with the expression boundary moving rapidly towards the posterior to reach a stable position into nc13 (Fig. 1C). For all synthetic reporters, the mean onset time of transcription following mitosis (Fig. 1D) showed a dependence on position along the AP axis, as observed for hb-P2 (Lucas et al., 2018). Thus, Bcd concentration is a rate-limiting factor for the expression of all reporters. As indicated by the distributions of onset time in the anterior (~-30 %EL), the first transcription initiation time at high Bcd concentration were not statistically different (p-values > 0.5) for all synthetic reporters (B6, B9 or B12) and hb-P2 (Fig. 1E). In contrast, in the middle of the axis where Bcd is limiting, transcription dynamics of the various reporters was quite diverse (Fig. 1F): the time it took for the hb-P2 reporter to reach the final decision to position its boundary (converging time, Table S2) was only 225 ± 25 s while it took about twice as much time for B6 (425 ± 25 s) or B9 (475 ± 25 s) and slightly less for B12 (325 ± 25 s).
For all reporters, the fraction of nuclei with MS2 signal during the cycle exhibited a sigmoid-like pattern along the AP axis reaching 100% in the anterior and 0% in the posterior (Fig. 1G). We fitted these patterns with sigmoid functions of position along the AP axis and extracted (see Materials & Methods, section C for derivation) quantitative values for the position and width of the expression boundary (Fig.  1E). Increasing the number of Bcd BS from 6 to 9, shifted the expression boundary towards the posterior and decreased the width of the boundary (Fig. 1H) whereas increasing the number of Bcd sites from 9 to 12 did not significantly change the boundary position nor the boundary width. Of note, B9, B12 and hb-P2 expression boundaries were at almost identical positions while the width of the hb-P2 boundary was smaller than the width of the B9 or the B12 boundaries (Fig. 1H).
Thus, even though 6 Bcd BS have been described in the hb-P2 promoter, having only 6 Bcd BS alone in a synthetic reporter is not sufficient to recapitulate the hb pattern. Increasing this number up to 9 is sufficient to recapitulate almost all spatial features of the hb-P2 pattern except for the steepness of the expression boundary. Of note, the Bcd-only reporters take much longer than the hb-P2 reporter to reach the final decision for boundary positioning suggesting that binding of additional transcription factors in the hb-P2 promoter likely contribute to speeding-up the process.

Bicoid-dependent transcription is bursty at steady state even in excess of Bicoid
To study the kinetics of transcription induced by Bcd, we compared the dynamics of transcription of the hb-P2 and the Bcd-only reporters at steady state (in the time window of 600-800s). From the time trace of MS2 activity in each nucleus, the fluctuation of the transcription process (burstiness) at a given position along the AP axis was featured by , the average fraction of the cycle length during which fluorescent spots were observed ( Fig. 2A). In the anterior (~-30 %EL), increased when increasing the number of Bcd BS in synthetic reporters from 6 to 9, with ( 6) = 0.47 ± 0.02 and ( 9) = 0.80 ± 0.07.
(hb-P2) = 0.84 ± 0.008 was as high as for B9 or B12 ( ( 12)= 0.76 ± 0.07). These values were all smaller than the fraction of expressing nuclei (~ 1, Fig. 1G). This indicated bursty transcription activity in individual nuclei for all reporters, as confirmed by their individual MS2 traces in this region. Interestingly, for all Bcd-only reporters reach a plateau in the anterior where the Bcd concentration is in excess ( Fig. 2A). Thus, as in this region the Bcd BS on those reporters are likely to be always occupied by Bcd molecules, the burstiness observed is not caused by the binding/unbinding of Bcd to the BS array but by downstream processes. Meanwhile, the mean intensity of the MS2 signals ( ) in the anterior region did not vary between reporters (Fig. 2B), suggesting that the number of bound Bcd molecules does not regulate the RNAP firing rate within transcription bursts.

A model to recapitulate expression dynamics from Bicoid-only synthetic reporters
To explain the observed dynamics of the expression patterns ( Fig. 1C) and bursty transcription in regions with excess Bcd ( Fig. 2A), we built a model for transcription regulation of the Bcd-only synthetic reporters (Fig. 2,. In this model, regulation occurs in two steps: first, nuclear Bcd molecules can bind to and unbind from the Bcd BS on the promoter (Fig. 2C) and second, bound Bcd molecules can activate transcription (Fig. 2D).
In step 1 (Fig. 2C), the binding and unbinding of Bcd to an array of identical Bcd BS were modeled explicitly, as in (Estrada et al., 2016;. In our model, the binding state is denoted by , with the number of bound Bcd molecules ( ≤ ). The binding rate constants depend on the number of free BS ( − + 1) and the Bcd search rate for a single BS . The unbinding rate constants were varied to account for various degrees of Bcd-DNA complex stability and binding cooperativity. In step 2 (Fig. 2D), we expanded this model to account for the burstiness in transcription uncoupled with Bcd binding/unbinding (Fig. 2B). The promoter dynamics was modeled as a two-state model, ON and OFF, to account for the observed bursts of transcription with a moderate time scale (between 10 s and 100 s (Bothma et al., 2014;Desponds et al., 2016;Lammers et al., 2019)). The turning ON rate ( ) was modulated by the number of bound Bcd molecules. When the BS arrays had less than Bcd molecules ( ≥ 0), transcription could not be activated ( ( ) = 0). To account for the uncoupling between the burstiness of transcription and the Bcd binding and unbinding, the turning OFF rate did not depend on the Bcd BS state. When the promoter is ON, RNAP could initiate transcription and be fired at rate . At any given time and nuclei position along the AP axis, we could calculate the probability for the promoter to be in the ON state (see Materials & Methods section D).
In this model, each kinetic parameter could be tuned independently to control the measured transcription dynamics features: Bcd binding rate constants ( , ) controlled the pattern boundary position, Bcd unbinding rate constants ( ) controlled the pattern steepness (Estrada et al., 2016;, the activation/deactivation rates ( , ) controlled the fraction of active loci during steady state ( ), and the RNAP firing rate ( ) controlled the mean loci intensity ( ). The fit of the model for each Bcd-only reporter with the data is shown Fig. 2 (E-G). By comparing features of the transcription dynamics among Bcd-only reporters, we identified which parameters and consequently which processes were dependent on the number of Bcd BS. We started with the parameters allowing a best fit of the model with the B6 data and allowed each of the parameters to vary, either alone or in combination, to fit the B9 and B12 data. These simulations indicated that very good fits could be obtained for B9 and B12 by allowing only 3 of the parameters to vary ( , & ) ( Fig. S4) while the other parameters remained those identified for B6. Also, as they have more Bcd BS than B6, the fitting of the B9 and B12 data to the model generated new parameters (i.e. ( )) for N>6).
Given that the expression patterns of hb-P2 and all Bcd-only reporters reached a plateau in the anterior where Bcd concentration is likely in excess, we compared the activation rates ( ) of the promoter when = 6, 9 or 12 Bcd BS were occupied. Assuming that the number of bound Bcd proteins did not affect the switch OFF rate , we found a fold change of ~4.51 between ( ) and ( ). This fold change is 3 times greater than the ratio of the Bcd BS numbers between B9 and B6, arguing against independent activation of transcription by individual bound Bcd TF. This analysis detailed in Materials & Methods section G provides evidence for synergistic effects between several bound Bcd molecules.

Hunchback reduces the burstiness of Bicoid-dependent transcription
Despite the same number of Bcd BS in the B6 and hb-P2 reporters, their expression pattern and dynamics were very different ( Fig. 1 and Fig. 2A). To determine whether this difference could be explained by the presence of BS for other TFs in the hb-P2 promoter (Fig. 1A), we used our synthetic approach to decipher the impact on the various features highlighted in our model when adding to the reporters BS for the two major partners of Bcd, Hb and Zld also present in hb-P2 promoter (Fig. 3A). As our goal was to determine to which mechanistic step of our model each of these TF contributed, we purposefully started by adding BS numbers that are much higher than in the hb-P2 promoter.
We first analyzed the impact of combining the Bcd BS with Hb BS. A H6 reporter containing only 6 Hb BS did not exhibit any MS2 signal from nc11 to nc13 (not shown). This indicated that the Hb protein alone, even with an abundance of Hb sites, could not activate transcription on its own. When combining 6 Hb BS with the 6 Bcd BS of B6 (henceforth named the H6B6 reporter, Fig. 3A), expression was detected in a similar domain to that of the B6 reporter, albeit with much higher fraction of active loci at any given time during the cycle (Fig. 3B, middle panel). Across the embryo AP axis, the mean onset time of transcription after mitosis with the H6B6 reporter was not changed (p-values > 0.5) when compared to B6 (Fig. 3C) and in the anterior region (excess of Bcd), the cumulative distribution of onset time was the same (Fig. 3D). Interestingly, at their respective boundary positions where the Bcd concentration is limiting, the presence of Hb BS reduced to 325 ± 25 s the time required for the synthetic H6B6 reporter to reach the final decision to position its boundary (purple dashed line in Fig. 3E and converging time, Table S2) when it was 425 ± 25 s for B6 (yellow dashed line in Fig. 3I and converging time, Table S2). For H6B6, the fraction of nuclei with MS2 signal during the cycle exhibited a sigmoid-like pattern ( Fig. 3F) with, when compared to B6, a boundary slightly (only one nucleus length) shifted towards the posterior and a width reduced by half (Fig. G). The kinetics of transcription regulation by the Hb protein was inferred from the fraction of the loci's active time ( ) at steady state. In the anterior region, this fraction was always near-saturation for the H6B6 reporter (~0.95 to 1) (Fig. 3H), with very few nuclei exhibiting bursty expression. This non-bursty behavior of the H6B6 reporter contrasts with the highly bursty expression of B6 reporter. Meanwhile, in the anterior region, the mean fluorescence intensity of active H6B6 loci was at least twice higher than that of all synthetic Bcd-only reporters (Fig. 3I).
To model H6B6 activity, the same formalism, as applied to B9 and B12 reporters, was used starting from the parameter values imposed from the fitted model of B6 and then varying those parameters, either alone or in combination, to fit the H6B6 data. The simulations indicated that a moderate fit to the data was obtained when varying only the and parameters while varying in addition the 3 of the parameters ( , & ) allowed very good fitting of the model to the data (Fig. S4).
Altogether, this suggests that Hb binding to the promoter accelerates the measurement of positional information by Bcd by improving both the unbinding kinetics of Bcd to its BS, which is consistent with the half reduction of the boundary steepness ( Fig. 3G) and the kinetics of activation/deactivation transcription rates, consistent with reduced burstiness (Fig. 3H).

Zelda lowers the Bcd threshold required for expression
As the hb-P2 promoter also contains Zld BS, we used our synthetic approach to investigate the role of Zld in the Bcd system. a reporter with only 6 Zld BS (Z6) was strongly expressed along the whole AP axis (Fig. S3A), we had to reduce the number of Zld BS in our synthetic approach to analyze Zld effect. A Z2 reporter containing only 2 Zld BS did not exhibit any MS2 signal (not shown). The Z2B6 reporter, combining 2 Zld BS with 6 Bcd BS (Fig. 3A), exhibited a very different expression pattern when compared to B6 (Fig. 3B, right panel). This expression pattern also varied with the nuclear cycles likely because of drastic changes in Zld transcriptional activity ( Fig. S3) but for simplicity, we focused here on nc13. The onset time of the Z2B6 reporter was similar in the anterior to those of the B6, H6B6 and hb-P2 reporters (Fig. 3, C-D) but unlike B6, H6B6 and hb-P2 it did not vary along the AP axis (Fig. 3C). This suggest that Zld binding can accelerate Bcd-dependent transcription when Bcd is rate-limiting but has no effect when Bcd is in excess (Fig. 3C). As observed with H6B6, the presence of Zld BS reduced to 300 ± 25 s the time required for the synthetic Z2B6 reporter to reach the final decision to position its boundary (blue dashed line in Fig. 3E and converging time, Table S2) when it was 425 ± 25 s for B6 (yellow dashed line in Fig. 3E and converging time, Table S2).
The most striking feature of the Z2B6 reporter was the drastic posterior shift of its expression boundary (17.5 %EL) when compared to B6 ( Fig. 3C-D). It indicates that the threshold of Bcd concentration required for activation is lowered when two Zld BS are present in the promoter together with 6 Bcd BS. Added to this, the pattern boundary width (Fig. 3G) and in the anterior, both the active loci fraction (Fig. 3H) and the loci intensity ( Fig. 3I) were very similar for the Z2B6 and B6 reporters. Therefore, we hypothesize that adding 2 Zld sites can accelerate and facilitate Bcd binding when Bcd is rate-limiting (i.e. increasing or ) without affecting the remaining parameters ( , , , ). Consistent with this hypothesis, simulations for best fitting of the model to the data, starting from the parameters imposed by B6, indicate that a very good fit of the model to the Z2B6 data is obtained when only varying the Bcd binding rate (Fig. S4).
Altogether, this suggests that Zld binding to the promoter accelerates the measurement of positional information by Bcd by facilitating Bcd binding when it is rate-limiting through an increase of the Bcd binding rate kb, without affecting the kinetics of activation/deactivation transcription rates.
A Bcd-activity gradient different from the Bcd-protein gradient Even though Hb is asymmetrically distributed along the AP axis, the presence of Hb BS in our synthetic reporters had only a very mild effect on the positioning of the expression boundary and thus on measurements of positional information by Bcd. Meanwhile, Zld is homogeneously distributed along the AP axis. Therefore, neither Zld nor Hb could provide strong sources of positional information along the AP axis, leaving the Bcd gradient as the principal source. Our synthetic Bcd-only reporters, which exclusively responded to Bcd, provided the tools to quantitatively and directly test if the position of their expression boundary was exclusively dependent on specific thresholds of Bcd concentration. For this, we reduced the amount of the Bcd protein by half in embryos from females, which were heterozygous for a CRISPR-induced deletion of the bcd gene (bcd) (see Materials & Methods, section A). As the amount of Bcd protein is produced from each bcd allele independently of any other allele in the genome (Liu et al., 2013) and as changing the genetic dosage of bcd in the female leads to proportional changes in both mRNA and protein number in the embryo (Petkova et al., 2014), we assumed that embryos from wild-type females (2X) express quantitatively twice as much Bcd proteins as embryos from bcd/+ females (1X). In such Bcd-2X and Bcd-1X embryos, we compared the fraction of nuclei expressing the B9 reporter along the AP axis as modeled at the top of Fig. 4A. For simplicity, we denoted ( ) the expression pattern in Bcd-2X embryos and ( ) the expression pattern in Bcd-1X embryos, with being the nuclei position along the AP axis.
To quantify the effects of perturbing the Bcd gradient, we first extracted from the experimental data the shift in position Δ( ) between two nuclei columns with the same expression distribution in Bcd-2X embryos ( ( )) and in Bcd-1X embryos ( ( − Δ( )), such that 4B). As all the expression patterns are noisy, we calculated the probability distribution of seeing a given shift (Δ( )| ) for each given position from our data and used a grey-scale log-probability map as a function of and Δ to present our results. An example of the log-probability map for the shift expected if the Bcd concentration was reduced by half at each position is shown at the bottom of Fig. 4A. As expected, the prediction of Δ( ) is most reliable in the boundary region (see Materials & Methods section E for a detailed formulation). From the log-probability map of the shift Δ( ) obtained from expression data of the B9 reporter in Bcd-2X vs Bcd-1X embryos, we observed that the shift in this zone can be described by a constant value Δ( ) = Δ and thus the Bcd gradient measured in this region was exponential. From the data, the best fit value of Δ was found to be 10.5 ± 1.0 %EL (cyan dashed line in Fig. 4C). If the threshold of Bcd concentration measured by the Bcd-only reporters was based on the concentration gradient with a decay length ~20% (Gregor et al., 2007b;Houchmandzadeh et al., 2002), the expected shift of the patterns would be |Δ | = . 2 ~14 %EL (red dashed lines in Fig. 4C). Thus, the fitted shift |Δ | ~10.5 %EL extracted from our data is significantly smaller than the expected shift from our knowledge of the Bcd gradient. Of note, the shift obtained at nc13 was larger than the shift obtained at nc12. However, the short length of nc12 (shorter than the time required for the Bcd-only reporters to reach steady state) likely introduces a bias in those measurements (Fig. S5A).
Since the synthetic Bcd-only reporters are expected to position their boundary at the same threshold of active Bcd concentration in Bcd 2X vs Bcd 1X embryos, we propose that they learn about the nuclei position from a gradient of active protein (effective Bcd gradient) rather than from a gradient of concentration detected by immunofluorescence or with GFP-tagged protein, which might include both active and inactive detectable proteins. The effective gradient highlighted by our analysis would be exponential with an effective decay length = |Δ |/ 2 = 15 ± 1.4 %EL. We projected the decay length for this effective gradient in the model accounting for the pattern dynamics of B9 in Bcd-2X embryos to predict its pattern in Bcd-1X embryos. The predicted patterns from the model (black dashed curves) match well with the data (yellow curves) (Fig. 4B).
Lastly, the comparison of hb-P2 patterns in Bcd-2X and Bcd-1X embryos indicated a shift of 11.0 ± 0.5 %EL of the expression boundary (Fig. 4E). As this value was indistinguishable from the shift |Δ | obtained with data of the B9 reporter which only responds to Bcd (Fig. 4E), we concluded that the measure of positional information by the hb-P2 promoter is based entirely on the effective Bcd gradient with ~15 ± 1.4 %EL and does not involve input from other TF binding to the hb-P2 promoter.
Of note, the shift obtained for the hb-P2 MS2 reporter was significantly larger than the shift of 8% EL described for hb in previous studies using the bcd E1 amorphic allele (Houchmandzadeh et al., 2002;Porcher et al., 2010). To understand this discrepancy, we measured the shift in the boundary positions of our hb-P2 and synthetic MS2 reporters in embryos from wild-type vs bcd E1 /+ females and confirmed that in this genetic background the shift of boundary position was 8% EL ( Fig. 4F and S5B). As the molecular lesion in the bcd E1 allele introduces a premature stop codon downstream of the homeodomain (Struhl et al., 1989), these results suggest that the bcd E1 allele likely allows the expression of a weakly functional truncated protein.
The hb-P2 pattern steepness can be explained by an equilibrium model of concentration sensing Assuming that nuclei extract positional information from an effective Bcd gradient with decay length < , we reassessed the Hill coefficient (denoted as ), which reflects the cooperativity of hb regulation by Bcd (Estrada et al., 2016;Gregor et al., 2007a). For this, we fitted the pattern of expressing nuclei by the hb-P2 reporters (Fig. 1D) to a sigmoid function. We transformed the fitted sigmoid function of position to a Hill function describing the transcription regulation function of hb-P2 by the Bcd protein concentration (see Materials & Methods section C for detailed formulations). Given the sigmoid function obtained from the data, the inferred Hill coefficient is proportional to assumed decay length (black line in Fig. 5A). Taking the observed effective decay length = 15 % EL, we obtain H ~ 5.2 compared to H ~ 6.9 for a decay length = 20 %EL. As the hb-P2 promoter contains only 6 known Bcd BS, the value of = 5.2 for the Hill coefficient inferred assuming the decay length is now within the limit of concentration sensing with 6 Bcd BS ( = 6, dashed horizontal line in Fig. 5A) while the former value ~ 6.9 was not achievable without energy expenditure (Estrada et al., 2016) or positive feedback from Hb protein (Lopes et al., 2011).
To verify whether both the dynamics and sharpness of the hb-P2 expression pattern can be sufficiently explained by an equilibrium model of Bcd concentration sensing via =6 Bcd BS, we fitted our model ( Thus, lowering the decay length of the Bcd gradient to its effective value allows a more reliable fit of the model to the data and places back the Bcd system within the physical limits of an equilibrium model for concentration sensing.

Discussion
In this work, we used a synthetic approach with the focus on cis-regulatory elements, to decipher the kinetics of transcription regulation downstream of the Bcd morphogen gradient in early Drosophila embryogenesis. The number and placement of TF BS in our MS2 reporters are not identical to those found on the endogenous hb promoter. Nevertheless, this synthetic approach combined with quantitative analyses and modeling sheds light on the mechanisms of transcription regulation by Bcd, Hb and Zld. Based on this knowledge from synthetic reporters, we built the equilibrium model of transcription regulation which agrees with the data from the hb-P2 reporter expression. Therefore, we expect that the mechanisms uncovered are relevant to the patterning of the endogenous hb gene.
Expression from the Bcd-only synthetic reporters indicate that increasing the number of Bcd BS from 6 to 9 shifts the transcription pattern boundary position towards the posterior region. This is expected as an array with more BS will be occupied faster with the required amount of Bcd molecules. Increasing the number of Bcd BS from 6 to 9 strongly increases the steepness of the boundary indicating that cooperativity of binding, or more explicitly a longer time to unbind is supported by our model fitting and likely to be at work in this system. Adding 3 more BS to the 9 Bcd BS has very limited impact, indicating that either Bcd molecules bound to the more distal BS may be too far from the TSS to efficiently activate transcription or that the system is saturated with a binding site array occupied with 9 Bcd molecules. In the anterior with excess Bcd, the fraction of time when the loci are active at steady state also increases when adding 3 Bcd BS from B6 to B9. By assuming a model of transcription activation by Bcd proteins bound to target sites, we found that the activation rate increases by much greater fold (~4.5 times) than the number of BS (1.5 to 2 times) suggesting a synergistic effect in transcription activation by Bcd.
The burstiness of the Bcd-only reporters in regions with saturating amounts of Bcd, led us to build a model in two steps. The first step of this model accounts for the binding/unbinding of Bcd molecules to the BS arrays. It is directly related to the positioning and the steepness of the expression boundary and thus to the measurement of positional information. The second step of this model accounts for the dialog between the bound Bcd molecules and the transcription machinery. It is directly related to the fluctuation of the MS2 signals including the number of firing RNAP at a given time (intensity of the signal) and bursting (frequency and length of the signal). Interestingly, while the first step of the process is achieved with an extreme precision (10% EL) (Gregor et al., 2007a;Porcher et al., 2010), the second step reflects the stochastic nature of transcription and is much noisier (150%) Little et al., 2013). Our model therefore also helps to understand and reconcile this apparent contradiction in the Bcd system.
As predicted by our original theoretical model (Lucas et al., 2018), 9 Bcd BS in a synthetic reporter appear sufficient to reproduce experimentally almost entirely the spatial features of the early hb expression pattern i.e. measurements of positional information. This is unexpected as the hb-P2 promoter is supposed to only carry 6 Bcd BS and leaves open the possibility that the number of Bcd BS in the hb promoter might be higher, as suggested previously (Park et al., 2019). In addition, while the B9 reporter recapitulates most spatial features of the hb-P2 reporter, it is twice slower to reach the final boundary position in the cycle than the hb-P2 reporter. This suggested that other maternally provided TFs binding to the hb-P2 promoter contribute to fast dynamics of the hb pattern establishment. Among these TFs, we focused on two known maternal partners of Bcd: Hb which acts in synergy with Bcd Simpson-Brose et al., 1994) and Zld, the major regulator of early zygotic transcription in fruit fly (Liang et al., 2008). Interestingly, adding Zld or Hb sites next to the Bcd BS array reduces the time for the pattern to reach steady state and modify the promoter activity in different ways: binding of Zld facilitates the recruitment of Bcd at low concentration, making transcription more sensitive to Bcd and initiate faster while the binding of Hb affects strongly both the activation/deactivation kinetics of transcription and the RNAP firing rate. Thus, these two partners of Bcd contribute differently to Bcd-dependent transcription. Consistent with an activation process in two steps as proposed in our model, Zld will contribute to the first step favoring the precise and rapid measurements of positional information by Bcd without bringing itself positional information.
Meanwhile, Hb will mostly act through the second step by increasing the level of transcription through a reduction of its burstiness and an increase in the polymerase firing rate. Interestingly, both Hb and Zld binding to the Bcd-dependent promoter allow speeding-up the establishment of the boundary, a property that Bcd alone is not able to achieve. Of note, the hb-P2 and Z2B6 reporters contain the same number of BS for Bcd and Zld but they have also very different boundary positions and mean onset time of transcription following mitosis when Bcd is limiting. This is likely due to the fact that the two Zld BS in the hb-P2 promoter of hb are not fully functional: one of the Zld BS is a weak BS while the other Zld BS has the sequence of a strong BS but is located too close from the TATA Box (5 bp) to provide full activity (Ling et al., 2019).
Zld functions as a pioneer factor by potentiating chromatin accessibility, transcription factor binding and gene expression of the targeted promoter (Foo et al., 2014;Harrison et al., 2011). Zld has recently been shown to bind nucleosomal DNA (McDaniel et al., 2019) and proposed to help establish or maintain cis-regulatory sequences in an open chromatin state ready for transcriptional activation (Eck et al., 2020;Hannon et al., 2017). In addition, Zld is distributed in nuclear hubs or microenvironments of high concentration (Dufourt et al., 2018;Mir et al., 2018). Interestingly, Bcd has been shown to be also distributed in hubs even at low concentration in the posterior of the embryo (Mir et al., 2017). These Bcd hubs are Zld-dependent (Mir et al., 2017) and harbor a high fraction of slow moving Bcd molecules, presumably bound to DNA (Mir et al., 2018). Both properties of Zld, binding to nucleosomal DNA and/or the capacity to form hubs with increased local concentration of TFs can contribute to reducing the time required for the promoter to be occupied by a sufficient number of Bcd molecules for activation. In contrast to Zld, our knowledge on the mechanistic properties of the Hb protein in the transcription activation process is much more elusive. Hb synergizes with Bcd in the early embryo (Simpson-Brose et al., 1994) and the two TF contribute differently to the response with Bcd providing positional and Hb temporal information to the system . Hb also contributes to the determination of neuronal identity later during development (Hirono et al., 2017). Interestingly, Hb is one of the first expressed members of a cascade of temporal TFs essential to determine the temporal identity of embryonic neurons in neural stem cells (neuroblasts) of the ventral nerve cord. In this system, the diversity of neuronal cell-types is determined by the combined activity of TFs specifying the temporal identity of the neuron and spatial patterning TFs, often homeotic proteins, specifying its segmental identity. How spatial and temporal transcription factors mechanistically cooperate for the expression of their target genes in this system is not known. Our work indicates that Hb is not able to activate transcription on its own but that it strongly increases RNAP firing probability and burst length of a locus licensed to be ON. Whether this capacity will be used in the ventral nerve cord and shared with other temporal TFs would be interesting to investigate.
The Bcd-only synthetic reporters also provided an opportunity to scrutinize the effect of Bcd concentration on the positioning of the expression domains. This question has been so far difficult to investigate with endogenous hb and hb-P2 reporters due to the presence of other TF BS in the promoter region which might mediate indirect effects of Bcd or Bcd-independent effects. When comparing the transcription patterns of the B9 reporter in Bcd-2X flies and Bcd-1X flies, we detected a shift of ~10.5 ± 1 %EL of the boundary position. Surprisingly, this shift revealed a gradient of Bcd activity with an exponential decay length of ~15 ± 1.4 %EL (~75 µm), significantly smaller than the value observed directly (~100 µm) with fluorescent-tagging (Gregor et al., 2007b;Liu et al., 2013) or immunofluorescence (Durrieu et al., 2018;Houchmandzadeh et al., 2002) of Bcd proteins. Importantly, reducing the decay length of the gradient also allows our model to fit significantly better the data. We propose that this decay length corresponds to a sub-population of "active" or "effective" Bcd distributed in a much steeper gradient than the Bcd gradient observed using immunofluorescence or GFP tagging which include all Bcd molecules. Bcd molecules have been shown to be heterogenous in intranuclear motility, age and spatial distributions but to date, we do not know which population of Bcd can access the target gene and activate transcription (Tran et al., 2020). The existence of two (or more) Bicoid populations with different mobilities Fradin, 2017;Mir et al., 2018) obviously raises the question of the underlying gradient for each of them. Also, the dense Bcd hubs persist even in the posterior region where the Bcd concentration is low (Mir et al., 2017). As the total Bcd concentration decreases along the AP axis, these hubs accumulate Bcd with increasing proportion in the posterior, resulting in a steeper gradient of free-diffusing Bcd molecules outside the hubs. At last, the gradient of newly translated Bcd was also found to be steeper than the global gradient (Durrieu et al., 2018). Our work calls for further experiments to better characterize the roles of the different Bcd subpopulations. Importantly, reducing by half the Bcd concentration in the embryo induced a similar shift in the position of the hb-P2 reporter boundary as that of the Bcd-only reporters. This further argues that the effective gradient of Bcd activity, with a decay length ~15 %EL, is a principal and direct source of positional information for hb expression.
The steeper effective Bcd gradient found here rekindles the debate on how a steep hb pattern can be formed in the early nuclear cycles. With the previous value of Bcd gradient decay length of =20 %EL, the Hill coefficient inferred from the fraction of loci's active time at steady state is ~6.9, beyond the theoretical limit of the equilibrium model of Bcd interacting with 6 target BS of the hb promoter (Estrada et al., 2016;Hopfield, 1974). This led to hypotheses of energy expenditure in Bcd binding and unbinding to the sites to overcome this limit (Estrada et al., 2016), hb promoters containing more than 6 Bcd sites (Lucas et al., 2018;Park et al., 2019) or additional sources of positional information . The effective decay length ~15% EL, found here with a Bcd-only reporter but also hb-P2, corresponds to a Hill coefficient of ~5.2, just below the physical limit of an equilibrium model of concentration sensing with 6 Bcd BS alone. However, a smaller decay length also means that the effective Bcd concentration decreases faster along the AP axis. In the Berg & Purcell limit (Berg and Purcell, 1977), the time length to achieve the measurement error of 10% at hb-P2 expression boundary with =15 %EL is ~2.1 times longer than with =20 %EL (see Materials & Methods section H where we show the same argument holds regardless of estimated parameter values). This points again to the trade-off between reproducibility and steepness of the hb expression pattern, as described in .  A) Arrangement of the binding sites for Bcd (yellow), Hb (purple) and Zld (blue) upstream of the TATA box (red) and the TSS (broken arrow) of each reporter. B) The MS2 reporters express the iRFP coding sequence followed by the sequence of the 24 MS2 stem loops. In the hb-P2 reporter, the hb-P2 promoter, 5'UTR sequence of the endogenous hb and its intron are placed just upstream of the iRFP sequence. In the synthetic reporters, the minimal promoter of the hsp70 gene was used. Of note, replacing the minimal promoter of hsp70 in B6 by the hb minimal promoter leads to a reporter with lower activity (Fig. S1, F-G). C) Kymographs of mean fraction of active loci (colormap on the right) as a function of time (Y axis in s) and nuclei position along the AP axis (X axis in %EL) at nc13. D) Along the AP axis (%EL), mean time of first spot appearance T0 (s) with shaded standard error of the mean and calculated only for loci with observed expression. E) Cumulative distribution function of T0 (s) in the anterior (-30 ± 2.5 %EL). F) Boundary position (%EL) of fraction of nuclei with MS2 signal along AP axis, with shaded 95% confidence interval, as a function of time. The dash vertical lines represent the time to reach the final decision boundary position (±2 %EL). G) Fraction of nuclei with any MS2 signal, averaged over n embryos, with shaded standard error of the mean, along the AP axis (%EL), at nc13. H) Boundary position and width were extracted by fitting the patterns (fraction of expressing nuclei, G) with a sigmoid function. Bar plots with 95% confidence interval for boundary position and width as the grey region placed symmetrically around the boundary position. Average values and confidence intervals are indicated in the adjacent table. D-H) reporter data are distinguished by color: hb-P2 (orange, n=5 embryos), B6 (yellow, n=5 embryos), B9 (cyan, n=6 embryos) and B12 (green, n=4 embryos).       shaded errors) and Δbcd/+ (Bcd-1X, solid yellow lines with shaded errors) females. Prediction of Bcd-1X pattern from the Bcd-2X pattern assuming a fitted constant shift Δ = 10.5±1.0 %EL for B9 and Δ = 11.0±0.5 %EL for hb-P2 (black dashed line). For B9 (B) data were from n2x = 6 embryos and n1x = 6 embryos and for hb-P2, data were from n2x = 5 embryos and n1x= 3 embryos. C & E) Log-probability map (log (Δ( )) of the shift Δ( ) (in %EL) at a given nuclei position in Bcd-2X embryos ( , in %EL), extracted from the experimental data. The horizontal cyan dashed line represents Δ the best fit value of Δ for B9 (panel C) and for hb-P2 (panel E). The horizontal red dashed line represents the expected shift Δ = -[ ] 2 given [ ] =20 %EL the detected decay length of the Bcd gradient (Gregor et al., 2007b;Houchmandzadeh et al., 2002). F) Comparison of the fitted shift, with 95% confidence interval, in nuclei position from wild-type embryos to embryos from Δbcd/+ females (left bars) and from wild-type embryos to embryos from bcd E1 /+ females (right bars) with B9 and hb-P2 reporters.  ) (solid red) in the steady state window (600s-800s into nc13 interphase) calculated numerically from the best fitted models. Given that the pattern sharpness = / is measured to be 0.34 from the hb-P2 data, the observed Hill coefficient as a function of  is given by the black line with shaded error. The physical limit of equilibrium sensing model with 6 BS ( = 6, black dashed line). B) Log-likelihood of the best fitted models (solid red). The dashed line corresponds to the log-likelihood thresholds for a significantly worse fit (p-value=0.05).
Further information and requests for raw data, resources and reagents should be directed to and will be fulfilled by the lead contact, Nathalie Dostatni (Nathalie.Dostatni@curie.fr).

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Drosophila stocks Embryos were obtained from crosses between females carrying the MCP-NoNLS-eGFP (Garcia et al., 2013) and the His2Av-mRFP (Bloomington # 23561) transgenes both on the second chromosome with males expressing MS2 reporters. Embryos with half dose of Bcd were obtained from females which were in addition heterozygotes for the bcd E1 amorphic allele (bcd 6 ). Unless otherwise specified, all MS2 reporters were inserted at the vk33 docking site (Bloomington # 9750) via φC31 mediated integration system (Venken et al., 2006) by BestGene. The site of insertion was chosen because the transcription dynamics of the original hb-P2 reporter (Lucas et al., 2018) inserted at this site was indistinguishable from the transcription dynamics of two randomly inserted siblings ( Fig. S1A-C). All fly stocks were maintained at 25°C.

MS2 reporters
The hb-P2 MS2 reporter was obtained by cloning the 745bp (300bp upstream of the transcription start site to 445 bp downstream, including the hb intron) located just upstream start codon of Hunchback protein (Drosophila melanogaster) from the previously used hb-MS2ΔZelda reporter (Lucas et al., 2018) into the attB-P[acman]-Cm R -BW plasmid. The synthetic MS2 reporters were created by replacing the hb region in hb-P2 MS2 by the hsp70Bb promoter and the synthetic sequences containing specific combinations of binding sites. GGGATTA was used as a Bcd binding site, CAGGTAG as a Zld binding site and TCAAAAAATAT or TCAAAAAACTAT as Hb binding sites. The sequences of these promoters are in the table S1.

Generation of the bcd mutant by CRISPR
The Δbcd molecular null allele was generated by CRISPR/Cas9 genome editing using the scarless strategy described in (Gratz et al., 2015). gRNA sequences were designed to induce Cas-9 dependent double strand DNA hydrolysis 460 pb upstream of the bcd gene TATA box and 890 bp downstream of the Bcd stop codon. For this, double stranded oligonucleotides (sequences available in Table S1: Oligo for 5' cut Fw, Oligo for 5' cut Rv, Oligo for 3' cut Fw, Oligo for 3' cut Rv) were inserted into the Bbs-1 restriction site of pCFD3-dU6_3gRNA vector (Port et al., 2014). The two homology arms flanking the cleavage sites were amplified from genomic DNA by PCR using the NEB Q5 high fidelity enzyme and specific oligonucleotides (sequences available in Table S1: Bcdnull_5HR_fw, Bcdnull_5HR_rv, Bcdnull_3HR_fw, Bcdnull_3HR_rv). The scarless-DsRed sequence was amplified by with Q5 from the pHD-ScarlessDsRed vector using specific oligonucleotides (sequences available in Table S1: Bcdnull_DsRed_fw, Bcdnull_DsRed_rv). The three PCR amplified fragments were mixed in equimolar ratio with the 2835 bp SapI-AarI fragment of pHD-ScarlessDsRed for Gibson assembly using the NEBuilder system. Injections and recombinant selection based on DsRed expression in the eye were performed by BestGene. For transformants, sequences at the junctions between deletion break point and the inserted dsRed marker were amplified by PCR and verified by sequencing.

METHOD DETAILS
A. Live Embryo imaging Sample preparation and live imaging of transcription was performed as in (Perez-Romero et al., 2018). Briefly, embryos were collected 30 minutes after egg laying, dechorionated by hand and mounted on coverslips covered in heptane-dissolved glue and immersed in 10S Voltatef oil (VWR). All embryos were imaged between nc10 and nc14 at stable temperature (25 o C) on a LSM780 confocal microscope equipped with a 40x (1.4 NA Plan-Apochromat) oil immersion objective. For each embryo, a stack of images is acquired (0.197µm pixel size, 8 bit per pixel, 0.55µs pixel dwell time, confocal pinhole diameter of 92µm, 0.75µm distance between consecutive images in the stack, ~1200x400 pixels image size). GFP and RFP proteins were imaged with 488nm and 561nm lasers, respectively, with appropriate power output. Embryo size and position of the imaged portion is calculated through imaging and measurement of a tiled image of the sagittal plane of the embryo.

B. Data extraction
Data extraction from MS2 movies was performed as in (Lucas et al., 2018) using the LiveFly toolbox (Tran et al., 2018b). In brief, nuclei were segmented in a semi-automatic manner based on His2Av-mRFP channel. The active MS2 loci detection was performed in 3D using a thresholding method. The pixel values at the detected loci location were then fitted with a gaussian kernel to obtain the MS2 loci intensity. The expression data containing each nucleus' position along AP axis and MS2 loci intensity trace over time was exported.
From each time trace in individual nuclei in nc13, we extracted three features: the detection of MS2 expression during the nuclear interphase (taking 0 or 1 values), the time of first MS2 spot appearance and the probability of detecting a spot during the steady state window. This window is defined to be 600-800s into nc13 due to a transient "surge" in transcription activity with the hb-P2 reporter inserted in the vk33 (Fig. S1).

C. Quantifying pattern sharpness, boundary width and Hill coefficient
To quantify the sharpness of the transcription patterns of ms2 reporters, we fit the patterns along the AP axis to a sigmoid function: (1) In Eq. 1, is the position of the nuclei. is the maximum expression level at the anterior pole ( = ( = −∞)). is the expression boundary position ( ( ) = /2). is the scaling coefficient of the AP axis. also corresponds to the pattern sharpness as it is the derivative of the sigmoid function at the boundary position divided by .
From the fitted pattern, we define the boundary width as the distance between two nuclei columns with ( ) of 5% and 95% of maximum expression level .
We assume an exponential Bcd gradient with the decay length ([ ] = / ), where is Bcd concentration at =0. We replace = − log ([ ]/ ) and find the gene expression pattern from the promoter: (3) One should note that in Eq. 3 the Hill function with the Hill coefficient depends on both the pattern sharpness and the decay length of the Bcd gradient : = .
(4) In Eq. 4, the pattern sharpness can be extracted directly from the ms2 movies. Therefore, the assumptions on the decay length will determine the inferred Hill coefficient and consequently the requirements of Bcd binding cooperativity and energy expenditure to achieve such coefficients (Estrada et al., 2016;Tran et al., 2018a).

D. Simulating the model of transcription regulation by Bcd
For a model of transcription regulation with Bcd binding sites, the system can be in 2x +2 states, that consists of ( +1) binding array states ( to ) and 2 transcriptional states ( and ).
We define the state probability vector = [ , , … * ( ) ] T (Σ = 1), where: At = 0, we assume that all the Bcd binding sites are free, and the promoter is OFF. Therefore, we have the initial condition: As all the transitions between promoter and binding array states are first-order reactions, the system dynamics can be described by a 2 +2 by 2 +2 transition matrix with elements .
The binding reaction (from to ) is modeled by: The unbinding reaction (from to ) is modeled by: The activation of the promoter (from OFF state to ON state) is given by: The deactivation of the promoter (from ON state to OFF state) is given by: with ≤ + 1.
(9) To balance the transition matrix, we set the diagonal elements as: (10) With the transition matrix defined, we can calculate the probability vector at a given time : (11) The probability of the promoter to be ON is given by: in which = [ ] is the emission vector specifying the ON state with = 1 with > + 1 and 0 otherwise.

E. Calculating the shift in pattern along AP axis
In this section, we present a framework to quantify the shift of the ms2 expression patterns along the AP axis from Bcd-2x to Bcd-1x flies.
We define two random variables and as the gene expression level in Bcd-2X and Bcd-1X embryos, respectively. Given that the experiments in Bcd-2X and Bcd-1X are independent, the probability distribution of and are independent. and are the nuclei position in Bcd-2X and Bcd-1X respectively, and are described by a uniform distribution from -50 %EL to 50 %EL. At any given position and in the embryo, | ( | = ) and | ( | = ) are assumed to be Gaussian distributed, with the mean and standard deviation equal to the mean and standard error of expression levels (either fraction of expressing nuclei (Fig. 1D).
We call the nuclei position in Bcd-1X that has the same expression level as in the nuclei at position in Bcd-2X. The conditional distribution | ( | = ) is given by: In Eq. 13, ( − ) is a Dirac delta function taking non-zero value only when the expression level in Bcd-2X and Bcd-1X are equal = . As the nuclei position is uniformly distributed from -50 %EL to 50 %EL, ( ) is constant and equal to (% ) .
Eq. 13 can be rewritten as (suppressing the explicit x and s dependence): with the term given by: From Eq. 15 and Eq. 16, we can analytically calculate the probability | ( | ) for each pair of position and from the mean and standard error of gene expression level at a given position.
(17) Given this probability distribution ( ) (Δ) as a function of , we can find a constant value of the shift Δ( ) = Δ that best describes the observed shift for all positions within :

F. Fitting the models of transcription regulation by Bcd
We fit the models of Bcd binding/unbinding to binding sites and activation of transcription, each with a different value of the Bcd gradient decay length , to the transcription dynamics by the synthetic reporters ( Fig. 2 and Fig. S4) and hb-P2 reporter (Fig. 5).

Data
The data used for the fitting is the fraction of active loci ( , , ) at a given time and position along the AP axis (observed from -30 %EL to 20 %EL) in the j th embryo (of total n embryos). ranges from 600 to 800s into nc13 for synthetic reporters and from 0 s to 800 s into nc13 for hb-P2.
We discretize the time axis to increments of 10 s and nuclei position to increments of 1 %EL.

Parameter constraints i) Number of binding sites
The model used in the fitting is described in the Results section and Fig. 2D-E of the main text. The number of Bcd binding sites is = 6 for B6, H6B6, Z2B6 and hb-P2 (Driever and Nusslein-Volhard, 1988), = 9 for B9 and = 12 for B12.

ii) Bcd search rate constant for a target binding site
The Bcd gradient follows an exponential gradient with a predefined decay length : with chosen to be 140 nM or 84 molecules/μm 3 . , which is -35 %EL, is the nuclei position where the Bcd gradient plateaus (Gregor et al., 2007b;Houchmandzadeh et al., 2002;Liu et al., 2013).
Unless otherwise specified, we assume that Bcd molecules independently diffuse in the 3D nuclear space in search for the target sites of the reporters and that the binding to the sites is diffusion limited. Therefore, the binding rate of Bcd to the binding site array depends on the number of free binding sites: (20) Here, is the binding rate constant of individual Bcd molecule to a single target site and equal ( . [ ]) (see Eq. 33 below). The diffusion coefficient is taken to be =4.6 μm 2 /s and the target size =3nm . Therefore, =5.52 μm 3 /s.

iii) Bcd unbinding rate constants
There are no constraints on the unbinding rate constants of Bcd from the target binding array. However, for simplicity, we assume only two forms of Bcd binding cooperativity: bi-cooperativity and 6-cooperativity. It is found from (Tran et al., 2018a) that these two forms of cooperativity are sufficient to generate regulation functions of any order (Hill coefficient between 1 and 6) within the physical limit without extra energy expenditure (Estrada et al., 2016). Therefore, we use three free parameters , and to describe . Assuming at state ( ) , bound Bcd can unbind independently from the BS array, the intermediate binding rates are given by: iv) Promoter switching rates ( ) and value of the fold change in the activation rates. Thus, our conclusions regarding the synergistic activation by bound Bcd molecules still hold. Note that, in the alternative model, where the Bcd-DNA complex regulates the deactivation rate ( ) but not the activation rate . The ratio between and ( ) still holds, as in Eq. 28 and 29. The fold change in the deactivation rates between B6, B9 and B12 is: H. Bicoid search time for its target sites We recalculate the Bcd search time for its target site at the critical hb expression boundary  with different values of the Bcd gradient decay length . The Bcd gradient follows an exponential gradient with a fixed concentration at the anterior : ] is the Bcd concentration at a given nuclei position of distance away from the anterior pole ( ). As Bcd concentration peaks at -35 %EL (Gregor et al., 2007a;Houchmandzadeh et al., 2002), is set to be -35 %EL.
We assume that each Bcd molecule searches for its target sites via 3D diffusion in the nuclear space, with a diffusion coefficient . Each target site has a size of . If the binding of Bcd to the target site is diffusion limited, the Bcd search time for individual target sites at a given position is given by: From Eq. 33, it can be seen that the search time is longer with decreasing .
In practice, there are many uncertainties regarding the exact values of , and , making it difficult to estimate (See Table S4 for details).
If the target gene can sense the Bcd concentration via =6 identical and independent binding sites, we calculate the minimum observation time it takes to achieve the prediction error [ ]/[ ] = 10 % of the Bcd concentration at the hb-P2 boundary ( = -5 %EL) in the Berg and Purcell limit (Berg and Purcell, 1977): We compare the target site search time and the minimum required readout time with equal 12 %EL and 20 %EL in Table S5. We consider different combinatory values of , and found in Table  S4.   Table S1. Promoter sequences of hb-P2, synthetic MS2 reporters and oligonucleotides required for generating the bcd molecular null allele. The sequences are shown in modules, as arranged in Fig. 1 in the main text. The TATA boxes of hb-P2 and HSBG promoter are highlighted in bold. In the binding array sequences, the binding sites for each protein are highlighted in grey (Bcd), green (Hb) and yellow (Zelda). For clarity, the binding array sequences are shown up to the TATA box of HSBG promoter. 19.6 ± 6.9 0.29 ± 0.12 Table S2. Position and width of the gene expression boundary based on fraction of expression nuclei feature ( Fig. 1H and Fig. 3G in Main text) and fraction of active loci feature ( Fig. 2A and Fig. 3H in Main text) in the steady window (600-800s into nc13) for hb-P2 and synthetic reporters in Bcd-2X, shown with 95% confidence interval. For the fraction of expression nuclei, also shown is the time to reach the final activation decision boundary (±2 %EL) starting from the detection the first spot (~225 s) after mitosis.  Table S3. Expected fold change in the activation rates in the anterior region (saturating Bcd concentration) ( ) between B9 and B6 and between B12 and B6. The fold changes are shown for different schemes of activation and values of K (Fig. S2). The fold changes above the value calculated from the data (~4.51) are made bold.  Table S4. Estimated values of Bcd concentration at the anterior ( ), diffusion coefficient ( ) and the size of the target for binding ( ) from previous work.  Figure S1 Figure S1. Comparison of transcription dynamics with hb-P2 reporters inserted at the vk33 site on the 3 rd chromosome and randomly in the 2 nd chromosome and 3 rd chromosome (Lucas et al., 2018). A) Kymographs representing the fraction of nuclei with active ms2 loci (represented as by the colormap) as a function of time and nuclei position along the AP axis. The white dashed horizontal lines represent mitoses between nuclear cycles. B) Fraction of nuclei with any ms2 expression in nc13, averaged over multiple embryos, with shaded error of the mean as a function of nuclei position along AP axis (%EL). C) Bar plots with 95% confidence interval for the expression boundary position for ms2 reporters, based on the fraction of expressing loci (panel B). For each reporter, also shown is the boundary width as the grey region placed symmetrically around the boundary position. D) Time evolution of the fraction of active ms2 loci near the anterior pole (from -35 to -25 %EL) in nc11, nc12 and nc13. E) Time evolution of the ms2 locus intensity per nucleus at the anterior pole (from -35 to -25 %EL) in nc11, nc12 and nc13. In panel B, E and D, data is shown for the hb-P2 reporter inserted at the vk33 site (blue, n=5 embryos), randomly in the 2 nd chromosome (red, n=3) and in the 3 rd chromosome (green, n=6). F) Fraction of nuclei with ms2 expression, averaged over multiple embryos, with shaded standard error of the mean, during nc13 along the AP axis (%EL). G) Fraction of ms2 loci active time at steady state, averaged over multiple embryos, with shaded standard error of the mean, during nc13 along the AP axis (%EL). In panel F and G, data shown for B6 reporters with TATA box of HSBG promoter (solid orange) and of hb-P2 promoter (dashed red). ( ) = . The activation rate with binding state is proportional to the number of bound Bcd: ( )/ ( ) = . C-D) Formation of a stable activating "K-mer": Bound Bcd molecules (orange balls) can randomly form protein complexes (cyan balls) containing molecules (hence called "K-mers"), which can activate transcription at rate . There is no activation of transcription with less than K bound Bcd ( ( ) = 0). In panel C, The "K-mer" is constantly formed and degraded at rates and respectively. Assuming that these rates are fast enough (so that the formation and degradation of "Kmers" are always in equilibrium) and that ≪ , the fold change in the activation rate is given by the number of ways to choose a subset of K molecules from i bound Bcd:    A) Boundary position of fraction of expressing nuclei along AP axis as a function time in nc12 (red) and nc13 (purple), shown for hb-P2 and synthetic reporters, shown with 95% confidence interval. The hb-P2 reporter reach position of expression boundary very rapidly therefore, the position of the boundary is the same at nc12 and nc13 for the hb-P2 reporter even though nc12 is very short. In contrast, positions of the Bcd-only reporters are different at the end of nc12 and nc13. This also true for H6B6 and Z2B6. B) Expression patterns of hb-P2 and synthetic reporters in embryos from wild-type (Bcd-2X, solid green lines with shaded errors) and bcd E1 /+ (solid yellow lines with shaded errors) females. Projection of bcd E1 /+ pattern (black dashed) from the Bcd-2X pattern assuming a fitted constant shift Δ = 9.0±0.5 %EL for hb-P2 (n2x = 5 embryos, nbcdE1/+ = 4 embryos), Δ = 8.5±1.0 %EL for B6 (n2x = 5, nbcdE1/+ = 4), Δ = 8.0±0.5 %EL for B9 (n2x = 6, nbcdE1/+ = 6), Δ = 7.0±0.5 %EL for B12 (n2x = 4, nbcdE1/+ = 6), Δ = 7.0±0.5 %EL for H6B6 (n2x = 7, nbcdE1/+ = 5), Δ = 7.5±1.0 %EL for Z2B6 (n2x = 3, nbcdE1/+ = 5). C) Comparison of the best fitted shift constant from Bcd-2X to bcd E1 /+ flies for hb-P2 and synthetic reporters.