Canalization of gene expression in the Drosophila blastoderm by gap gene cross regulation
In the previous section, it was shown that the model correctly predicts the positional variance of six borders across a 50% EL region of the embryo. In this section we use the model to characterize the gene interactions responsible for this canalizing behavior of the embryo. This was done in two steps. First, we determined the regulators that set each of the correctly predicted domain borders. The basic idea of this method is that at a border, the relative activation of a gene g(ua) changes from a value close to zero to a value close to one as a result of changes in regulatory input. Those regulatory inputs responsible for changing the state of a gene from off to on can be read from an appropriate graph [33–36] (see Section S3 in Protocol S1 for further details). In the second step, we ascertained how the regulation of a gene at its domain border is affected when subjected to Bcd variation by comparing the strengths of the regulatory inputs at the border position under differing individual Bcd profiles.
The results of the regulatory analysis are shown in Figure 4A–4C for the posterior borders of the hb anterior domain, the central Kr domain, and the posterior kni domain; and in Figure 5A–5C for the anterior border of the posterior kni domain and both borders of the posterior gt domain. This analysis is consistent with earlier results [31,36]. We briefly summarize the results here; see the captions of these figures and previous studies [31,36] for details. In general, these borders are set up by an activator (either bcd or cad) and two repressors. There are two tiers of repression. Two pairs of genes with complementary domains in the modeled region, Kr and gt as well as kni and hb (Figure S4B in Protocol S1), have strong mutual repression. Genes with overlapping domains in the modeled region, gt and kni for example, have weak mutual repression. One border, the anterior border of the posterior gt domain, is set by a single repressor, Kr. This is well supported by the very large anterior shift of this border in Kr – embryos [58].
We studied the effects of varying the Bcd profile in the gap gene circuit as follows. The position of each border was calculated in the 88 simulations using Bcd profiles from individual embryos. With the exception of the borders of the posterior gt domain, the simulations were pooled into 1% EL bins according to the position of each border. We then averaged the Bcd activation and gap repression levels at the border position in each group. To understand how the regulation of a border is affected under Bcd variation, we plotted the pooled Bcd activation and gap repression levels together for each border under consideration (Figures 4D–4I and 5D–5G).
The low positional variation of the posterior border of the anterior hb domain has been studied extensively [9,12,14–16,59]. Our analysis shows that this border forms by Bcd activation and repression from Kr and Kni (Figure 4A). The results of varying the Bcd profile on the regulation of the hb border are presented in Figure 4D and 4G. It shows that both activation and repression levels remain correlated for different positions of the hb border, with embryos having hb borders in more posterior positions having greater activation by Bcd and greater repression by Kr and Kni. This correlation suggests an epistatic balance between activation and repression at the hb border.
These results are in apparent contradiction to a previous report indicating that the standard deviation of the location of the hb border was unchanged in single mutants for either Kr or kni [9]. We found the same result for single mutants (unpublished data) but in double mutants, uncharacterized previously, the standard deviation of the location of the anterior hb border doubles and that of the posterior border of the third anterior gt domain increases significantly (Figure 6). This experimental result strongly supports the model and rules out a picture in which Bcd provides the sole input to hb [16].
The same epistatic balance is seen in the posterior borders of the Kr and kni domains (Figure 4E, 4H, 4F, and 4I). Of the other correctly predicted domain borders, the anterior border of kni is set by repression from Kr and Hb (Figure 5A). Reduced levels of Bcd in the circuit cause this border to form at more anterior positions (Figure 5D); however, the shift is limited by increased levels of Kr and Hb repression toward the anterior (Figure 5E). Both borders of the posterior gt domain show epistatic balance (Figure 5F and 5G), but it does not correlate with position, because Bcd is not a morphogen in this region and does not provide positional information.
We have shown that the reduction in the variation of gene expression in the gap gene system takes place because of cross regulation between zygotic gap genes. Cross regulation is an intrinsic property of the system, and hence we have demonstrated that the system canalizes. More specifically, the experimentally supported analysis of the Results section “Variance Reduction by Gap Gene Cross Regulation and Experimental Verification” shows that the Bcd gradient is not sufficient to generate gap gene borders with the observed low variation.
Our results, together with those of other investigators, rule out three classes of theoretical models that have been invoked to settle certain questions associated with the positional variation of hb. The first two classes of models invoke either active transport [13] or a second gradient at the posterior pole of the embryo [12,14,15] as a mechanism for the reduction in variation of the hb border. With one exception [15], these studies only considered the hb border, while the work reported here correctly predicts the variation of six borders, some of which are not under the control of bcd. Moreover, the border positions in the model reported here scale with egg length (Figure 3E) without the need to invoke an undiscovered gradient. A third type of model [21] invokes transient behavior of the Bcd gradient and its associated thresholds during cycles 10–13 as a mechanism of variance reduction. This model has been invalidated by the experimental demonstration that intranuclear Bcd concentrations are in fact constant during this period of time [44]. Furthermore, threshold concentrations of Bcd propagate toward the posterior of the embryo during the establishment of the Bicoid gradient, but gap domains in the posterior, in fact, shift in an anterior direction [35]. A third contradiction with experiment is that the model used by the authors of [21] produces incorrect mutually exclusive gap domains in the presence of diffusion.
Many seemingly puzzling experimental results concerning the variance of expression borders driven by Bcd can be understood in a simple and unified manner by considering our results together with two recent studies of the Bcd gradient, one in vivo [16] and the other in fixed tissue [18]. Our comparison of Bcd-GFP data with fixed tissue hb data (Figure 1) shows that the range of x Bcd in vivo is 10% EL, but that of x Hb in fixed tissue data is 5% EL. This comparison finds support in the fixed-tissue study [18] in which a modified staining protocol was used and data were acquired in a single imaging cycle to minimize experimental error. The Bcd data thus obtained had the same level of error as Bcd-GFP data [16], yet a direct comparison between the Bcd and Hb data showed that the posterior border of the anterior hb domain was half as variable. The observation that the range of x Bcd is twice that of x Hb is extremely important, because there are many experimental manipulations that can double the range or standard deviation of x Hb, making it equivalent to the range or standard deviation of x Bcd.
First, these results indicate that the increased variance of x Hb in Kr;kni double mutants (Figure 6A caption) can be interpreted as arising from the unmodified action of Bcd without reduction in variance mediated by other genes. This is a key point, because the larger standard deviation of x Bcd observed in data from fixed tissue here and elsewhere [9] appeared to imply that variance reduction still took place in double mutants, albeit at a reduced level. Second, a 9% range of border positions has been observed from reporters driven by a fragment of the hb promoter believed to contain only Bcd binding sites (see Figure 4 in [21]), suggesting that this expression reflects the underlying variation of x Bcd. Third, when chromosome arm 3L is completely removed [9], the standard deviation of x Hb doubles. We suggest that all three of these experimental manipulations uncover unmodified action of Bcd on hb by removing additional modulators respectively in trans, in cis, and by a currently uncharacterized mechanism.
Other important experimental results can be understood by recalling that the maximum accuracy of the Bcd gradient is inversely proportional to the square root of the Bcd concentration [16,19,60], suggesting that borders under unmodified Bcd control will have lower variation at more anterior positions. This fact explains the observation that a synthetic reporter construct driven only by Bcd has the same standard deviation as x Hb [59]. The synthetic construct forms its border at 28.6% EL, where the in vivo data (see Figure 2B in [16]) indicate that the Bcd concentration is about three times greater than at the native hb border. If the standard deviation of 2.2% EL seen here in Kr;kni double mutants is the unmodified response to Bcd at 47% EL, then the expected standard deviation at 29% EL is 1.3%, which is close to the 1.6% observed (see Figure 5C in [59]). Finally, the fact that the standard deviation of the posterior border of the anterior gt domain increases by less than that of hb in Kr;kni mutant embryos can be understood by noting that it is located at about 40% EL. At this position, the in vivo data appear to indicate that the Bcd concentration is 1.69 times greater than that at the hb border. This concentration ratio predicts that the standard deviation of this border in Kr;kni double mutants will be 1.7% EL, while the actual value is 1.9% EL.
In recently published work, He et al. [18] report that the Bcd threshold at which the posterior border of the anterior hb domain forms is not constant but varies in a correlated manner with the slope of the Bcd profile. Figure 4D shows that in our simulations, the Bcd level at which the hb border forms changes as the Bcd profile is varied. Figure 4G shows that this effect is a consequence of the dynamic adjustment of the repression levels, arising from gap gene cross regulation. He et al. speculate that the stability of the Bcd protein might affect its ability to activate hb. However, our results suggest that the dynamic adjustment of the activation level is in fact an indirect consequence of gap gene cross regulation.
One remaining question concerns the reason for the correct prediction of gap gene domain border variation presented here, as it is clear that the fixed tissue Bcd data used for the theoretical study exaggerate the standard deviation of x Bcd [16]. We believe that the high degree of canalization observed (Figure 3A) makes the theoretical model quite insensitive to initial variation, as is true of the biological system itself. In a companion work [61], we show that robust insensitivity to initial conditions is an inherent mathematical property of the gene circuit equations for this system. Moreover, the observed in vivo standard deviation of x Bcd was obtained in an approximately isogenic stock grown at fixed temperature, while natural populations have much more variable genotypes and develop in the presence of temperature fluctuations. These factors alter the profile of the Bcd gradient but have little effect on the border position of downstream targets [9,62]. We believe that the fixed-tissue Bcd profiles used in the present study capture this extra level of variability characteristic of natural populations.
This analysis of canalization has focused on the reduction of variation of border positions of the gap genes with respect to the threshold positions of the maternal determinant Bcd (Results sections “Simulation of Bcd Variation and Size Variation” and “Variance Reduction by Gap Gene Cross Regulation and Experimental Verification”). The gap gene borders are specified at a low enough molecular number of Bcd [16] that intrinsic fluctuations are a significant source of variation [16,60]. The buffering of phenotypic variation against intrinsic fluctuations in molecular number, while not considered by Waddington, is likely to be an important aspect of canalization in biological systems. It is noteworthy that we were able to model canalization using a deterministic model with what were, essentially, stochastic initial conditions. Future work will need to treat the role of intrinsic fluctuations in the formation of maternal gradient systems [63,64]. While stochastic interactions undoubtedly occur among zygotic gap gene products, understanding the functional role of stochastic processes in the segmentation system will require the formulation of a fully stochastic model. Such a model must also produce the observed reduction in variation in the gap gene system, although the numerical details will differ.
Waddington hypothesized that canalization arises from gene interactions [1,65–68]. We have identified specific regulatory interactions among the gap genes responsible for canalization. It is noteworthy that variation reduction is an emergent property of the system and does not require special built-in mechanisms such as a posterior gradient. The mutual repression between the gap genes that is required for correct spatiotemporal expression is also responsible for canalization, providing support for the “intrinsic” [67] or “single-mode” [69] canalization hypothesis.
Our results make it possible to characterize in a more precise way how canalization works. The reduction of gap gene expression variation over time [6] suggests that canalization occurs because the possible developmental trajectories of the system tend to converge over time, thus reducing variation. In an accompanying manuscript [61], we show that the trajectories of the gap gene system in individual nuclei do tend to converge and are stable against perturbations. This convergence and stability are a consequence of the fact that the regulatory interactions of the gene circuit cause the state of the system to approach a stable steady state [70,71], also called an attractor. A given nucleus can potentially approach one of many attractors, a property called multistability, and a particular attractor is chosen based on the initial state of the nucleus. Waddington inferred canalization, in part, from the observation that there are discrete tissue types rather than a continuous blend. Multistability provides a rigorous mathematical justification for the connection between canalization and the discreteness of tissue type. If the initial state is perturbed by a small amount, the system is attracted by the same attractor, providing variation reduction. However, if the perturbation is large enough, or if attractors are created or destroyed, the trajectory can switch to a different attractor providing a discrete response. Our analysis [61] suggests that switching between attractors occurs in the gap gene system and is the mechanism by which the posterior border of the anterior hb domain forms. These results show that the gap gene network is a system with the intrinsic property of reducing variation, a form of error correction, and that this canalization property is a consequence of the mutual regulatory actions of the gap gene network.
The methods used to obtain and characterize the quantitative data are as described in earlier work [6]. All gene expression levels are on a scale of 0–255, chosen to maximize dynamic range without saturation. The Kr;kni double mutant embryos were made by crossing Kr 1 and Df(3L)ri-79c flies. 3–4-hours-old embryos were fixed and stained as described [72] using guinea pig anti-Kr and anti-Kni antibodies with Alexa Fluor 488 conjugated secondary antibody, rabbit anti-Gt and rat anti-Hb antibodies with Alexa Fluor 555 and 647 conjugated secondary antibodies (Molecular Probes), respectively. Embryos null for Kr and Kni showed no signal on the first channel. Double mutants in lateral orientation were scanned using a Leica TCS SP2 confocal microscope and the images were then processed as described [73]. Domain border positions for wild type and mutants were determined by calculating the A–P position at which the relative concentration is half its maximum value in the domain. The local maximum in a domain was determined using a quadratic spline approximation [74].
The numerical implementation of the gene circuit equations is as described [35,36] with the addition of time varying external inputs (see Section S1.2 in Protocol S1). The gap gene circuit was fit to integrated gap gene data [6] using Parallel Lam Simulated Annealing (PLSA) [28,29]. PLSA minimizes the root mean squared (RMS) difference between model output and data. The minimized RMS difference for a circuit is a measure of the quality of fit, and is called the RMS score. The RMS score indicates the average difference between the model and the data. Search spaces, penalty function, and other annealing parameters were as described [27,36]. The circuit analyzed in detail had an RMS score of 10.76, corresponding to a proportional error in expression residuals of about 4–5%. Protocol S1 (Section S1.3) contains further details.
We thank T. Gregor, W. Bialek, D. Tank, and E. Wieschaus for generously sharing the numerical data [16] used to plot Figure 1A. We thank J. P. Gergen, J. Jaeger, K. Krishan, C. Martinez, and A. M. Samsonov for valuable discussions and helpful comments on the manuscript.