Morphogen gradients infer cell fate as a function of cellular position. Experiments in Drosophila embryos have shown that the Bicoid (Bcd) gradient is precise and exhibits some degree of scaling. We present experimental results on the precision of Bcd target genes for embryos with a single, double or quadruple dose of bicoid demonstrating that precision is highest at mid‐embryo and position dependent, rather than gene dependent. This confirms that the major contribution to precision is achieved already at the Bcd gradient formation. Modeling this dynamic process, we investigate precision for inter‐embryo fluctuations in different parameters affecting gradient formation. Within our modeling framework, the observed precision can only be achieved by a transient Bcd profile. Studying different extensions of our modeling framework reveals that scaling is generally position dependent and decreases toward the posterior pole. Our measurements confirm this trend, indicating almost perfect scaling except for anterior most expression domains, which overcompensate fluctuations in embryo length.
Bicoid (Bcd) is a well studied morphogen involved in the patterning of the anterior–posterior (AP) axis of the Drosophila embryo (Driever and Nüsslein‐Volhard, 1988a, 1988b). Zygotic downstream genes read out this gradient and their expression domains determine the basic body plan of the embryo along this axis. The positions of these domains are remarkably insensitive to fluctuations in the external environment (Houchmandzadeh et al, 2002; Crauk and Dostatni, 2005; Lucchetta et al, 2005) and their relative proportions are maintained across embryos of different sizes. The latter feature is referred to as scaling and occurs within a single species (Houchmandzadeh et al, 2002; Lott et al, 2007) and also across different species (Gregor et al, 2005, 2008; Lott et al, 2007).
Recent experiments also show that the Bcd gradient itself is rather precise and that its length scale correlates to some extent to the embryo size (Gregor et al, 2007a; He et al, 2008). These new findings suggest that the precision and scaling of Bcd target genes may, at least in part, be attributed to that of the morphogen gradient itself. Therefore, we focus here on single morphogen models that aim at explaining the precision and scaling of Bcd target genes at the level of the morphogen gradient. However, it is likely that other mechanisms (e.g. Bcd interactions with the staufen gene (Aegerter‐Wilmsen et al, 2005), gap genes interactions (Jaeger et al, 2004a, 2004b, 2007; Jaeger and Reinitz, 2006; Manu et al, 2009a, 2009b), Bcd interactions with maternal Hunchback and the terminal system (Ochoa‐Espinosa et al, 2009) or bistability (Lopes et al, 2008)) also contribute to further increase robustness in AP patterning.
In this work, we assess the spatial precision of expression domains for the gap genes Krüppel (Kr), Giant (Gt) and Hunchback (Hb), as well as the pair‐rule gene Even‐skipped (Eve) in more than 150 staining images of embryos with a single, double and quadruple dose of bcd. Our data indicate that the precision is maximal at mid‐embryo as well as position dependent rather than gene dependent. This provides independent support for the conclusions drawn from direct measurements of Bcd (Gregor et al, 2007a; He et al, 2008) that the morphogen gradient itself is the main contributor to the precision of the target genes. It motivates our subsequent analytical investigation of noise propagation during Bcd gradient formation within a single morphogen modeling framework. Our analysis shows that fluctuations in the parameters affecting morphogen production, degradation, diffusion as well as nuclear trapping can give rise naturally to highest precision at mid‐embryo, provided the Bcd gradient is decoded at pre‐steady‐state. Finally, we investigate the scaling of the morphogen gradient in a variety of models that couple the embryo size to parameters affecting gradient formation. To single out the most likely scenario, we measure scaling at the level of the target genes. We find that expression domain boundaries posterior to ∼40% embryo length (L) exhibit almost perfect scaling, whereas in the more anterior region we observe over‐compensation to fluctuations in embryo size. This effect appears to be position dependent rather than gene dependent. Interestingly, it is in good agreement with a model that explicitly includes the nuclear trapping of Bcd. We conclude that the formation of the Bcd gradient itself is likely to be a main contributor to robust patterning along the AP axis and that pre‐steady‐state decoding and nuclear trapping are efficient means to increase robustness.
The gap and pair‐rule gene expression domains are precise
We studied staining images of 154 Drosophila melanogaster embryos at cleavage cycle 14 (see Methods in Bergmann et al, 2007). Using image processing tools described in ‘Materials and methods’, we measured the relative positions x/L of the protein domains of Gt, Hb, Kr and Eve. We screened both wild‐type embryos with two copies of bcd and mutant strains with one or four bcd copies resulting in shifted expression domains. We observed that in mid‐embryo, the standard deviation of the relative domain localizations, σ(x/L), is between 1 and 2% and it tends to be higher (2–4%) toward the anterior and posterior poles (Figure 1).
These results are in good agreement with previous results on the Hb domain (Houchmandzadeh et al, 2002). Importantly, they show that the precision of the target genes has the same magnitude and positional trend as that of Bcd according to Gregor et al (2007a). Moreover, our experimental analysis shows that precision is more position dependent than gene dependent. Indeed, when expression domains are shifted because of different bcd mRNA dosages, the precision of the domain seems to change according to the new domain position. Overall, all the investigated domains follow similar precision trends. This suggests that a major contribution to the precision of the Bcd target genes can be attributed to the morphogen gradient itself. We note that position dependent precision may also arise as an experimental artifact when analyzing embryos with different orientations and varying ages (times classes T5–T8 of nuclear cycle 14) or because of imperfect scaling (cf. Supplementary Text S3). However, estimating the maximal contributions of such effects we still find that higher precision at mid‐embryo is statistically significant (cf. Supplementary Text S1; Supplementary Figure S4).
Modeling precision in a single morphogen framework
Given the aforementioned pieces of evidence for achieving precision at the level of the gradient formation, we consider a French‐flag model (Wolpert, 1969) in which a single morphogen gradient induces expression domain boundaries in a concentration‐dependent manner (see Materials and methods). During the first stages of embryonic development, 13 nuclear divisions occur about every 10 min. Gregor et al (2007b) have recently shown that nuclear concentrations stabilize at cycle 10. From a robustness perspective, this means that Bcd read‐out is based on a similar number of molecules and so the associated stochastic noise does not change significantly from cycle 10 on. Therefore, we focus here on the noise propagated from the external Bcd gradient and investigate the fluctuations in domain localization when perturbing the various parameters characterizing the gradient formation (i.e. production, degradation, diffusion and nuclear trapping rates).
We want to assess the variation in domain position Δx, which is induced by embryo‐to‐embryo fluctuations of magnitude Δq in a parameter q affecting the morphogen gradient formation. Δx can be estimated analytically to first order in Δq
where we used that the threshold concentration is fixed, yielding dM=(∂M/∂x)dx+(∂M/∂q)dq=0. This estimate can be computed for any morphogen distribution M(x, t) with explicit dependence on q. Here, we consider the time‐dependent solution of a model in which Bcd is produced at the anterior pole with production rate s0, diffuses according to a uniform diffusion constant D, is degraded at uniform rate α and is trapped and released by the nuclei in the embryo at rates kn, k−n, respectively (see Materials and methods). In Figure 2, we plotted the imprecision measure Δx for small fluctuations (5%) in each of these parameters (Figure 2A–D) as well as a combination of fluctuations in the production and diffusion rates (Figure 2E) (numerical simulations gave very similar results, data not shown). We choose ti=100 min as the time when patterning is initiated, but qualitatively our results only depend on the relation between ti and the Bcd decay time τ=1/α (for ti≪τ the gradient has not yet reached steady‐state). We find that fluctuations in s0 (Figure 2A) give rise to decreasing Δx toward the posterior pole if the pre‐steady‐state gradient is read out (i.e. for small α). In contrast, for fluctuations in D or α, Δx increases toward the posterior pole (Figure 2B and C). Thus, in the case of Bcd pre‐steady‐state decoding (bluish curves) and fluctuations both in the production and degradation or diffusion rates (Figure 2E), maximal precision of the Bcd gradient around mid‐embryo can arise naturally, as was indeed observed directly by Gregor et al (2007a) and indirectly by our analysis (cf. Figure 1). In contrast, for high decay rates (reddish curves) the gradient rapidly equilibrates such that fluctuations in the production rate give rise to uniform noise in Bcd concentrations, preventing the possibility of minimal noise in the central region. Interestingly, fluctuations in the nuclear trapping rate alone also yield higher precision around mid‐embryo (Figure 2D) if the Bcd gradient is decoded early in pre‐steady‐state. However, the position with highest precision is sensitive to the decoding time, shifting more toward the posterior pole for later decoding.
Our modeling approach provides a proof of principle that maximal precision at mid‐embryo can arise if the gradient is decoded before steady‐state. However, we cannot rule out that other mechanisms may yield such a pattern of precision even for a steady‐state Bcd gradient (e.g. in the wing disc, it was shown that cell‐to‐cell variability in the production, diffusion and degradation rates can yield higher precision around mid‐field, Bollenbach et al, 2008).
In our analysis, we assumed the classical French‐flag model (Wolpert, 1969) in which domain boundaries are determined from critical morphogen concentrations. This is a simplification and recent work by the Reinitz group (Jaeger et al, 2004a, 2004b, 2007; Jaeger and Reinitz, 2006; Surkova et al, 2008; Manu et al, 2009a, 2009b) showed that the gap gene expression domains are established by a highly dynamic process characterized by drifting domain boundaries. In our modeling, we neglected the gap gene dynamics and rather focused on how precision can be achieved already at the level of the Bcd gradient, as suggested by our data (as well as those reported in Gregor et al, 2007a; He et al, 2008). Thus, we assume that the positional information is transmitted relatively early to the gap genes, around cycle 10 (as suggested in Lucchetta et al, 2008). The precision and proper scaling of this read‐out is then maintained (and possibly refined) by the gap genes independently of Bcd, because of their cross‐ and auto‐regulation (Bergmann et al, 2007). Meanwhile, the Bcd gradient may continue to evolve and eventually decays at the onset of gastrulation (e.g. by activation of its PEST domain, Niessing et al, 1999).
We next sought to extend our modeling framework to accommodate also scaling. To quantify scaling analytically at the morphogen level, we define a scaling coefficient S(x, t) for a position x and time t
Perfect scaling corresponds to S=1. In this case, fluctuations in embryo length, dL/L, are exactly compensated by fluctuations in position, dx/x, implying perfectly conserved proportions. We use the terms hypo‐ and hyper‐scaling to refer to S<1 and S>1, respectively. A position that hypo‐scales does not compensate enough for a change in embryo size, meaning that in a bigger embryo the absolute position is not shifted enough posteriorward to keep the correct proportions. In contrast, hyper‐scaling is the tendency to overcompensate for a change in embryo size.
Scaling is a property of the external gradient, which is then transmitted to the nuclear concentrations (cf. equation (7)). Assuming that the threshold concentration is fixed (implying dM=(∂M/∂x)dx+(∂M/∂L)dL=0) it follows from equation (2) that the scaling coefficient is
The above definition of scaling is generic and can be computed for any morphogen distribution M(x, t) with explicit dependence on L. Here, we consider a model in which the embryo length impacts the generation of the morphogen gradient through the nuclei density. Specifically, assuming that at decoding time all embryos have the same number of nuclei independent of their size (which agrees with the deterministic doubling of nuclei at each cycle), the nuclei density N depends on the embryo size like N∝L−n, where n∈[1;3] (n=3 corresponds to a uniform distribution of nuclei, whereas n=2 is true if nuclei are distributed on a shell with a fixed width, Gregor et al, 2007b). In this scenario, we find that scaling is time dependent and position dependent (Figure 3B; see Supplementary Figure S9 for the dependence of scaling on N/K and n). Specifically, anterior domain boundaries hyper‐scale, in particular if decoding occurs relatively late, whereas posterior domains show very good scaling for a wide range of decoding times.
We also investigated three other models that give rise to scaling by coupling the nuclear degradation, cytoplasmic degradation or the morphogen production rate to the embryo length (cf. Supplementary Text S2). We find that their scaling behavior is qualitatively similar (in the sense that S decreases toward the posterior pole) but more sensitive to the Bcd decoding time (Supplementary Figure S8).
The gap and pair‐rule gene expression domains scale with embryo size
The model we investigated above predicts that scaling should be position dependent rather than gene dependent and close to perfect except for anterior most domains. To test these two predictions experimentally on our collection of staining images, we adapted the continuous definition of the scaling coefficient S(x, t) in equation (2) for discrete measurements as follows:
where is the estimated slope from a linear regression x=α+βL of the domain positions xi (with mean value ) onto their respective embryo sizes Li (with mean value ) (see Figure 3A). In our data, the relative embryo size fluctuations were of the order of 10–15%. The results of our scaling analysis for Gt, Hb, Kr and Eve in embryos with single, double and quadruple bcd dosage are presented in Figure 3B. We observe that anterior domains indeed tend to hyper‐scale, whereas mid‐embryo and posterior domains show good scaling (P‐value <0.00021, cf. Supplementary Text S1; Supplementary Figure S5). Moreover, the magnitude of the scaling coefficient depends mainly on the position of the respective domain (boundary) rather than the associated gene. This is in good agreement with the model discussed above, whereas in the other models we investigated (cf. Supplementary Text S2) it is difficult to obtain close to perfect scaling both at mid‐embryo and toward the posterior pole.
Precision and scaling are important to achieve robust pattern formation. Our new position‐dependent measures provide a unified quantification of scaling for both functional descriptions of gradient profiles derived from models (equation (2)) and experimentally determined expression domains (equation (4)). These measures will be useful to address scaling also in different systems like the wing‐disk. Importantly, our measures clearly identify perfect scaling (S=1) corresponding to the correct preservation of proportions, which is not necessarily equivalent to perfect correlation (see Supplementary Text S4 for a detailed discussion).
Interestingly, our experimental analysis using embryos with different bcd dosages showed that both precision and scaling are more position dependent than gene dependent, suggesting that the morphogen gradient itself has an important function to set robust domain boundaries. As our study showed, gradient formation of a single morphogen can naturally ensure maximal robustness only in a limited patterning domain. For Bcd, this appears to be mid‐embryo, where most of its targets are expressed. Outside of this domain, other systems and mechanisms may have evolved to cooperate in maintaining and potentially increasing robustness. Considering precision and scaling as position‐ (and possibly time‐) dependent features will allow to better characterize and develop models for these systems.
Materials and methods
We use an interface developed in MATLAB (Supplementary Figure S1) to extract and analyze expression profiles from the staining images of different Bcd target genes, yielding quantitative information on the positions of each protein domain (Bergmann et al, 2007). A total of 154 staining images were analyzed using a semi‐automated analysis tool where the position of the anterior and posterior poles were marked 50 times for each embryo. On the basis of this input, the software extracted a rectangular region from the image from which it generated protein concentration profiles and automatically determined the positional information of the gap and pair‐rule gene expression domains. Thus, for each embryo, its total length and the domain localizations were characterized by a mean value and a standard deviation computed from the 50 markings. The right boundary of the anterior Hb expression domain was defined as the position at which the decline in the Hb concentration is the steepest. Gt and Kr domains were described similarly by their left and right boundaries, whereas the position with maximal staining intensity characterized their center. As the Eve stripes are not so broad, we only extracted the position at which the Eve intensity was maximal. (All the staining images are available for download at http://www.unil.ch/cbg/morphogen/Images_paper.zip)
We also developed a second, almost fully automated analysis tool that only requires human input for selecting which of the two poles is the anterior one (while their position was detected automatically). This tool was used to analyze the expression domain positions of a set of 119 embryos, which included 78 (51%) from the first analysis. Precision and scaling results turned out to be very similar for the two analysis approaches, suggesting that measurement errors and inter‐embryo variation (e.g. because of different embryo orientations) can be safely neglected (Supplementary Figure S6). We also note that part of the observed fluctuations in embryo size may be induced by the fixation procedure. However, computing precision and scaling by considering embryos on different slides separately yielded similar results (Supplementary Figure S7), suggesting that fixation or other batch effects have no positional bias toward either of the embryo poles.
We consider a model in which Bcd is produced at the anterior pole with production rate s0, diffuses according to a uniform diffusion constant D, is degraded at uniform rate α and is trapped and released by the nuclei in the embryo at rates kn, k−n, respectively. The corresponding coupled partial differential equations for the free morphogen concentration M(x, t) and the nuclear concentration Mn(x, t) read (Bergmann et al, 2007; Coppey et al, 2007; Gregor et al, 2007b):
where N is the nuclei density and s0δ(x) the source term localized at x=0 (anterior pole). We consider zero‐flux (Neumann) boundary conditions at the posterior pole and that there is no morphogen at the initial time (t=0). Assuming that nucleo‐cytoplasmic exchanges occur rapidly (∂Mn/∂t∣M≈0), as demonstrated by the experiments reported in Gregor et al (2007b), the following effective diffusion equation holds for the free morphogen (Bergmann et al, 2007; Coppey et al, 2007; Gregor et al, 2007b):
where =D/(1+(N/K)), and (with K=k−n/kn). At time t=150 min (corresponding to nuclear cycle 14), the morphogen profile is numerically fitted to an exponentially decaying gradient with length scale λ=0.2L. Thus, given a degradation rate, the diffusion constant is adjusted accordingly (see caption of Figure 2 for numerical values). We account for the presence of the nuclei by setting N/K=1 (i.e. to be specific, we assume that the probability that external Bcd is trapped equals the probability that nuclear Bcd is released; yet our qualitative results for precision are robust with respect to the exact choice of N/K).
where νn is the nuclear volume. Importantly, according to equation (7) the fluctuations in the nuclear concentrations cn are proportional to those in the external gradient, provided we can neglect variability in νn and K.
We thank Naama Barkai, Ben‐Zion Shilo and Zvi Tamari for careful reading of the manuscript, as well as all the members of the Bergmann Computational Biology Group and David Morton de Lachapelle for their support and comments during the whole project. In particular, we thank Toby Johnson, Karen Kapur and Sascha Dalessi for their help with the theoretical analyses and are very grateful to Zoltán Kutalik who implemented the fully automated version of the profile extraction tool. We are also grateful for financial support from the Giorgi‐Cavaglieri Foundation, the Swiss National Science Foundation (Grant #3100AO‐116323/1), SystemsX.ch (through the WingX project) and the European Framework Project 6 (through the AnEuploidy project).
Conflict of Interest
The authors declare that they have no conflict of interest.
Precision and Scaling in Morphogen Gradient Read‐out
* Image analysis and control experiments:
We discuss the image analysis tools and limitations as well as the validity of our results.
* Scaling models
We discuss other models that could yield scaling of Bicoid as well as the advantages of our measure of scaling compared to correlations.
Dataset S1. Numerical values for the quantitative analysis of precision (relative position x/L, imprecision s(x/L), error, gene name, co‐staining (for Eve domains), bcd dosage).
Dataset S2. Numerical values for the quantitative analysis of scaling (relative position x/L, scaling coefficient, 68% CI, gene name, co‐staining (for Eve domains), bcd dosage).
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2010 EMBO and Macmillan Publishers Limited