Shadow enhancers can suppress input transcription factor noise through distinct regulatory logic

To test our hypothesis that the separation of inputs between Kruppel’s (Kr) shadow enhancers provides them with noise-buffering capabilities, we needed to first test the ability of each enhancer to act independently. In previous work, we found that the individual enhancers in the shadow enhancer pair are controlled by different activating TFs (Wunderlich et al., 2015). These experiments established that the enhancers responded differently to perturbations in key TFs, indicating that each enhancer uses a distinct regulatory logic. The proximal enhancer is activated by Hunchback (Hb) and Stat92E, and the distal enhancer is activated by Bicoid (Bcd) and Zelda (Zld) (Figure 1A). Given this separation of inputs, the shadow enhancer pair could provide a form of noise buffering if variability in gene expression is driven primarily by fluctuations in upstream factors. Conversely, variability in upstream regulators may be low enough in the developing embryo that these fluctuations are not the primary driver of downstream expression noise. If this were the case, the separation of inputs is unlikely to be a key requirement of shadow enhancer function.

To investigate these possibilities, we measured and compared the correlation of allele activity in homozygous or heterozygous embryos that carry two reporter genes. Proximal homozygotes contained the proximal enhancer driving a reporter, inserted in the same location on both homologous chromosomes, and distal homozygotes similarly had the distal enhancer driving reporter expression on both homologous chromosomes (Figure 1B). We also made heterozygous embryos, called shadow heterozygotes, which had one proximal and one distal reporter, again in the same location on both homologous chromosomes. To measure live mRNA dynamics and correlations in allele activity, we used the MS2-MCP reporter system (Figure 1C,D). This system allows the visualization of mRNAs that contain the MS2 RNA sequence, which is bound by an MCP-GFP fusion protein (Bertrand et al., 1998). In the developing embryo, only the site of nascent transcription is visible, as single transcripts are too dim, allowing us to measure the rate of transcription (Garcia et al., 2013; Lucas et al., 2013). In blastoderm-stage embryos with two MS2 reporter genes, we can observe two distinct foci of fluorescence corresponding to the two alleles (Figure 1D; Videos 1, 2, 3, 4, 5, 6), in line with previous results that suggest there are low levels of transvection at this stage (Lim et al., 2018; Fukaya and Levine, 2017). To confirm our ability to distinguish the two alleles, we imaged transcription in embryos hemizygous for our reporter constructs, which only show one spot of fluorescence per nucleus. Our counts of active transcription sites in homozygous embryos correspond well to the expected value calculated from hemizygous embryos (Figure 1—figure supplement 1). Therefore, we are able to measure the correlation of allele activity, although we cannot identify which spot corresponds to which reporter.

Video 1

Download asset

Video 2

Download asset

Video 3

Download asset

Video 4

Download asset

Video 5

Download asset

Video 6

Download asset

We predicted that if variability in gene expression is driven by fluctuations in input TFs, we would observe lower correlations of allele activity in shadow heterozygotes than in either the proximal or distal homozygotes. However, if global factors affecting both enhancers dominate, there would be no difference in allele activity correlations. During the ~1 hr of nuclear cycle 14 (nc14) we found that allele activity is more than twice as correlated in both proximal and distal homozygotes than in shadow heterozygote embryos at 47–57% egg length, which encompasses the central region of Kr expression during this time period (Figure 1). The difference in our ability to measure allele correlation in the more anterior and posterior regions of the embryo stems from the slightly different expression patterns driven by the proximal and distal enhancers (Figure 1—figure supplements 1,2). The lower allele correlation in shadow heterozygote embryos indicates not only that the individual member enhancers of the shadow enhancer pair can act largely independently in the same nucleus, but that differential TF inputs are likely the primary determinants of transcriptional bursts in this system. Notably, heterozygotes still show marginal allele correlation, indicating that some correlation is induced by either shared input TFs or factors that affect transcription globally. The independence of individual Kr enhancers allows for the possibility that shadow enhancers can act to buffer noise by providing distinct inputs to the same gene expression output.

To explore the conditions needed for the two Kr enhancers to act nearly independently within the same nucleus, we generated a simple model of enhancer-driven dynamics. We considered an enhancer E_i that interacts with a transcription factor T_i, which together bind to the promoter to form the active promoter-enhancer complex C_i (Figure 2A). When the promoter is bound by the enhancer, it drives the production of mRNA. Since the MS2 system only allows us to observe mRNA at the site of transcription, we modeled the diffusion of mRNA away from the transcription site as decay. The transcription factor T_i is produced in bursts of n_i molecules at a time, and it degrades linearly. For simplicity, the transcription factor T_i is an abstraction of the multiple activating TFs that interact with the enhancer, and T_i corresponds to a different set of TFs for the proximal and distal enhancer. This nonlinear model generalizes the linear model by Bothma et al., 2015 by explicitly taking into account the presence of TFs.

Figure 2 with 4 supplements see all

Download asset Open asset

We estimated some model parameters directly from experimental data and others by fitting using simulated annealing. The mRNA degradation parameter α and production parameter r were measured directly from fluorescence data without any input from the model (see Materials and methods for details). The remaining parameters were first estimated using mathematical analysis, then fine-tuned using simulated annealing. We found separate parameter sets for the proximal and distal enhancers that, when used to simulate transcription, fit the experimentally measured characteristics of the transcriptional traces, including transcription burst size, frequency, and duration, as well as the total mRNA produced (Figure 2—figure supplement 1).

We hypothesized that a model that lacks fluctuations in the input TFs could not recapitulate the high correlation of transcriptional activity in homozygotes versus the low correlation in heterozygotes. To test this hypothesis, we generated another model of TF production. We call our original model described above bursting TFs. The other model is one in which TF numbers are constant over time, which we call constant TFs and is equivalent to the model in Bothma et al., 2015. If the difference in transcription correlation between homozygotes and heterozygotes is due to fluctuating numbers of TFs, we expected that the bursting TFs model will recapitulate this behavior, while the constant TFs model will not. However, if the constant TFs model is also able to recapitulate the observed difference in correlations, then the correlations are likely a consequence of the identical enhancers simply being regulated by the same set of TFs.

For each model, we used the 10 best parameter sets to simulate transcriptional activity in homozygote and heterozygote embryos and analyzed the resulting allele correlations. We found that the bursting TFs model always produced results in which both homozygote allele correlations are significantly higher than the heterozygote, which qualitatively mirrors the experimental observations (Figure 2B). None of the best fitting parameter sets for the constant TF model were able to produce the experimentally-observed behavior and always resulted in near zero correlations for both the homozygote and heterozygote embryos (Figure 2C). Notably, using the bursting TFs model, all the simulated allele correlations were lower than the experimentally observed values, for example the simulated heterozygote allele correlation was near zero, while the experimental value was 0.14 at the embryo’s midpoint. We hypothesized that this discrepancy was because the model assumes complete independence of the proximal and distal enhancer input TFs, while in reality, there may be some degree of shared inputs, either of known TFs or a general component of the transcriptional machinery. To test this hypothesis, we generated a model that added a common TF to the bursting TFs model and attempted to fit the model parameters. Some of the best parameter sets recapitulated the nonzero correlation of the heterozygote embryos, indicating that a shared factor may play a role in this system; however, this behavior was inconsistent from one parameter set to the next (Figure 2—figure supplement 2). Therefore, we concluded that the simpler bursting TFs model, which consistently recapitulated the key features of the allele correlation data, was more suitable for subsequent analysis.

In conclusion, in our minimalist model of enhancer-driven transcription, the presence of TF fluctuations is required for the observed differences in allele correlation. These results also demonstrate the advantage of using a single generic TF for each enhancer. By abstracting away TF interactions, we reduced the complexity and number of parameters in the model, which allowed us to explore the relationship between TF production and allele correlation.

Both the experimental measurements of allele correlation and the computational model suggest that input TF fluctuations are an appreciable source of noise for enhancer activity. Further, previous experimental work (Wunderlich et al., 2015) and the low correlation of transcriptional activity in heterozygotes (Figure 1E) indicates that each individual Kr enhancer receives different TF input signals. This suggests that the shadow enhancer pair will be less sensitive to an input TF fluctuation than a single enhancer, because the shadow enhancer pair’s activity is dependent on a broader range of TF inputs. To directly observe the relationship between input TF levels and enhancer output, we simultaneously tracked Bcd levels and enhancer activity in individual nuclei (Figure 3; Supp. Videos 7–8). We measured this relationship for both the distal enhancer, which is activated by Bcd, and the shadow enhancer pair, and predicted that the distal enhancer’s transcription dynamics are more strongly influenced by fluctuations in Bcd levels.

Figure 3

Download asset Open asset

Video 7

Download asset

Video 8

Download asset

To allow for tracking of both Bcd levels and enhancer activity, we crossed female flies that express eGFP-tagged Bcd in the place of endogenous Bcd (called Bcd-GFP from here on; Gregor et al., 2007) and MCP-mCherry with male flies homozygous for either the shadow pair or distal enhancer reporter. As the Bcd-GFP transgene was inserted in a Bcd null background, the resulting embryos should receive roughly WT levels of Bcd. The females flies were heterozygous for the maternally deposited Bcd-GFP, and therefore, we estimate that roughly half of the Bcd proteins were labeled. Given the previous work demonstrating the normal function and expression levels of tagged Bcd, we expect the Bcd-GFP levels to be a representative sample of total Bcd (Gregor et al., 2007).

Higher activator TF levels increase enhancer activity. We therefore measured the correlation of nuclear Bcd-GFP levels to the slope of MS2 signal. When the enhancer is active, MS2 signal has a positive slope, when the enhancer is inactive, slope is negative. If the shadow enhancer pair is less sensitive than the distal enhancer to fluctuations in Bcd levels, we would predict higher correlation between Bcd-GFP levels and the activity of the distal enhancer than that of the shadow enhancer pair. We find that the transcription dynamics driven by the distal enhancer are indeed significantly more correlated to nuclear Bcd-GFP levels (median r = 0.18) than the dynamics driven by the shadow pair (median r = 0.14; Figure 3F; p-value=6.1×10⁻³), although both correlations are modest (see Discussion). The lower correlation indicates that transcription driven by the shadow pair is less sensitive to Bcd level fluctuations than is the distal enhancer. Our modeling recapitulates this finding, showing that the separated TF inputs of the shadow pair are sufficient to explain the observed decreased sensitivity to TF fluctuations (median r = 0.14 for the distal enhancer; 0.11 for the shadow pair; p-value=2.2×10⁻²; Figure 3G). These findings indicate that the shadow enhancer pair is better able to buffer fluctuations in a single activating TF than a single enhancer, likely due to the shadow enhancer pair’s separation of TF inputs.

We wanted to further test whether the shadow enhancer pair drives less noisy gene expression output than a simple enhancer duplication. We compared the noise in expression driven by the shadow enhancer pair to that driven by two copies of either the distal or proximal enhancer (Figure 4). If the shadow enhancer pair drives lower noise, this suggests that having two independently acting enhancers is a critical feature of shadow enhancers’ ability to reduce variability and mediate robustness. Alternatively, if duplicated enhancers drive similar levels of expression noise, this suggests that enhancer independence is not critical for shadow enhancer function and that shadow enhancers mediate robustness through a different mechanism, such as ensuring a critical threshold of expression is met (Lam et al., 2015; Perry et al., 2011).

Figure 4 with 1 supplement see all

Download asset Open asset

We tracked transcriptional activity in embryos expressing MS2 under the control of the shadow enhancer pair, a duplicated proximal enhancer, or a duplicated distal enhancer (Figure 4). To measure noise associated with each enhancer, we used these traces to calculate the coefficient of variation (CV) of transcriptional activity across nc14. CV is the standard deviation divided by the mean and provides a unitless measure of noise to allow comparisons among our enhancer constructs. We then grouped these CV values by the embryo position of the transcriptional spots and found the average CV at each position for each enhancer construct. All the enhancer constructs display the lowest expression noise at the embryo position of their peak expression (Figure 4A), in agreement with previous findings of an inverse relationship between mean expression and noise levels (Dar et al., 2016; Figure 4—figure supplement 1). The shadow enhancer pair’s expression noise is ~30% or 15% lower, respectively, than that of the duplicated proximal or distal enhancers in their positions of maximum expression (Figure 4C).

If the primary function of shadow enhancers is only to ensure a critical threshold of expression is reached, we would not expect to also see the lower expression noise associated with the shadow enhancer pair compared to either duplicated enhancer. Furthermore, this decreased expression noise is not simply a consequence of higher expression levels, as the shadow enhancer pair produces less mRNA than the duplicated distal enhancer during nc14 (Figure 4B). The lower expression noise associated with the shadow enhancer pair suggests that it is less susceptible to fluctuations in upstream TFs than multiple identical enhancers.

To explore which factors drive the difference in CVs between the duplicated and shadow enhancer constructs, we extended our model to have a single promoter controlled by two enhancers (Figure 5A). To do so, we assumed that either or both enhancers can be looped to the promoter and drive mRNA production. The rate of mRNA production when both enhancers are looped is the sum of the rates driven by the individual enhancers. We assumed that some parameters, for example the TF production rates and mRNA decay rate, are the same as the single enhancer case. We allowed the parameters describing the promoter-enhancer looping dynamics (the k_on and k_off values) to differ, depending on the enhancer’s position in the construct relative to the promoter and whether another enhancer is present. To fit the k_on and k_off values, we used the medians of the 10 best single enhancer parameter sets as a starting point and performed simulated annealing to refine them.

Figure 5 with 1 supplement see all

Download asset Open asset

This approach allowed us to examine how the model parameters that describe promoter-enhancer looping dynamics change when two enhancers are controlling the same promoter. We compared the k_off and k_on values for each enhancer in the two enhancer constructs to their values from the single enhancer model. We generally found that k_off values increased and k_on values decreased (Figure 5B). The effect is most pronounced in the duplicated distal enhancer, with large changes to the k_off and k_on values for the enhancer in the position far from the promoter (position 2). Given that our model assumes that enhancers act additively and only allows for changes in the k_off and k_on values, these observed effects may indicate that either the presence of a second enhancer interferes with promoter-enhancer looping or that the promoter can be saturated. Our model cannot distinguish between these two possibilities, but these observations are consistent with our (Figure 5—figure supplement 1) and previous results indicating that the Kr enhancers can act sub-additively (Scholes et al., 2019). Additionally, the dramatic changes in k_off and k_on values in the duplicated distal enhancer are consistent with a previous assertion that enhancer sub-additivity is most pronounced in cases of strong enhancers (Bothma et al., 2015).

We used these models to simulate transcription and predict the resulting CVs from the duplicated enhancer and shadow pair constructs. In line with experimental data, we found the model predicts that the shadow pair construct drives lower noise than the duplicated distal or duplicated proximal enhancer constructs in the middle of the embryo (Figure 5C). This is particularly notable, as we did not explicitly fit our model to reproduce the experimentally observed CVs. There is only one fundamental difference between the shadow pair and duplicated enhancer models, namely the use of separate TF inputs for the shadow pair. Therefore, in our simplified model, we can conclude that the separation of input TFs is sufficient to explain the lower noise driven by the shadow enhancer pair construct.

To further understand the sources of noise the shadow enhancer pair is able to buffer, we compared the extrinsic and intrinsic noise associated with the shadow enhancer pair to that associated with either single or duplicated enhancers. To do so, we measured the transcriptional dynamics of embryos with two identical reporters in each nucleus and calculated noise sources using the approach of Elowitz et al., 2002. Intrinsic noise corresponds to sources of noise, such as TF binding and unbinding, that affect each allele separately. It is quantified by the degree to which the activities of the two reporters in a single nucleus differ. Extrinsic noise corresponds to global sources of noise, such as TF levels, that affect both alleles simultaneously. It is measured by the degree to which the activities of the two reporters change together. Intrinsic and extrinsic noise are defined such that, when squared, their sum is equal to total noise², which corresponds to the CV² of the two identical alleles in each nucleus in our system (see Materials and methods). Because our data do not meet one key assumption needed to measure extrinsic and intrinsic noise with the two-reporter approach (see Discussion; Figure 6—figure supplement 1), we use the terms inter-allele noise and covariance in place of intrinsic and extrinsic noise.

Based on our separation of inputs hypothesis and CV data, we expected the total noise associated with the shadow enhancer pair to be lower than that associated with the duplicated enhancers. We predicted that the shadow enhancer pair will mediate lower total expression noise through lower covariance, as the two member enhancers are regulated by different TFs. Given the complexity of predicting inter-allele noise from first principles (Materials and methods; Figure 6—figure supplement 2), we predicted that constructs with two enhancers will have lower inter-allele noise than single enhancer constructs but did not have a strong prediction regarding the relative inter-allele noise among the different two-enhancer constructs. Comparisons of noise between the single and duplicated enhancer constructs would further allow us to discern whether reductions in noise are generally associated with two-enhancer constructs or whether this is a particular feature of the shadow enhancer pair.

Neither the duplicated proximal nor distal enhancers drive significantly lower total noise than the corresponding single enhancers, indicating that the addition of an identical enhancer is not sufficient to reduce expression noise in this system (Figure 6A). The shadow enhancer pair drives lower total expression noise than either single or duplicated enhancer, consistent with the temporal CV data in Figure 4. The median total expression noise associated with the duplicated distal and duplicated proximal enhancers is 1.4 or 2.4 times higher, respectively, than that associated with the shadow enhancer pair (Figure 6A). Note that for measurements of noise, our distal construct places the enhancer at the endogenous spacing from the promoter, as we wanted to control for positional effects on expression and noise (Scholes et al., 2019; Figure 6—figure supplement 3).

Figure 6 with 3 supplements see all

Download asset Open asset

In line with our expectations, the shadow enhancer pair has significantly lower covariance levels than either single or duplicated enhancer (Figure 6B). The shadow enhancer pair also has lower inter-allele noise than all of the other constructs, though these differences are only marginally significant (p=0.13) when compared to the duplicated distal enhancer. Covariance makes a larger contribution to the total noise for the duplicated distal enhancer and the shadow enhancer pair, while inter-allele noise is the larger source of noise for the single distal enhancer and the single or duplicated proximal enhancers (Figure 6B).

The lower total noise and covariance of the shadow enhancer pair support our hypothesis that, by separating regulation of the member enhancers, the shadow enhancer pair can buffer against upstream fluctuations. The lower inter-allele noise associated with the shadow enhancer pair warrants further investigation. A simple theoretical approach predicts that two enhancer constructs will have lower inter-allele noise (Figure 6—figure supplement 2). Given that this is not universally observed in our data, this suggests that there is still much to discover about how inter-allele noise changes as additional enhancers control a gene’s transcription.

We showed the Kr shadow enhancer pair drives expression with lower total noise than either single or duplicated enhancer, yet previous studies have generally found individual member enhancers of a shadow enhancer set are dispensable under ideal conditions (Frankel et al., 2010; Perry et al., 2011; Osterwalder et al., 2018). However, in the face of environmental or genetic stress, the full shadow enhancer group is necessary for proper development (Frankel et al., 2010; Osterwalder et al., 2018; Perry et al., 2011). We therefore decided to investigate whether temperature stress causes significant increases in expression noise and whether the shadow enhancer pair or duplicated enhancers can buffer these potential increases in noise.

Similar to our findings at ambient temperature (26.5°C), the shadow enhancer pair drives lower total noise than all other tested enhancer constructs at 32°C (Figure 7B). At 32°C, the duplicated distal and duplicated proximal enhancers display 35% or 52%, respectively, higher total noise than the shadow enhancer pair. At 17°C, the shadow enhancer pair has approximately 46% lower total noise than either the single or duplicated proximal enhancer, 21% lower total noise than the single distal enhancer, and is not significantly different than the duplicated distal enhancer (Figure 7A). As seen by the variety of shapes in the temperature response curves (Figure 7C), temperature perturbations have enhancer-specific effects, suggesting input TFs may differ in their response to temperature change. The low noise driven by the shadow enhancer pair across conditions is consistent with previous studies showing shadow enhancers are required for robust gene expression at elevated and lowered temperatures (Frankel et al., 2010; Perry et al., 2010).

Figure 7

Download asset Open asset

When measured in fixed embryos, the TFs used in Drosophila embryonic development show remarkably precise expression patterns, displaying errors smaller than the width of a single nucleus (Dubuis et al., 2013; Gregor et al., 2007; Little et al., 2013; He et al., 2008). It therefore was unclear whether fluctuations in these regulators play a significant role in transcriptional noise in the developing embryo. By measuring the temporal dynamics of the individual Kr enhancers, each of which is controlled by different transcriptional activators, we show that TF fluctuations do significantly contribute to the noise in transcriptional output of a single enhancer. Within a nucleus, expression controlled by the two different Kr enhancers is far less correlated than expression driven by two copies of the same enhancer, indicating that TF inputs, as opposed to more global factors, are the primary regulators of transcriptional bursting in this system. Our current findings leave open the possibility that additional mechanisms, such as differences in 3D nuclear organization between different reporters, may also contribute to the differences in noise that we see.

We also showed that activity driven by the Kr shadow enhancer pair is less sensitive to levels of a single TF than is activity driven by an individual Kr enhancer. While prior work has shown that changes in TF levels precede changes in target transcription (Bothma et al., 2018), the sensitivity of individual enhancers to changes in TF levels had not been previously quantified. The correlation between Bcd levels and activity of the distal enhancer is modest, and we expect that this reflects both the influence of additional TF inputs and nuclear heterogeneity that causes the local Bcd levels available to the enhancer to differ from total nuclear levels (Mir et al., 2018). We suspect that the correlation between the activity of the distal enhancer and Bcd levels in the microenvironment surrounding the enhancer is higher than what we were able to measure here. New and emerging technologies will likely allow for live measurements of multiple TF inputs at higher spatial resolution, enabling further insights into the dynamics of expression regulation.

The finding that the Kr shadow enhancer pair is less sensitive to TF levels helps reconcile our finding that the individual Kr enhancers are influenced by fluctuations in input TFs with previous studies showing that endogenous Kr expression patterns are rather reproducible (Little et al., 2013). Previous work has cited the role of spatial and temporal averaging, which buffers noisy nascent transcriptional dynamics to generate more precise expression levels. Shadow enhancers operate upstream of this averaging, driving less noisy nascent transcription than either single enhancers or enhancer duplications.

We developed a stochastic mathematical model of Kr enhancer dynamics and mRNA production that recapitulates our main experimental results. This model is based on that by Bothma et al., 2015, but it is expanded to include the dynamics of a TF that regulates each enhancer. We placed a strong emphasis on the simplicity of this model, for example by using a single abstract TF for each enhancer. This choice both avoids a combinatorial explosion of parameters and makes the model results and parameters easier to interpret. One of the most notable features of the model is that it recreates the differences in noise between shadow and duplicated enhancer constructs without any additional fitting, indicating that these differences in the model system are a direct result of the separation of input TFs to the proximal and distal enhancers.

Future versions of this model can include refinements. For example, in the current model, we do not include the influence of repressive TFs or consider the multiple modes of action used by activating TFs. Future experiments and models can also be designed to identify the mechanism of enhancer non-additivity: changes in promoter-enhancer looping, saturation of the promoter, or other mechanisms.

In our investigation of sources of noise, we decomposed total noise into extrinsic and intrinsic components as in Elowitz et al., 2002. In that study, the authors showed that the activity of one reporter did not inhibit expression of the other reporter, and therefore their calculations assumed no negative covariance between the reporters’ expression output. In our system, we found a small amount of negative covariance between the activity of two alleles in the same nucleus (Figure 6—figure supplement 1). For this reason, we called our measurements covariance and inter-allele noise. The negative covariance we observed indicates that activity at one allele can sometimes interfere with activity at the other allele, suggesting competition for limited amounts of a factor necessary for reporter visualization. The two possible limiting factors are MCP-GFP or an endogenous factor required for transcription. If MCP-GFP were limiting, we would expect to see the highest levels of negative covariance at the center of the embryo, where the highest number of transcripts are produced and bound by MCP-GFP. Since the fraction of nuclei with negative covariance is highest at the edges of the expression domain (Figure 6—figure supplement 1), the limiting resource is likely not MCP-GFP, but instead a spatially-patterned endogenous factor, like a TF.

Currently, the field largely assumes that adding reporters does not appreciably affect expression of other genes. However, sequestering TFs within repetitive regions of DNA can impact gene expression (Liu et al., 2007; Janssen et al., 2000), and a few case studies show that reporters can affect endogenous gene expression (Laboulaye et al., 2018; Thompson and Gasson, 2001). If TF competition is responsible for the observed negative covariance between reporters, a closer examination of the effects of transgenic reporters on the endogenous system is warranted. In addition, TF competition may be a feature, not a bug, of developmental gene expression control, as modeling has indicated that molecular competition can decrease expression noise and correlate expression of multiple targets (Wei et al., 2019).

There are likely several features of shadow enhancers selected by evolution outside of their noise-suppression capabilities. Preger-Ben Noon, et al. showed that all shadow enhancers of shavenbaby, a developmental TF gene in Drosophila, drive expression patterns in tissues and times outside of their previously characterized domains in the larval cuticle (Preger-Ben Noon et al., 2018). This suggests that shadow enhancers, while seemingly redundant at one developmental stage, may play separate, non-redundant roles in other stages or tissues. Additionally, a recent study investigating shadow enhancer pairs associated with genes involved in Drosophila embryonic development found that CRISPR deletions of the individual enhancers result in different phenotypes, suggesting each plays a slightly different role in regulating gene expression (Dunipace et al., 2019). In several other cases, both members of a shadow enhancer pair are required for the precise expression pattern generated by the endogenous locus (El-Sherif and Levine, 2016; Perry et al., 2012; Dunipace et al., 2011; Perry et al., 2011; Yan et al., 2017). These sharpened expression patterns achieved by a shadow enhancer pair may reflect enhancer dominance or other forms of enhancer-enhancer interaction and are likely another important function of shadow enhancers (El-Sherif and Levine, 2016).

In the case of Kr, the endogenous expression pattern is best recapitulated by the shadow enhancer pair, with the individual enhancers driving slightly more anterior or posterior patterns of expression (Figure 1—figure supplement 1; El-Sherif and Levine, 2016). Additionally, the early embryonic Kr enhancers drive observable levels of expression in additional tissues and time points, but these expression patterns overlap those driven by additional, generally stronger, enhancers, suggesting that the primary role of the proximal and distal enhancers is in early embryonic patterning (Hoch et al., 1990). Therefore, while we cannot rule out the possibility that the proximal and distal enhancers perform separate functions at later stages, it seems that their primary function, and evolutionary substrate, is controlling Kr expression pattern and noise levels during early embryonic development.

Here, we have investigated the details of shadow enhancer function for a particular system, and we expect that some key observations may generalize to many sets of shadow enhancers. Shadow enhancers seem to be a general feature of developmental systems (Cannavò et al., 2016; Osterwalder et al., 2018), but the diversity among them has yet to be specifically addressed. While we worked with a pair of shadow enhancers with clearly separated TF activators, shadow enhancers can come in much larger groups and with varying degrees of TF input separation between the individual enhancers (Cannavò et al., 2016; Osterwalder et al., 2018). To discern how expression dynamics and noise driven by shadow enhancers depend on their degree of TF input separation, we are investigating these characteristics in additional sets of shadow enhancers with varying degrees of differential TF regulation. Our current results combined with data gathered from additional shadow enhancers will inform fuller models of how developmental systems ensure precision and robustness.

The single, duplicated, or shadow enhancers were each cloned into the pBphi vector, upstream of the Kruppel promoter, 24 MS2 repeats, and a yellow reporter gene as in Fukaya et al., 2016. We defined the proximal enhancer as chromosome 2R:25224832–25226417, the distal enhancer as chromosome 2R:25222618–25223777, and the promoter as chromosome 2R:25226611–25226951, using the Drosophila melanogaster dm6 release coordinates. The precise sequences for each reporter construct are given in Supplementary file 4. For the allele correlation experiments, each enhancer was cloned 192 bp upstream of the Kr promoter, separated by the endogenous sequence found between the proximal enhancer and the promoter. For transcriptional noise experiments, the distal enhancer was placed at its endogenous spacing, 2835 bp upstream of the promoter, and the proximal enhancer sequence was replaced by a region of the lambda genome that is predicted to have few relevant TF-binding sites. In the shadow enhancer pair or duplicated enhancer constructs, the two enhancers were separated by the sequence separating the proximal and distal enhancers in the endogenous locus.

Using phiC31-mediated integration, each reporter construct was integrated into the same site on the second chromosomes by injection into yw; PBac{y[+]-attP-3B}VK00002 (BDRC stock #9723) embryos by BestGene Inc (Chino Hills, CA). To produce embryos with biallelic expression of the MS2 reporter, female flies expressing RFP-tagged histones and GFP-tagged MCP (yw; His-RFP/Cyo; MCP-GFP/TM3.Sb) were crossed with males containing one of the enhancer-MS2 reporter constructs. Virgin female F1 offspring were then mated with males of the same parental genotype, except in the case of shadow heterozygous flies, which were mated with males containing the other single enhancer-MS2 reporter.

Live embryos were collected prior to nc14, dechorionated, mounted on a permeable membrane, immersed in Halocarbon 27 oil, and put under a glass coverslip as in Garcia et al., 2013. Individual embryos were then imaged on a Nikon A1R point scanning confocal microscope using a 60X/1.4 N.A. oil immersion objective and laser settings of 40uW for 488 nm and 35uW for 561 nm. To track transcription, 21 slice Z-stacks, at 0.5 um steps, were taken throughout the length of nc14 at roughly 30 s intervals. To identify the Z-stack’s position in the embryo, the whole embryo was imaged after the end of nc14 at 20x using the same laser power settings. Later in the analysis, each transcriptional spot’s location is described as falling into one of 42 anterior-posterior (AP) bins, with the first bin at the anterior of the embryo. Unless otherwise indicated, embryos were imaged at ambient temperature, which was on average 26.5°C. To image at other temperatures, embryos were either heated or cooled using the Bioscience Tools (Highland, CA) heating-cooling stage and accompanying water-cooling unit.

Tracking of nuclei and transcriptional spots was done using the image analysis Matlab pipeline described in Garcia et al., 2013. For every spot of transcription imaged, background fluorescence at each time point is estimated as the offset of fitting the 2D maximum projection of the Z-stack image centered around the transcriptional spot to a gaussian curve, using Matlab lsqnonlin. This background estimate is subtracted from the raw spot fluorescence intensity. The resulting fluorescence traces across the time of nc14 are then subject to smoothing by the LOWESS method with a span of 10%. The smoothed traces were used to measure transcriptional parameters and noise. Traces consisting of fewer than three time frames were removed from calculations. To calculate transcription parameters, we used the smoothed traces to determine if the promoter was active or inactive at each time point. A promoter was called active if the slope of its trace (change in fluorescence) between one point and the next was greater than or equal to the instantaneous fluorescence value calculated for one mRNA molecule (F_RNAP, described below). Once called active, the promoter is considered active until the slope of the fluorescence trace becomes less than or equal to the negative instantaneous fluorescence value of one mRNA molecule, at which point it is called inactive until another active point is reached. The instantaneous fluorescence of a single mRNA was chosen as the threshold because we reasoned that an increase in fluorescence greater than or equal to that of a single transcript is indicative of an actively producing promoter, while a decrease in fluorescence greater than that associated with a single transcript indicates transcripts are primarily dissociating from, not being produced from, this locus. Visual inspection of fluorescence traces agreed well with the burst calling produced by this method (Figure 2—figure supplement 3).

Using these traces and promoter states, we measured burst size, frequency and duration. Burst size is defined as the integrated area under the curve of each transcriptional burst, from one ‘ON’ frame to the next ‘ON’ frame, with the value of 0 set to the floor of the background-subtracted florescence trace (Figure 2—figure supplement 3 panel C). Duration is defined as the amount of time occurring between the frame a promoter is determined active and the frame it is next determined inactive (Figure 2—figure supplement 3 panel F). Frequency is defined as the number of bursts occurring in the period of time from the first time the promoter is called active until 50 min into nc14 or the movie ends, whichever is first (Figure 2—figure supplement 3 panel E). The time of first activity was used for frequency calculations because the different enhancer constructs showed different characteristic times to first transcriptional burst during nc14. For these, and all other measurements, we control for the embryo position of the transcription trace by first individually analyzing the trace and then using all the traces in each AP bin (anterior-posterior; the embryo is divided into 41 bins each containing 2.5% of the embryo’s length) to calculate summary statistics of the transcriptional dynamics and noise values at that AP position.

All Matlab codes used for burst calling, noise measurements, and other image processing are available at the Wunderlich Lab GitHub (Waymack, 2020; copy archived at https://github.com/elifesciences-publications/KrShadowEnhancerCode).

To compare the sensitivity of the activity of the shadow pair and distal enhancer to Bcd levels, we tracked the fluorescence of Bcd-GFP and MCP-mCherry in individual nuclei across the time of nc14. To obtain embryos for simultaneous tracking, we crossed female flies heterozygous for Bcd-GFP and MCP-mcherry with male flies homozygous for either the shadow pair or distal enhancer reporter. Bcd-GFP and MCP-mCherry are maternally deposited and thereby allow us to measure levels of Bcd and enhancer activity in individual nuclei of the resulting embryos. Embryo collection and preparation was performed as described above. The same microscope, objective, and Z-step profile were used as described above, but laser settings were switched to 40uW for 561 nm and 35uW for 488 nm. Analysis of transcriptional activity was performed as described above. Time traces of Bcd-GFP levels in individual nuclei were subjected to background correction by subtracting the average fluorescence of the regions of the image not containing a nucleus at each time point from the raw Bcd-GFP fluorescence. The resulting Bcd-GFP time traces were then subjected to smoothing by the MATLAB smooth function, using the LOWESS method with a span of 10%. To measure the sensitivity of enhancer activity to Bcd levels, we correlated the slope of MS2 traces to the corresponding Bcd-GFP levels in the same nucleus. Slope was calculated between the MS2 values at consecutive time points and compared to the Bcd-GFP value at the earlier of the two time points. This process was done for all time points through 50 min into nc14.

To put our results in physiologically relevant units, we calibrated our fluorescence measurements in terms of mRNA molecules. As in Lammers et al., 2018, for our microscope, we determined a calibration factor, α, between our MS2 signal integrated over nc13, F_MS2, and the number of mRNAs generated by a single allele from the same reporter construct in the same time interval, N_FISH, using the hunchback P2 enhancer reporter construct (Garcia et al., 2013). Using this conversion factor, we can calculate the integrated fluorescence of a single mRNA (F₁) as well as the instantaneous fluorescence of an mRNA molecule (F_RNAP). With our microscope, F_RNAP is 379 AU/RNAP and F₁ is 1338 AU/RNAP∙min. With these values, we are able to convert both integrated and instantaneous fluorescence into total mRNAs produced and number of nascent mRNAs present at a single time point, by dividing by F₁ and F_RNAP, respectively.

To calculate the temporal CV each transcriptional spot i, we used the formula:

where m_i(t) is the fluorescence of spot i at time t.

We also decomposed the total noise experienced in each nucleus to inter-allele noise and co-variance, analogous to the approach of Elowitz et al., 2002.

Inter-allele noise is calculated one nucleus at a time. It is the mean square difference between the fluorescence of the two alleles in a single nucleus:

where m₁(t) is the fluorescence of one allele in the nucleus at time t, and m₂(t) is the fluorescence of the other allele in the same nucleus and the angled brackets indicate the mean across the time of nc14.

Covariance is the covariance of the activity of the two alleles in the same nucleus across the time of nc14:

The inter-allele and covariance values are defined such that they sum to give the total transcriptional noise displayed by the two alleles in a single nucleus.

This total noise value is equal to the coefficient of variation of the expression of the two alleles in a single nucleus across the time of nc14.

To determine any significant differences in total noise, covariance, or inter-allele noise values between the different enhancer constructs, we performed Kruskal-Wallis tests with the Bonferroni multiple comparison correction.

We constructed a model of enhancer-driven transcription based on the following chemical reaction network:

where E is an enhancer that interacts with a transcription factor T, which together bind to the promoter at a rate k_on to form the active promoter-enhancer complex C. When the promoter is in this active form, it leads to the production of mRNA denoted by R, which degrades by diffusion from the gene locus at a rate α. Transcription is interrupted whenever the complex C disassociates spontaneously at a rate k_off. In the bursting TFs model, the transcription factor T appears at a rate β₁ and degrades at a rate β₋₁. To recapitulate Kruppel expression patterns, the value of β₁ was assumed to be given by

where x is the percentage along the length of the egg and c is a scaling constant. Since Kruppel activity peaks near the center of the egg, we chose µ = 50, while c and σ were fitted along with the other parameters. Lastly, n was assumed to be fixed across the length of the egg.

We also generated a constant TF model, which is an adaptation of the model in Bothma et al., 2015. This model implicitly assumes that TF numbers are constant and, therefore, are incorporated into the value of k_on as described by the reactions

In this case, the value for T was fitted for each bin in a similar way to β₁, that is the constant number of TFs was assumed to be described by Equation (1) (values were rounded to the nearest integer).

To simulate the transcriptional traces, we implemented a stochastic approach. Individual chemical events such as enhancer-promoter looping take place at random times and are influenced by transcription factor numbers. Individual trajectories of chemical species over time were calculated using the Gillespie algorithm (Gillespie, 1976), and these trajectories are comparable to the experimentally measured transcriptional traces. Since the enhancer is either bound or not bound to the promoter, we imposed the constraint that C + E = 1 when simulating model dynamics.

To yield a starting estimate for the k_on and k_off parameters, we defined the start and end of a burst as the time when the reactions ${\displaystyle E\to C}$ and ${\displaystyle C\to E}$ occur, respectively. The length of the ith burst was defined as the range of [b_i,p_i] where b_i corresponds to the time of the i^th instance of the reaction ${\displaystyle E\to C}$ and p_i to the time of the i^th instance of the reaction ${\displaystyle C\to E}$. The time between the i^th burst and the i + 1^th burst is [p_i,b_i+1]. The Gillespie algorithm dictates that the time spent in any given state is determined by an exponentially distributed random variable with a rate parameter equal to the product of two parts: the sum of rate constants of the outgoing reactions, and the number of possible reactions. If the enhancer is either bound or unbound, we have that C = 1 or E = 1, respectively. Therefore, by letting t_b be the average time between bursts and t_d be the average duration of a burst, we can write

and

where N is the number of bursts for spot j, b_ij and p_ij denote the beginning and end of burst i in spot j respectively, and M denotes the total number of spots in the egg. The right-hand sides are given by the expected value of the exponential distribution and the assumption that, on average, T is close to 1. While this may not be the case for T, the assumption provides a convenient upper bound for the average time between bursts, which is likely not to have a much smaller value for a lower bound. (A low enough value of t_b would imply nearly constant fluorescence intensity instead of bursts.) Finally, the average duration of a burst t_d can be calculated directly from the data and used to obtain k_off by calculating 1/t_d. Similarly, the average time between bursts t_b is readily available from the data giving us k_on ≈ 1/t_b.

We were able to directly estimate mRNA production and degradation rates from the experimental data. To estimate α, we focused on periods of mRNA decay; that is periods where no active transcription is taking place and are thus described by

which in turn can be solved to be

where c is a constant of integration. Taking the derivative of Equation 2 yields

which corresponds to the slope of the decaying burst. We define the interval of decay of the i^th burst as [p_i,b_i+1]. For some point t₀ ∈ (p_i,b_i+1), let ${\displaystyle R_0=R(t_0)=c{e}^{-t_0\alpha }}$. Solving this expression for c gives that ${\displaystyle c=R_0{e}^{t_0\alpha }}$. Substituting for c in Equation 3 evaluated at t₀ results in ${\displaystyle R^{\prime }(t_0)=-\alpha R_0{e}^{t_0\alpha }{e}^{-t_0\alpha }=-\alpha R_0}$. Then, it follows that

In other words, the rate of decay of mRNA fluorescence can be calculated from any trace by taking the ratio of the slope during burst decay and its intensity at a given time t₀ ∈ (p_i,b_i+1).

Adjacent measurements of fluorescence intensity from the single enhancer systems were used to approximate the slope at each point in the traces. Then, Equation 4 was applied to each point. A histogram of all calculated values was generated (Figure 2—figure supplement 4). In this figure, there was a clear peak, which provided us with an estimate of α ≈ 1.95.

The estimation of r was done for periods of active transcription, which are also accompanied by simultaneous mRNA decay. By noting that C = 1 during mRNA transcription, we can approximate these periods as the zeroth order process

The differential equation associated with this system is given by

and has steady state R^∗=r/α. Equation 5 can be solved explicitly for R by choosing

where c is a constant of integration. For two adjacent measurements at times t₁ and t₂ we can write their respective measured amounts of mRNA as

and

Solving for c₁ and c₂ gives

The short-term fluctuations of mRNA from R₁ to R₂ between two adjacent discrete time points in the stochastic system can be approximated by Equations 6 and 7. This implies that

which in turn gives

Therefore, the estimation of r can be computed given two adjacent measurements of fluorescence and the time between them. Finally, we used a similar approach as done with α to calculate values of r from fluorescence data. However, unlike α, r was calculated for each bin to account for differences in transcriptional efficiency across the length of the embryo.

Simulations and parameter fitting were done with MATLAB. Optimization in fitting was done by minimizing the sum of squared errors (SSE) between the normalized vectors of burst properties and allele correlations of the experimental and simulated data. In particular, a vector y of experimental data was created by concatenating the following vectors: burst size, integrated fluorescence, frequency, duration, and allele correlation across the length of the embryo. The vector y was subsequently normalized by dividing each burst property by the largest element in their respective vectors (except correlation which by definition is unitless between −1 and 1). A vector x was created in an analogous fashion to y but using simulated instead of experimental data. However, x was normalized using the same elements that were used to normalize y. Then, the discrepancy between the experimental and simulated data was measured by:

We used a high-performance computing cluster to compute 200 independent runs of parameter fitting with simulated annealing for each model variant. The algorithm requires an initial guess of the parameter set P₀, an initial temperature Γ₀, a final temperature Γ’, the number of iterations per temperature N, and a cooling factor µ. Then, each iteration is as follows:

1. If the current iteration i is such that i > N, then update the current temperature Γ_k = µ^kΓ₀ to µ^k+1Γ₀ and set i = 0. Otherwise, set i to i + 1.

2. Check if Γ_k < Γ’. If so, return the current parameter set P_j and terminate.

3. Choose a parameter randomly from P_j and multiply it by a value sampled from a normal distribution with a mean equal to 1. The standard deviation of such distribution should be continuously updated to be Γ_k. The result of this step is the newly generated parameter set P_j+1.

4. Calculate ∆E as the difference in SSE between the data generated by P_j and that generated by P_j+1. Update P_j to P_j+1 if ∆E < 0 or with probability p<e^∆E/Γk where p is a uniformly distributed random number.

5. Repeat all steps until termination.

To generate our results, we chose Γ₀ = 1, Γ’ = Γ₀/10, N = 30, and µ = 0.8. We observed an improvement in the quality of the fittings by using analysis-derived parameter values as initial guesses instead of values given through random sampling. The sampled space ranged from 10⁻³ to 10³ for all parameters, except n, which was sampled from 10⁰ to 10², and σ, which was randomly chosen to be an integer between 1 and 20. Equal numbers of parameter values were sampled at each order of magnitude. The analysis in the section above was used to estimate the parameters in P₀. Parameters that were not estimated in the previous section were given the following initial guesses: n = 10, β₋₁ = 1, σ = 6, and c = 40. Initial guesses for c and σ were based on the experimental observation that there is little transcription outside of 20–80% egg length. Based on this observation, simulations were limited to this egg length range, as well. For the constant TFs model, both analysis-derived and random initial parameter values were used to maximize the likelihood of finding any parameter set capable of recapitulating the observed allele correlation.

Parameter sets resulting from fitting were sorted in ascending order based on their sum of squared errors, and the 10 lowest error parameter sets are what we called the 10 best parameter sets. For all figures, we simulated 80 spots per bin and simulated each bin five times to generate error bars. Data for the distal enhancer at the proximal location was used to reproduce simulated allele correlations in all cases.

Gillespie simulations update the counts of each chemical species at random time intervals. However, for ease of parameter fitting and to better recapitulate the experiments, we generated data in two distinct timescales: one consisting of 30 second intervals after which mRNA counts were recorded, and another consisting of random time intervals generated by the algorithm after which chemical counts were updated. The former one was used for all parameter fitting rounds and generation of figures.

To explore two enhancer systems, we expanded our previous model to include an additional enhancer. First, we considered duplicated enhancer systems, which consist of either two proximal or two distal enhancers. Enhancers were denoted by E₁ and E₂, which correspond to two identical enhancers that exist in different locations relative to the promoter. They are activated by the same transcription factors as described by the reactions

Without loss of generality, we used E₁ to denote the enhancer at the proximal location and E₂ to denote the enhancer at the distal location. This model describes independent enhancer dynamics; that is the behavior of one enhancer does not affect the behavior of the other, and, as such, both enhancers can be simultaneously looped to the promoter. Consequently, to account for potential enhancer interference or competition for the promoter, we assumed distinct k_on and k_off values for each enhancer in the duplicated enhancer constructs. We also used distinct values of r for each distal enhancer in the duplicated distal construct since fluorescence data was available for this enhancer at the proximal and endogenous location. For proximal enhancers, we assume r₁ = r₂.

To describe the dynamics of the shadow enhancer pair, we denoted the activators for E₁ (the proximal enhancer) and E₂ (the distal enhancer) by T₁ and T₂, respectively:

The production rate of T₂, γ₁, was calculated in the same way as production rate of T₁, β₁, but differed in the values of c and σ. The two enhancer models were also used to calculate allele correlation between homozygotes and heterozygotes because a distinction between the mRNA produced by C₁ and C₂ was made. This approach works because, for example when considering the homozygote embryos, each single enhancer resides in the same nucleus and is therefore affected by the same fluctuating TF numbers. In the duplicated enhancer model, each enhancer E₁ or E₂ is affected by the same fluctuations in the number of transcription factor T. An analogous logic applies to the heterozygotes.

To fit the two enhancer models to experimental data, we retained several parameters from the single enhancer models. Parameters r and α were directly calculated from the data, and, as such, did not vary across models. We assume that parameters concerning transcription factors (β₁, β₋₁, γ₁, γ₋₁, n₁, and n₂) are not affected by the presence of an additional enhancer. Therefore, in our model, only k_on and k_off are free to change. To fit the values of k_on1, k_on2, k_off1, and k_off2, we set the other model parameters to the median values of the 10 best parameter sets in the respective single enhancer model. We then used a similar simulating annealing approach to fit the k_on and k_off values. We used the resulting values to simulate transcriptional traces and to calculate the predicted CV values shown in Figure 5.

To make a prediction about the expected change in inter-allele noise between single and two enhancer reporter constructs, we used the theory put forth in Sánchez and Kondev, 2008; Sanchez et al., 2011. This formalism can be used to calculate the expected mean and variance of the transcriptional output of a promoter, given the possible states of the promoter, transition rates between states, and the rate of transcription resulting from each state. In these papers, the authors apply their formalism to different promoter architectures. Here, we generate a simpler model, in which we abstract away the individual transcription factor (TF) binding configurations, which would be numerous and poorly parametrized, and simply define states by whether an enhancer is looped to the promoter and activating transcription. Since these models do not account for fluctuations that would contribute to extrinsic noise, for example fluctuations in TF or RNA polymerase levels, they can predict the dependence of intrinsic noise on enhancer arrangement.

To apply this model to our system, we used these parameters: γ, degradation rate of mRNA; p, production rate of mRNA; k, on rate for enhancer-promoter looping; l, off rate for enhancer-promoter looping. We generated five models that represent different configurations of either one or two enhancers controlling a single promoter. Key assumptions are that the parameters describing this system are independent of both the position of the enhancer relative to the promoter and the presence of a second enhancer controlling the same promoter. In Model 1, there is a single enhancer controlling one promoter. There are two states, when the enhancer and promoter are not looped (mRNA production rate of 0), and when the enhancer and promoter are looped (mRNA production rate of p). In Model 2 (OR model), there are two enhancers controlling one promoter, transcription is activated if either enhancer is looped, and both enhancers can’t be bound at the same time. In Model 3 (additive model), there are two enhancers controlling one promoter, transcription is activated if either enhancer is looped, and, if both enhancers are bound, transcription occurs at twice the rate of single enhancer looping states. In Model 4 (synergistic model), there are two enhancers controlling one promoter, transcription is activated if either enhancer is looped, and, if both enhancers are bound, transcription occurs at three times the rate of single enhancer looping states. In Model 5 (XOR model), there are two enhancers controlling one promoter, transcription is activated if either enhancer is looped, and, if both enhancers are bound, no transcription occurs. Results of these models are shown in Figure 6—figure supplement 2.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

The authors thank Hernan Garcia for providing advice on MS2 imaging, the MCP-GFP fly line, and a plasmid containing the MS2 reporter construct, Nick Lammers for assistance in the calibration of MS2 signal to mRNA counts, Clarissa Scholes for assistance in image processing, Lily Li for assistance in data analysis and comments on the text, Ceazar Nave for feedback on the text and figures, and Thomas Schilling and Rahul Warrior for useful suggestions for the project. We thank the reviewers and editor for their helpful suggestions during the review process.

© 2020, Waymack et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.