Canalization of gene expression in the Drosophila blastoderm by gap gene cross regulation

Developing embryos exhibit a robust capability to reduce phenotypic variations that occur naturally or as a result of experimental manipulation. This reduction in variation occurs by an epigenetic mechanism called canalization, a phenomenon which has resisted understanding because of a lack of necessary molecular data and of appropriate gene regulation models. In recent years, quantitative gene expression data have become available for the segment determination process in the Drosophila blastoderm, revealing a specific instance of canalization. These data show that the variation of the zygotic segmentation gene expression patterns is markedly reduced compared to earlier levels by the time gastrulation begins, and this variation is significantly lower than the variation of the maternal protein gradient Bicoid. We used a predictive dynamical model of gene regulation to study the effect of Bicoid variation on the downstream gap genes. The model correctly predicts the reduced variation of the gap gene expression patterns and allows the characterization of the canalizing mechanism. We show that the canalization is the result of specific regulatory interactions among the zygotic gap genes. We demonstrate the validity of this explanation by showing that variation is increased in embryos mutant for two gap genes, Krüppel and knirps, disproving competing proposals that canalization is due to an undiscovered morphogen, or that it does not take place at all. In an accompanying article in PLoS Computational Biology (doi:10.1371/journal.pcbi.1000303), we show that cross regulation between the gap genes causes their expression to approach dynamical attractors, reducing initial variation and providing a robust output.

These results demonstrate that the Bicoid gradient is not sufficient to produce gap gene borders having the low variance observed, and instead this low variance is generated by gap gene cross regulation. More generally, we show that the complex multigenic phenomenon of canalization can be understood at a quantitative and predictive level by the application of a precise dynamical model.

C. H. Waddington inferred developmental canalization of gene expression by noting that differentiation leads to the formation of discrete types of tissue, rather than a continuous blend, and that genotypic and environmental variation is suppressed at the phenotypic level in wild-type but not mutant organisms. These points led him to state that “developmental reactions, as they occur in organisms submitted to natural selection, are in general canalized. That is to say, they are adjusted so as to bring about one definite end-result regardless of minor variations in conditions during the course of the reaction” (p. 563 in [1]). These adjustments manifest themselves as a reduction in the variation of phenotypes. In this article, we make a detailed experimental and theoretical study of the canalization process in the gap gene system of Drosophila. Interestingly, this investigation puts us in contact with a number of recent investigations of variation in the gap gene system as we now explain.

Among contemporary geneticists and evolutionary biologists, the buffering of phenotypic variation to underlying genotypic variation in wild type is well known, and recent experimental studies have identified individual genes that are responsible for this aspect of canalization [2,3]. Theoretical studies have demonstrated that this phenotypic buffering is an intrinsic consequence of the adjustments of developmental trajectories postulated by Waddington [4], but direct evidence of developmental canalization was still lacking. Direct evidence of developmental canalization must meet two requirements. First, the system must show variation between individuals that decreases significantly over the course of development. Second, this decrease in variation must be demonstrated to be an inherent property of the system under study, rather than being imposed from a different part or from outside the organism. The first requirement implies the need for data on the dynamics of developmental determinants, some or all of which will be molecular. The demonstration that reduction in variance of a developmental system is inherent to that system requires precise understanding of the consequences of the interaction of many developmental determinants. Understanding such consequences requires a quantitative model.

The segmentation system of Drosophila melanogaster [5,6] is highly suitable for studies of canalization because it is already known that the first of the above requirements, decrease in phenotypic variability over time, is satisfied. In work reported elsewhere [6], we performed a quantitative analysis of the expression of segmentation genes expressed from the maternal genome only (bicoid), the maternal and zygotic genomes (caudal and hunchback), and the zygotic genome (Krüppel, giant, knirps, tailless, fushi tarazu, even skipped, runt, hairy, odd skipped, paired, and sloppy paired). These expression data have cellular resolution in space and about 6.5-min resolution in time, and comprise either the expression levels of even skipped (eve) and two other genes in individual embryos or the integrated averaged expression of all 14 protein gene products. These data show that segmentation gene expression is highly variable among individual embryos in cleavage cycle 13 and the early part of cycle 14. There is extensive variation in expression levels, locations of domain borders, and the time and order of the appearance of individual domains. However, the variation in the expression patterns reduces over time and is significantly lower at the onset of gastrulation than at earlier times.

In this article, we restrict our attention to a particularly important class of phenomena concerning the variation in the location of the boundaries of zygotic gene expression domains. These boundaries shift when maternal gradients are perturbed, and hence are at least in part under the control of these gradients [7,8]. Nevertheless, the variation in the boundary positions of gap and pair-rule expression domains is much lower than the variation in the maternal gradient of Bicoid (Bcd) protein [6,9,10]. Under the simplest model of specification, expression borders would form at a fixed threshold of its concentration [11]. The variation of the Bcd gradient can be measured in terms of the range (ρBcd) or standard deviation (σBcd) of the position (x Bcd) where it crosses the threshold concentration. In our data σBcd = 4.6% egg length (EL), whereas gap gene domain border positions have σ ∼ 1% EL [6]. Pair-rule stripe border positions also have similar low variation [6] compared to the Bcd threshold position.

Of the large number of gene expression borders that have lesser variation than Bcd, the posterior border of the anterior hunchback (hb) domain has received intense scrutiny [9,12–18]. Gregor and colleagues [16] measured the absolute concentration of Bcd in the nuclei at the hb border using a Bcd-GFP fusion rescue construct. This measurement led the authors to two important but mutually contradictory conclusions. First, using a result from bacterial chemotaxis theory [19], it was shown that noise due to small number of Bcd molecules is high enough so that the hb border cannot be accurately specified by Bcd alone. Second, by measuring Bcd profiles from several live embryos in parallel, it was found that σBcd ∼ 2% EL. The amount of variance reduction implied by this value of σBcd is significantly smaller than the estimates from fixed-tissue experiments [6,9], presumably because the GFP measurements do not have scaling error introduced in the setting of microscope gain in separate fixed-tissue experiments. Under the assumption that Bcd is the only regulator of hb, it was argued that these results can be reconciled if spatial averaging from the diffusion of Hb molecules is taken into account, and the authors draw the conclusion that Bcd alone is sufficient to specify the hb border accurately.

The assumption that hb is under the sole control of Bcd is incorrect, because its border position changes in embryos mutant for giant (gt) [9], Krüppel (Kr) [20], Krüppel;knirps (Kr;kni) [20], and in embryos lacking Nanos [9]. Apart from the average position, the variation of the position is also increased in embryos lacking chromosome arm 3L or Nanos to ∼2% EL and ∼1.6% EL respectively [9], close to the Bcd variation measured in live embryos. The Bcd-dependent response of a fragment of the hb promoter has variation close to Bcd variation (see Figure 4 in [21]) (ρBcd ∼ 9% EL; compare with Figure 1A) and higher than the variation of endogenous hb (Figure 1B). All these experiments suggest that endogenous hb expression has lower variation than Bcd, which is increased to the level of Bcd variation when other inputs to hb are removed.

To make the relationship between Bcd variation and hb border variation clear, we have plotted the data of Gregor et al. and compare it with fixed-tissue Hb data (Figure 1). The hb border is steep, and its position can be measured unambiguously in fixed tissue since it is not sensitive to scaling variation, unlike a Bcd threshold concentration. This comparison shows that the Bcd threshold position at which the hb border forms is about twice as variable as the border position. It is important to note that the hb border is just one among many that have low positional variance. As an example, consider the anterior border of the posterior kni domain, which is activated by Bcd [22]. This border is located ∼10% EL posterior to the hb border and has a positional variance of ∼1% EL [6]. At this position, the Bcd-GFP data have a variance of ∼4% EL (see Figure 5C in [16]), which is much larger than that of the kni border. This example clearly contradicts a picture of positional specification in which Bcd specifies all of its targets accurately.

Like hb, kni is not regulated by Bcd alone, but also by Caudal (Cad) [22,23], Hb [20,24], and Gt [25]. The examples of the hb and kni borders show that it is not possible to explain the variance of a single expression border position in isolation. It is necessary to take into account all the genetic interactions required for generating the expression patterns over space and time. Gene circuits [26,27] are dynamical models that can accurately reproduce observed gene expression patterns by reconstituting the required set of genetic interactions in silico. In the gene circuit method, (1) a gene circuit model, in which one real number characterizes the regulatory effect of one gene on another, is (2) fit to quantitative gene expression data by (3) parallel Lam simulated annealing [28,29] or some other nonlinear optimization method [30–32], and finally (4) biological conclusions are obtained. This method has been successfully used to analyze both the pair-rule and gap systems [27,33–36] and performed better than other models in a comparative study [31].

In what follows, we briefly describe the gene circuit used in the Results section “Gap Gene Circuits”; a full description is provided in Protocol S1 (Section S1). We then demonstrate in the Results section “Simulation of Bcd Variation and Size Variation” that a gene circuit model for the gap genes hb, Kr, gt, and kni under maternal control and including mutual gap–gap interactions is sufficient to explain the reduction in the variance of gap gene domain border positions, showing that it is not necessary to postulate an undiscovered gradient [12,14,15] or active transport [13] to explain this property. In the Results section “Variance Reduction by Gap Gene Cross Regulation and Experimental Verification,” we characterize the regulatory interactions underlying variance reduction. This analysis shows that the reduction in variance of the hb border is a consequence of repression by Kr and Kni, and hence predicts an increase in variance when these factors are removed. We present new experimental data from Kr;kni double mutant embryos that validate this prediction and show the correctness of the gene circuit approach. The validation of the gene circuit approach shows that the gap gene system also satisfies the second requirement for demonstrating canalization, which is that variance reduction arises from intrinsic properties of the system.

We constructed a gene circuit that models the expression patterns of the gap genes hb, Kr, gt, and kni in a region of the blastoderm extending from 35% EL to 92% EL along the anteroposterior (A–P) axis. The protein products of these genes are transcription factors [37–40] and localize to nuclei. The model computes the time course of nuclear protein concentrations in a time interval beginning at the start of cleavage cycle 13 and ending at the onset of gastrulation [41] (Figure S2 in Protocol S1). The initial conditions were specified using cleavage cycle 12 gene expression data [6] (Section S1 in Protocol S1). The 13th nuclear division occurs in the duration of the model (Figure S2 in Protocol S1) after which the nuclei in the modeled region are divided, and daughter nuclei inherit the state of the mother nucleus.

The gene circuit equations model the processes of protein synthesis, degradation, and diffusion (Section S1.1 in Protocol S1). During the mitosis preceding the 13th nuclear division, the protein concentrations are governed only by diffusion and degradation as synthesis shuts down [42]. During interphase, the synthesis rate of a protein depends on the regulation of the gene by the other gap proteins, the maternal gradient Bcd, and the upstream regulators Cad and Tailless (Tll). The regulatory effect of a protein is characterized by the product of its concentration and its regulatory strength. The concentration varies from nucleus to nucleus, but the regulatory strength, a parameter of the model, is invariant throughout the embryo, reflecting the fact that the zygotic genome is the same in each blastoderm nucleus. The sum of all the regulatory terms for a gene, called the total regulatory input (u), determines its synthesis rate through the sigmoidal regulation-expression function (Figure S1 in Protocol S1).

To fully specify the model the concentrations of Bcd, Cad and Tll must be provided. The concentrations of Cad and Tll were provided by the interpolation of average data [6] in time (Section S1.2 in Protocol S1). A representative profile of Bcd was provided as follows.

The Bcd gradient is essentially stationary during cleavage cycles 13 and 14 [6,43,44], and hence, its concentration v Bcd is assumed to be constant in time. It is known through antibody studies [6,9,12,43] and a recent GFP-Bcd [44] study that the Bcd profile is an exponential function of A–P position x, so that

The arithmetic mean of exponential curves is not exponential. Thus, it is not possible to generate a representative Bcd profile by taking an arithmetic average over embryos. Instead of averaging, a representative individual Bcd profile v Bcd(x) was chosen in the following manner. We obtained data from 88 cycle 13 embryos immunostained for Bcd and removed the background signal from the Bcd profiles as described [45]. Taking the logarithm of Equation 1, we get

The background-removed profiles were then fit by linear least-squares to Equation 2. This procedure yielded two parameters for each profile, γ (length scale), and A (concentration scale). Figure 2A is a scatter plot of log A with γ. A profile (“median” profile) was chosen such that its parameters lie in the middle of the scatter plot (Figure 2A). Since this study concerns only embryo-to-embryo variation and not nucleus-to-nucleus noise, the exponential fit of the profile was used in the model. Using background-removed profiles directly in the model yields circuits with the same properties as those of circuits with exponential fits (See section “Simulation of Bcd Variation and Size Variation” below and Table S4 in Protocol S1).

To compare with quantitative gene expression data [6], we calculate the solution of the model at multiple time points. The data are classified into time classes to give nine time points for comparison (see Table 1 and Figure S2 in Protocol S1).

We obtained specific gene circuits by performing a least-squares fit of the model to gene expression data from nine time points (see Methods). We used a median Bcd profile from an individual embryo as described in the previous section, together with averaged expression profiles from the other genes considered. Twenty-three circuits were obtained with almost identical network topology (Table S3 in Protocol S1) and expression patterns. One of these circuits was chosen for further study. Its parameters are given in Protocol S1 (Tables S1 and S2). The circuit's gap gene patterns (Figure S4A and S4B in Protocol S1) are consistent with data except for two minor defects (Section S1.3 in Protocol S1) and the circuit has the same network topology that was obtained previously [31,35,36].

We simulated the effects of embryo-to-embryo variation by running an ensemble of 88 simulations of the chosen circuit such that each simulation used a Bcd profile from a distinct individual embryo (Figure 2B). All parameters of the circuit were kept fixed, and only the Bcd input to the circuit was varied across the ensemble. The positional variance of the gap gene borders was compared with that of Bcd profiles. The position of a gap gene domain's border was calculated as the point where the concentration of the gene's protein product was half its maximum level in that domain (see Methods). The positional variance of a border is the standard deviation of its position in the 88 simulations. The positional variance of the family of Bcd profiles is the standard deviation (σBcd) of the position (x Bcd) at which each profile crosses a threshold concentration. The Bcd threshold for a particular gap gene border is the Bcd concentration at which the border forms in the circuit with the median Bcd.

The model correctly predicts the positional variance of six gap gene borders: the posterior border of the anterior hb domain, the posterior border of the central Kr domain, both borders of the posterior kni domain, and both borders of the posterior gt domain (see Figure 3A and 3C and Table 2). For example, the posterior border of the anterior hb domain has a standard deviation of 1.3% EL in the simulations and 1.0% EL in data (Figure 3B and 3D), both of which are much less than the standard deviation of the Bcd threshold position, 4.6% EL (Figure 2B). To further confirm that the gene circuit was reducing Bcd variation, we performed the in silico control experiment of turning off gap gene cross regulation of hb in the model. In this simulation, the variation of the posterior border of the anterior hb domain increased to 4.3% EL, close to the level of variation in the Bcd threshold position (Figure 3F).

Of the 23 circuits that had consistent network topology, 15 showed low positional variation for these six borders. The ability of the gene circuit to filter out Bcd's positional variation does not depend on the particular Bcd profile used to fit the circuit. Circuits produced with other profiles from the middle of the parameter scatter (black circles in Figure 2A) also showed accuracy when subjected to Bcd variation (see Table S4 in Protocol S1). The gene circuit has this filtration property irrespective of whether smoothed or unsmoothed profiles are used to simulate Bcd variation (Table S4 in Protocol S1).

Four borders show large variation: the posterior border of the third gt domain, the anterior border of the central Kr domain, and both the borders of the posterior hb domain (see Figure 3A). All these borders are at the edges of the modeled region, where the model is missing key regulators. In the anterior, the head gap genes orthodenticle, empty spiracles, and buttonhead [46–48] are excluded. They are known to influence the anterior pattern of trunk gap genes [25], and respond to Bcd in a concentration-dependent manner [49]. In the posterior terminal region huckebein (hkb) [50–52] has been omitted from the model.

Drosophila eggs vary in length by 10–20%, and different wild lines may have differing mean egg lengths [53]. Although the location of eve stripes [53] and the posterior border of the anterior hb domain scale with egg length, the Bcd gradient and hence its threshold positions do not [9]. To explain this scaling property, an as-yet-undiscovered gradient has been proposed [12,14]. These proposals specifically considered the scaling of the hb border, but one theoretical study suggests that the second gradient could scale multiple borders [15]. The hypothetical gradient would be set up from the posterior pole of the embryo and would directly regulate the gap genes or Bcd. The proposed posterior gradient would provide the target with a measure of distance from the posterior pole.

We show here that gap gene circuits exhibit scaling with egg length. In other words, it is not necessary to postulate an undiscovered posterior gradient. bcd mRNA is localized at the anterior tip of the embryo [54]. The stationary exponential profile of the Bcd gradient is consistent with the diffusion of the protein from the anterior tip and its degradation by a first-order process, and the gradient has been modeled as such [21,55,56]. Recently, an alternative model, in which degradation is replaced with the reversible trapping of the Bcd protein by the nuclei in the blastoderm, was proposed [57]. Both the trapping model and the diffusion-degradation model do not show scaling of the Bcd gradient with egg length and so we represented the scaling properties of the Bcd gradient using the simpler diffusion-degradation model.

Let L denote the length of an embryo and let z represent A–P position so that z varies between 0 and L. Using Equation 1, the concentration of Bcd is given by v Bcd(z) = A exp(–γz), where γ is the characteristic length scale of the profile and A is the concentration scale. Let x denote the relative position in % EL so that . The Bcd concentration can then be written in terms of x as

where γrel(L) = γL is the length scale in relative position units. Equation 3 may also be arrived at by explicitly rescaling the diffusion-degradation model [21,55,56].

Size variation can be simulated by varying γrel(L) in Equation 3. For an embryo with length L′ = L + ΔL

Using the parameters A and γrel(L) of the median Bcd profile (Figure 2B) we varied between −0.1 and 0.1 to generate Bcd profiles representing egg length variation of 20%, which corresponds to the range of egg lengths observed in data [9].

We simulated the gap gene expression patterns with these Bcd profiles and the chosen gene circuit (Section S1.3 in Protocol S1). The results of these simulations are shown in Figure 3E. The six gap gene border positions that have reduced variation in the simulations of Bcd variation also scale proportionally with egg length.