Shadow enhancers can suppress input transcription factor noise through distinct regulatory logic
Abstract
Shadow enhancers, groups of seemingly redundant enhancers, are found in a wide range of organisms and are critical for robust developmental patterning. However, their mechanism of action is unknown. We hypothesized that shadow enhancers drive consistent expression levels by buffering upstream noise through a separation of transcription factor (TF) inputs at the individual enhancers. By measuring the transcriptional dynamics of several Kruppel shadow enhancer configurations in live Drosophila embryos, we showed that individual member enhancers act largely independently. We found that TF fluctuations are an appreciable source of noise that the shadow enhancer pair can better buffer than duplicated enhancers. The shadow enhancer pair is also uniquely able to maintain low levels of expression noise across a wide range of temperatures. A stochastic model demonstrated the separation of TF inputs is sufficient to explain these findings. Our results suggest the widespread use of shadow enhancers is partially due to their noise suppressing ability.
eLife digest
In all higher organisms, life begins with a single cell. During the early stages of development, this single cell grows and divides multiple times to develop into the many different kinds of cells that make up an organism. This is a highly regulated process during which cells receive instructions telling them what kind of cell to become.
These instructions are relayed via genes, and a particular combination of activated genes determines the cell’s fate. Specific pieces of DNA, known as enhancers, act as switches that control when and where genes are active, while so-called shadow enhancers are found in groups and work together to turn on the same gene in a similar way.
Shadow enhancers are often active during the early stages of life to direct the formation of specialized cells in different parts of the body. But so far, it has been unclear why it is beneficial to the divide the role of activating genes across several shadow enhancers rather than a single one.
Here, Waymack et al. examined shadow enhancers around a gene called Kruppel in embryos of the fruit fly Drosophila melanogaster. Manipulating the shadow enhancers showed that they help to make gene activity more resistant to changes. Factors such as fluctuations in temperature have different effects on each shadow enhancer. Having several shadow enhancers working together ensures that, whatever happens, the right genes still get activated. For genes like Kruppel, which are key for healthy development, the ability to withstand unexpected changes is a valuable evolutionary benefit.
The study of Waymack et al. reveals why shadow enhancers are involved in the regulation of many genes, which may help to better understand developmental defects. Many conditions caused by such defects are influenced by both genetics and the environment. Genetic illnesses can vary in severity, which may be related to the roles of shadow enhancers. As such, studying shadow enhancers could lead to new approaches for treating genetic diseases.
Introduction
The first evidence that transcription occurred in bursts, as opposed to as a smooth, continuous process, was observed in Drosophila embryos. Electron micrographs showed that even highly transcribed genes had regions of chromatin lacking associated transcripts in between regions densely associated with nascent transcripts (McKnight and Miller, 1979). As visualization techniques have improved, it is increasingly clear that transcriptional bursting is the predominant mode of expression across organisms from bacteria to mammals (Chubb et al., 2006; Dar et al., 2012; Sanchez and Golding, 2013; Zenklusen et al., 2008; Bothma et al., 2014). These bursts of transcriptional activity, separated by periods of relative silence, have important implications for cellular function, as mRNA numbers and fluctuations largely dictate these quantities at the protein level (Csárdi et al., 2015; Hansen et al., 2018). Such fluctuations in regulatory proteins, like TFs and signaling molecules, can propagate down a gene regulatory network, significantly altering the expression levels or noise of downstream target genes (Blake et al., 2003).
Given the inevitable fluctuations in regulatory proteins like TFs, it is unclear how organisms establish and maintain the precise levels of gene expression seen during development, where expression patterns can be reproducible down to half-nuclear distances in Drosophila embryos (Dubuis et al., 2013; Gregor et al., 2007). Many mechanisms that buffer against expression noise, either inherent or stemming from genetic or environmental variation, have been observed (Lagha et al., 2012; Stapel et al., 2017; Raj et al., 2010). For example, organisms use temporal and spatial averaging mechanisms and redundancy in genetic circuits to achieve the precision required for proper development (Stapel et al., 2017; Raj et al., 2010; Erdmann et al., 2009; Lagha et al., 2012). Here, we propose that shadow enhancers may be another mechanism by which developmental systems manage noise (Barolo, 2012).
Shadow enhancers are groups of two or more enhancers that control the same target gene and drive overlapping spatiotemporal expression patterns (Hong et al., 2008; Barolo, 2012). Shadow enhancers are found in a wide range of organisms, from insects to plants to mammals, particularly in association with developmental genes (Cannavò et al., 2016; Osterwalder et al., 2018; Garnett et al., 2012; Bomblies et al., 1999). While seemingly redundant, the individual enhancers of a shadow enhancer group have been shown to be critical for proper gene expression in the face of both environmental and genetic perturbations (Frankel et al., 2010; Osterwalder et al., 2018; Perry et al., 2010). Such perturbations may exacerbate fluctuations in upstream regulators (Cheung and Ma, 2015; Chen et al., 2015). Although shadow enhancers are shown to be pervasive in developmental systems and necessary for robust gene expression, their precise mechanism of action is still unknown. One proposed mechanism is that having multiple enhancers controlling the same promoter reduces the effective ‘failure rate’ of the promoter and ensures a critical threshold of gene expression is reached (Lam et al., 2015; Perry et al., 2011). An alternative, but not mutually exclusive, possibility is that shadow enhancers ensure robust expression by buffering noise in upstream regulators. Several studies suggest that individual enhancers of a shadow enhancer group tend to be controlled by different sets of TFs, which we call a ‘separation of inputs’ (Wunderlich et al., 2015; Cannavò et al., 2016; Ghiasvand et al., 2011). We hypothesize that this separation of inputs allows shadow enhancers to buffer against fluctuations in upstream TF levels to drive more consistent expression levels.
The Drosophila gap gene Kruppel (Kr) provides a useful system in which to probe the mechanisms of action of shadow enhancers. During early embryogenesis, Kr expression is critical for thorax formation, and like the other gap genes in the Drosophila embryo, has quite low noise (Preiss et al., 1985; Dubuis et al., 2013). During this time, Kr is controlled by the activity of two enhancers, proximal and distal (Perry et al., 2011), that drive overlapping expression in the center of the embryo (Figure 1—figure supplement 1). We call the two individual enhancers together the shadow enhancer pair. Previous experiments have shown that each enhancer is activated by different TFs (Figure 1A; Wunderlich et al., 2015). Here, we focus on differences in activation, as the key repressors of Kr, Knirps and Giant, are likely to regulate both enhancers. By measuring live mRNA dynamics, we can use the Kr system in Drosophila embryos to assess whether and how shadow enhancers act to buffer noise and to identify the sources of noise in the developing embryo.
To test our hypothesis, we measured live mRNA dynamics driven by either single Kr enhancer, duplicated enhancers, or the shadow enhancer pair and compared the dynamics and noise associated with each. We showed that the individual Kr enhancers can act largely independently in the same nucleus, while identical enhancers display correlated activity. We constructed a simple mathematical model to describe this system and found that TF fluctuations are necessary to reproduce the correlated activity of identical enhancers in the same nucleus. The shadow enhancer pair drives lower noise than either duplicated enhancer, and using the model, we found that this is a natural consequence of the separation of TF inputs. Additional experiments, including simultaneous measurements of TF levels and expression and a decomposition of noise sources, further demonstrate that the shadow enhancer pair is less sensitive to fluctuations in TF levels than is a single enhancer. Additionally, the shadow enhancer pair is uniquely able to maintain low levels of expression noise across a wide range of temperatures. We suggest that this noise suppression ability is one of the key features that explains the prevalence of shadow enhancers in developmental systems.
Results
Individual enhancers in the shadow enhancer pair act nearly independently within a nucleus
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.
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.
Transcription factor fluctuations are required for the observed differences in the correlations of enhancer activity
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 Ei that interacts with a transcription factor Ti, which together bind to the promoter to form the active promoter-enhancer complex Ci (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 Ti is produced in bursts of ni molecules at a time, and it degrades linearly. For simplicity, the transcription factor Ti is an abstraction of the multiple activating TFs that interact with the enhancer, and Ti 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.
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.
The shadow pair’s activity is less sensitive to fluctuations in Bicoid levels than is the activity of a single enhancer
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.
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−3), 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−2; 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.
The shadow enhancer pair drives less noisy expression than enhancer duplications
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).
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.
Modeling indicates the separation of input TFs is sufficient to explain the low noise driven by the shadow enhancer pair
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 kon and koff 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 kon and koff values, we used the medians of the 10 best single enhancer parameter sets as a starting point and performed simulated annealing to refine them.
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 koff and kon values for each enhancer in the two enhancer constructs to their values from the single enhancer model. We generally found that koff values increased and kon values decreased (Figure 5B). The effect is most pronounced in the duplicated distal enhancer, with large changes to the koff and kon 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 koff and kon 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 koff and kon 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.
The shadow enhancer pair buffers against intrinsic and extrinsic sources of noise
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 noise2, which corresponds to the CV2 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).
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.
The shadow enhancer pair drives low noise at several temperatures
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).
Discussion
Fluctuations in the levels of transcripts and proteins are an unavoidable challenge to precise developmental patterning (Raser and O'Shea, 2005; Arias and Hayward, 2006; Hansen et al., 2018). Given that shadow enhancers are common and necessary for robust gene expression (Osterwalder et al., 2018; Frankel et al., 2010; Perry et al., 2010), we proposed that shadow enhancers may function to buffer the effects of fluctuations in the levels of key developmental TFs. To address this, we have, for the first time, extensively characterized the noise associated with shadow enhancers critical for patterning the early Drosophila embryo. By either tracking biallelic transcription or simultaneously measuring input TF levels and transcription, we tested the hypothesis that shadow enhancers buffer noise through a separation of TF inputs to the individual member enhancers. Our results show that TF fluctuations play a significant role in transcriptional noise and that a shadow enhancer pair is better able to buffer both extrinsic and intrinsic sources of noise than duplicated enhancers. Using a simple mathematical model, we found that fluctuations in TF levels are required to reproduce the observed correlations between reporter activity and that the low noise driven by the shadow enhancer pair may be a natural consequence of the separation of TF inputs to the member enhancers. Lastly, we showed that a shadow enhancer pair is uniquely able to buffer expression noise across a wide range of temperatures. Together, these results support the hypothesis that shadow enhancers buffer input TF noise to drive robust gene expression patterns during development.
Temporal fluctuations in transcription factor levels drive expression noise in the embryo
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.
A stochastic model underscores importance of transcription factor fluctuations
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.
Noise source decomposition suggests competition between reporters
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).
Additional functions of shadow enhancers and outlook
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.
Materials and methods
Generation of transgenic reporter fly lines
Request a detailed protocolThe 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.
Sample preparation and image acquisition
Request a detailed protocolLive 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.
Burst calling and calculation of transcription parameters
Request a detailed protocolTracking 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 (FRNAP, 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).
Simultaneous tracking of Bcd-GFP and enhancer activity
Request a detailed protocolTo 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.
Conversion of integrated fluorescence to mRNA molecules
Request a detailed protocolTo 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, FMS2, and the number of mRNAs generated by a single allele from the same reporter construct in the same time interval, NFISH, 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 (F1) as well as the instantaneous fluorescence of an mRNA molecule (FRNAP). With our microscope, FRNAP is 379 AU/RNAP and F1 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 F1 and FRNAP, respectively.
Calculation of noise metrics
Request a detailed protocolTo calculate the temporal CV each transcriptional spot i, we used the formula:
where mi(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 m1(t) is the fluorescence of one allele in the nucleus at time t, and m2(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.
Statistical methods
Request a detailed protocolTo 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.
Description of the single enhancer model and associated parameters
Request a detailed protocolWe 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 kon 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 koff. In the bursting TFs model, the transcription factor T appears at a rate β1 and degrades at a rate β−1. To recapitulate Kruppel expression patterns, the value of β1 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 kon as described by the reactions
In this case, the value for T was fitted for each bin in a similar way to β1, 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.
Estimation of model parameters from experimental data
Request a detailed protocolTo yield a starting estimate for the kon and koff parameters, we defined the start and end of a burst as the time when the reactions and occur, respectively. The length of the ith burst was defined as the range of [bi,pi] where bi corresponds to the time of the ith instance of the reaction and pi to the time of the ith instance of the reaction . The time between the ith burst and the i + 1th burst is [pi,bi+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 tb be the average time between bursts and td be the average duration of a burst, we can write
and
where N is the number of bursts for spot j, bij and pij 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 tb would imply nearly constant fluorescence intensity instead of bursts.) Finally, the average duration of a burst td can be calculated directly from the data and used to obtain koff by calculating 1/td. Similarly, the average time between bursts tb is readily available from the data giving us kon ≈ 1/tb.
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 ith burst as [pi,bi+1]. For some point t0 ∈ (pi,bi+1), let . Solving this expression for c gives that . Substituting for c in Equation 3 evaluated at t0 results in . 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 t0 ∈ (pi,bi+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 t1 and t2 we can write their respective measured amounts of mRNA as
and
Solving for c1 and c2 gives
The short-term fluctuations of mRNA from R1 to R2 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.
Parameter fitting with simulated annealing
Request a detailed protocolSimulations 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 P0, an initial temperature Γ0, a final temperature Γ’, the number of iterations per temperature N, and a cooling factor µ. Then, each iteration is as follows:
If the current iteration i is such that i > N, then update the current temperature Γk = µkΓ0 to µk+1Γ0 and set i = 0. Otherwise, set i to i + 1.
Check if Γk < Γ’. If so, return the current parameter set Pj and terminate.
Choose a parameter randomly from Pj 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 Pj+1.
Calculate ∆E as the difference in SSE between the data generated by Pj and that generated by Pj+1. Update Pj to Pj+1 if ∆E < 0 or with probability p<e∆E/Γk where p is a uniformly distributed random number.
Repeat all steps until termination.
To generate our results, we chose Γ0 = 1, Γ’ = Γ0/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−3 to 103 for all parameters, except n, which was sampled from 100 to 102, 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 P0. Parameters that were not estimated in the previous section were given the following initial guesses: n = 10, β−1 = 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.
Generation of simulated experimental data
Request a detailed protocolParameter 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.
Description of two enhancer model, parameter estimation, and fitting
Request a detailed protocolTo 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 E1 and E2, 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 E1 to denote the enhancer at the proximal location and E2 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 kon and koff 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 r1 = r2.
To describe the dynamics of the shadow enhancer pair, we denoted the activators for E1 (the proximal enhancer) and E2 (the distal enhancer) by T1 and T2, respectively:
The production rate of T2, γ1, was calculated in the same way as production rate of T1, β1, 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 C1 and C2 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 E1 or E2 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 (β1, β−1, γ1, γ−1, n1, and n2) are not affected by the presence of an additional enhancer. Therefore, in our model, only kon and koff are free to change. To fit the values of kon1, kon2, koff1, and koff2, 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 kon and koff values. We used the resulting values to simulate transcriptional traces and to calculate the predicted CV values shown in Figure 5.
Theoretical modeling of inter-allele noise
Request a detailed protocolTo 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.
Data availability
All data generated or analysed during this study are included in the manuscript and supporting files. Code for analyzing the transcriptional traces and for creating the computational models is available on GitHub: https://github.com/WunderlichLab/KrShadowEnhancerCode (copy archived at https://github.com/elifesciences-publications/KrShadowEnhancerCode).
References
-
Filtering transcriptional noise during development: concepts and mechanismsNature Reviews Genetics 7:34–44.https://doi.org/10.1038/nrg1750
-
Localization of ASH1 mRNA particles in living yeastMolecular Cell 2:437–445.https://doi.org/10.1016/S1097-2765(00)80143-4
-
Redundant Enhancers Mediate Transcriptional Repression of AGAMOUS by APETALA2Developmental Biology 216:260–264.https://doi.org/10.1006/dbio.1999.9504
-
Transcriptional pulsing of a developmental geneCurrent Biology 16:1018–1025.https://doi.org/10.1016/j.cub.2006.03.092
-
Accurate measurements of dynamics and reproducibility in small genetic networksMolecular Systems Biology 9:639.https://doi.org/10.1038/msb.2012.72
-
Stochastic gene expression in a single cellScience 297:1183–1186.https://doi.org/10.1126/science.1070919
-
Role of spatial averaging in the precision of gene expression patternsPhysical Review Letters 103:258101.https://doi.org/10.1103/PhysRevLett.103.258101
-
A general method for numerically simulating the stochastic time evolution of coupled chemical reactionsJournal of Computational Physics 22:403–434.https://doi.org/10.1016/0021-9991(76)90041-3
-
Mapping transgene insertion sites reveals complex interactions between mouse transgenes and neighboring endogenous genesFrontiers in Molecular Neuroscience 11:385.https://doi.org/10.3389/fnmol.2018.00385
-
Mechanisms of transcriptional precision in animal developmentTrends in Genetics 28:409–416.https://doi.org/10.1016/j.tig.2012.03.006
-
Functional sequestration of transcription factor activity by repetitive DNAJournal of Biological Chemistry 282:20868–20876.https://doi.org/10.1074/jbc.M702547200
-
Live imaging of bicoid-dependent transcription in Drosophila embryosCurrent Biology 23:2135–2139.https://doi.org/10.1016/j.cub.2013.08.053
-
Shadow enhancers foster robustness of Drosophila gastrulationCurrent Biology 20:1562–1567.https://doi.org/10.1016/j.cub.2010.07.043
-
Precision of hunchback expression in the Drosophila embryoCurrent Biology 22:2247–2252.https://doi.org/10.1016/j.cub.2012.09.051
-
Effect of promoter architecture on the cell-to-cell variability in gene expressionPLOS Computational Biology 7:e1001100.https://doi.org/10.1371/journal.pcbi.1001100
-
Location effects of a reporter gene on expression levels and on native protein synthesis in Lactococcus lactis and Saccharomyces cerevisiaeApplied and Environmental Microbiology 67:3434–3439.https://doi.org/10.1128/AEM.67.8.3434-3439.2001
-
Regulation by competition: a hidden layer of gene regulatory networkQuantitative Biology 7:110–121.https://doi.org/10.1007/s40484-018-0162-5
-
Single-RNA counting reveals alternative modes of gene expression in yeastNature Structural & Molecular Biology 15:1263–1271.https://doi.org/10.1038/nsmb.1514
Article and author information
Author details
Funding
Eunice Kennedy Shriver National Institute of Child Health and Human Development (R00-HD073191)
- Zeba Wunderlich
Eunice Kennedy Shriver National Institute of Child Health and Human Development (R01-HD095246)
- Zeba Wunderlich
Hellman Foundation
- Zeba Wunderlich
National Institute of Biomedical Imaging and Bioengineering (T32-EB009418)
- Alvaro Fletcher
ARCS Foundation
- Rachel Waymack
National Science Foundation (DMS1763272)
- German Enciso
Simons Foundation (594598)
- German Enciso
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
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.
Copyright
© 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.
Metrics
-
- 3,780
- views
-
- 476
- downloads
-
- 37
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Chromosomes and Gene Expression
- Developmental Biology
Transcription often occurs in bursts as gene promoters switch stochastically between active and inactive states. Enhancers can dictate transcriptional activity in animal development through the modulation of burst frequency, duration, or amplitude. Previous studies observed that different enhancers can achieve a wide range of transcriptional outputs through the same strategies of bursting control. For example, in Berrocal et al., 2020, we showed that despite responding to different transcription factors, all even-skipped enhancers increase transcription by upregulating burst frequency and amplitude while burst duration remains largely constant. These shared bursting strategies suggest that a unified molecular mechanism constraints how enhancers modulate transcriptional output. Alternatively, different enhancers could have converged on the same bursting control strategy because of natural selection favoring one of these particular strategies. To distinguish between these two scenarios, we compared transcriptional bursting between endogenous and ectopic gene expression patterns. Because enhancers act under different regulatory inputs in ectopic patterns, dissimilar bursting control strategies between endogenous and ectopic patterns would suggest that enhancers adapted their bursting strategies to their trans-regulatory environment. Here, we generated ectopic even-skipped transcription patterns in fruit fly embryos and discovered that bursting strategies remain consistent in endogenous and ectopic even-skipped expression. These results provide evidence for a unified molecular mechanism shaping even-skipped bursting strategies and serve as a starting point to uncover the realm of strategies employed by other enhancers.
-
- Chromosomes and Gene Expression
- Developmental Biology
About 70% of human cleavage stage embryos show chromosomal mosaicism, falling to 20% in blastocysts. Chromosomally mosaic human blastocysts can implant and lead to healthy new-borns with normal karyotypes. Studies in mouse embryos and human gastruloids showed that aneuploid cells are eliminated from the epiblast by p53-mediated apoptosis while being tolerated in the trophectoderm. These observations suggest a selective loss of aneuploid cells from human embryos, but the underlying mechanisms are not yet fully understood. Here, we investigated the cellular consequences of aneuploidy in a total of 125 human blastocysts. RNA-sequencing of trophectoderm cells showed activated p53 pathway and apoptosis proportionate to the level of chromosomal imbalance. Immunostaining corroborated that aneuploidy triggers proteotoxic stress, autophagy, p53-signaling, and apoptosis independent from DNA damage. Total cell numbers were lower in aneuploid embryos, due to a decline both in trophectoderm and in epiblast/primitive endoderm cell numbers. While lower cell numbers in trophectoderm may be attributed to apoptosis, aneuploidy impaired the second lineage segregation, particularly primitive endoderm formation. This might be reinforced by retention of NANOG. Our findings might explain why fully aneuploid embryos fail to further develop and we hypothesize that the same mechanisms lead to the removal of aneuploid cells from mosaic embryos.