Systems biology derived source-sink mechanism of BMP gradient formation
- University of Pennsylvania Perelman School of Medicine, United States
- Purdue University, United States
Abstract
A morphogen gradient of Bone Morphogenetic Protein (BMP) signaling patterns the dorsoventral embryonic axis of vertebrates and invertebrates. The prevailing view in vertebrates for BMP gradient formation is through a counter-gradient of BMP antagonists, often along with ligand shuttling to generate peak signaling levels. To delineate the mechanism in zebrafish, we precisely quantified the BMP activity gradient in wild-type and mutant embryos and combined these data with a mathematical model-based computational screen to test hypotheses for gradient formation. Our analysis ruled out a BMP shuttling mechanism and a bmp transcriptionally-informed gradient mechanism. Surprisingly, rather than supporting a counter-gradient mechanism, our analyses support a fourth model, a source-sink mechanism, which relies on a restricted BMP antagonist distribution acting as a sink that drives BMP flux dorsally and gradient formation. We measured Bmp2 diffusion and found that it supports the source-sink model, suggesting a new mechanism to shape BMP gradients during development.
https://doi.org/10.7554/eLife.22199.001eLife digest
Before an animal is born, a protein called BMP plays a key role in establishing the difference between the front and the back of the animal. Cells nearer the front of the embryo contain higher amounts of the BMP protein, whilst cells nearer the back have progressively lower levels of BMP. This gradient of BMP ‘concentration’ affects the identity of the cells, with the level of BMP in each cell dictating what parts of the body are made where.
The prevailing view among scientists is that the BMP gradient is created by an opposing gradient of another protein called Chordin, which is found at high levels at the back of the embryo and lower levels near the front. Chordin inhibits BMP and the interaction between the two proteins establishes the gradients that create order across the embryo.
Zinski et al. used computer models to investigate how the BMP gradient is created. Several possibilities were considered, including the effect of Chordin. Comparing the models to precise experimental measurements of BMP activity in zebrafish embryos suggested that a different mechanism known as a source-sink model, rather than the opposing Chordin gradient, may be responsible for the pattern of BMP found in the embryo. In this model, the BMP is produced at the front of the embryo and moves towards the back end by diffusion. At the back of the embryo, BMP is mopped up by Chordin, resulting in a constant gradient of BMP along the embryo.
Many other processes that control how animals grow and develop rely on the formation of similar protein gradients, so these findings may also apply to other aspects of animal development. Understanding how animals grow and develop may help researchers to develop strategies to regrow tissues and organs in human patients.
https://doi.org/10.7554/eLife.22199.002Introduction
Morphogen gradients pattern axonal pathways, the neural tube, the dorsal-ventral (DV) and anterior-posterior (AP) embryonic axes, as well as multiple organ systems (Bökel and Brand, 2013; Briscoe and Small, 2015; Cohen et al., 2013; Rogers and Schier, 2011; Rushlow and Shvartsman, 2012; Sansom and Livesey, 2009; Schilling et al., 2012; Tuazon and Mullins, 2015). Morphogens are defined as factors that form a spatially non-uniform distribution spanning multiple cell-lengths that instructs different cell fates at distinct levels. Their importance in specifying multiple cell fates in a gradient has spurred decades of research deciphering how they work. In 1970, Francis Crick proposed that such a gradient could be formed by a source of morphogen flowing to a sink that destroyed it (Crick, 1970). We now know that the mechanisms by which morphogen gradients are established are diverse and complex, and that understanding these mechanisms is paramount to understanding developmental biology (Briscoe and Small, 2015; Müller et al., 2013; Rogers and Schier, 2011). Bone Morphogenetic Proteins (BMPs) act as morphogens repeatedly during development, including in patterning the embryonic DV axis, the neural tube, and the Drosophila wing disc (Bier and De Robertis, 2015; Briscoe and Small, 2015; Rogers and Schier, 2011).
BMP morphogen systems are established by a network of extracellular regulators (Dutko and Mullins, 2011). A crucial class of these regulators is the BMP antagonists, defined by their ability to bind BMP ligand with high affinity, thereby blocking ligand-receptor interaction (Brazil et al., 2015). During axial patterning in zebrafish and Xenopus, three antagonists, Chordin, Noggin, and Follistatin play key roles in inhibiting BMP signaling to promote dorsal cell fate specification (Dal-Pra et al., 2006; Khokha et al., 2005; Schulte-Merker et al., 1997). These antagonists bind to the BMP ligand, preventing BMP from binding its receptors (Iemura et al., 1998; Piccolo et al., 1996; Zimmerman et al., 1996). In both zebrafish and frogs, Chordin differs from Noggin and Follistatin in its expression domain, its phenotype, and its interaction with the metalloprotease Tolloid (Tld). Chordin is expressed in a larger domain than Noggin and Follistatin (Dal-Pra et al., 2006; Khokha et al., 2005). While the loss of Chordin ventralizes the embryo (Hammerschmidt et al., 1996; Oelgeschläger et al., 2003), the depletion of Noggin or Follistatin alone or together does not (Dal-Pra et al., 2006; Khokha et al., 2005). Unlike Noggin and Follistatin, Chordin can be cleaved by the metalloprotease Tolloid, releasing bound BMP ligand and allowing it to signal (Blader et al., 1997; Piccolo et al., 1997).
Previous studies in Drosophila show that the Drosophila ortholog of Chordin, Sog, can act as both a BMP agonist and as an antagonist during DV patterning. To act as an agonist, Sog binds to and moves BMP ligand via facilitated diffusion to regions of Tolloid activity (Figure 1A). Tolloid then cleaves Sog, which releases BMP thus increasing peak BMP levels, a process altogether known as shuttling (Figure 1A) (Eldar et al., 2002; Marqués et al., 1997; Holley et al., 1996; Peluso et al., 2011; Shilo et al., 2013; Shimmi et al., 2005; Umulis et al., 2010). The shuttling mechanism is essential to Drosophila DV patterning, where Sog shuttles BMP ligand from lateral regions to dorsal regions (Figure 1A) (Eldar et al., 2002; Marqués et al., 1997; Holley et al., 1996; Peluso et al., 2011; Shilo et al., 2013; Shimmi et al., 2005; Umulis et al., 2010). This shuttling mechanism is required to steepen the BMP signaling gradient and specify the dorsal-most cell fates in the Drosophila embryo (Eldar et al., 2002; Marqués et al., 1997; Holley et al., 1996; Peluso et al., 2011; Shilo et al., 2013; Shimmi et al., 2005; Umulis et al., 2010). The shuttling of BMP ligand by Chordin has also been suggested to play a role in DV patterning in Echinoderms (Lapraz et al., 2009) and Nematostella (Genikhovich et al., 2015).
Figure 1
Potential Mechanisms of BMP Morphogen Gradient Formation.
(A) Cross-sectional view of the Drosophila embryo depicting Sog shuttling Dpp (the fly BMP ligand) dorsally. (B) Lateral view of the zebrafish embryo depicting Chordin (Chd) shuttling BMP ventrally. (C) Counter-Gradient: Chd diffuses ventrally to form a counter-gradient repressing BMP. (D) Shuttling: BMP bound to Chd is shuttled ventrally, where it is released by Tolloid cleavage. (E) Transcriptional: BMP stays where it is produced, mirroring the bmp expression gradient. (F) Source-sink: BMP diffuses from its source of ventral production to a sink of dorsal Chd.
It is unclear whether Chordin shuttles BMP in patterning vertebrate tissues. In Xenopus, the shuttling of a particular BMP ligand, ADMP, by Chordin was reported to play a role in DV axial patterning in the scaling of embryos (Ben-Zvi et al., 2008; Reversade and De Robertis, 2005). In the mouse, Chordin has been suggested to shuttle BMP ligand from where it is expressed in the intervertebral disc to its site of signaling in the vertebral body (Zakin et al., 2010). Mathematical models of zebrafish and Xenopus DV patterning have predicted that Chordin could shuttle BMP ligand (Ben-Zvi et al., 2008; Zhang et al., 2007). The transcriptional profiles of zebrafish BMP components at the onset of gastrulation resemble that of the Drosophila embryo (Dutko and Mullins, 2011; O'Connor et al., 2006). In Drosophila, sog is expressed ventral-laterally while the BMP ligand dpp is expressed dorsally (Figure 1A). Vertebrates have undergone a DV axis inversion with respect to arthropods (De Robertis and Sasai, 1996; Gerhart, 2000; Lacalli, 1995; Sander and Schmidt-Ott, 2004), thus chordin is expressed dorsally while bmp ligands are expressed ventrally (Figure 1B). However, whether Chordin acts as a BMP agonist by shuttling BMP ligand during DV patterning in zebrafish or other vertebrates has not been determined (Figure 1B).
In vertebrates, the mechanism by which the BMP ligands and antagonists shape this gradient is unclear. Several potential mechanisms have been proposed: 1) an inverse gradient of BMP antagonists imparts the shape of the BMP signaling gradient (Figure 1C) (Blitz et al., 2000; Connors et al., 1999; Little and Mullins, 2006; Thomsen, 1997), 2) BMP antagonists generate the peak BMP signaling levels by shuttling BMP ligand to these regions (Figure 1B,D) (Ben-Zvi et al., 2008; Shilo et al., 2013; Zhang et al., 2007), 3) the gradient shape mirrors the shape of the bmp expression domain (Figure 1E) (Ramel and Hill, 2013), and 4) the gradient is generated by BMP diffusing from its ventral source to a dorsal sink of BMP antagonists (Figure 1F). These mechanisms are not mutually exclusive and multiple may act in combination.
To identify the mechanism of BMP signaling gradient formation in the zebrafish embryo, we established a robust quantitative imaging method to directly measure the BMP signaling gradient. We integrated the results with a mathematical modeling approach, using the experiments to inform our model selection. The modeling then provided information on key parameters to measure to identify the mechanism by which the BMP signaling gradient is formed. We used phosphorylated Smad5 protein as a direct read-out for BMP signaling in both wild type (WT), chordin mutant, and chordin heterozygous embryos. We quantified nuclear phosphorylated-Smad5 (P-Smad5) fluorescent intensity across the entire embryo at single-cell resolution at different stages of development. Combining the P-Smad5 data with a computational screen of mathematical models showed that shuttling of BMP during DV patterning does not shape the gradient, and that a gradient of bmp transcription cannot account for the gradient of BMP signaling activity. From these results, we conclude that the signaling gradient patterning the vertebrate DV axis is generated by either a source-sink or counter-gradient mechanism. To discern between these mechanisms, we developed and measured the diffusion rate of a BMP2-Venus fusion protein in the zebrafish blastula and found that it is relatively mobile, which supports a large number of source-sink simulations, but far fewer counter-gradient simulations. Our results suggest that significant differences exist between the biophysical parameters of conserved proteins in zebrafish and Drosophila DV patterning. Through quantification and modeling, we present a new view of the mechanism that the BMP antagonists and ligand use to establish the BMP signaling gradient patterning the DV axis in zebrafish.
Results
Quantifying the Wild-Type signaling gradient
To measure the BMP signaling gradient, we quantified the levels of the BMP signal transducer P-Smad5 across the entire embryo at single cell resolution. Smad5 is directly phosphorylated by the BMP type I receptor in response to BMP signaling, and P-Smad5 concentration has been shown to linearly correlate with the concentration of BMP ligand in the Drosophila wing disc and S2 cells (Bollenbach et al., 2008; Serpe et al., 2008). Fixed embryos were whole-mount immunostained for P-Smad5 and imaged using a Line Scanning Confocal Microscope (Figure 2A–E). We developed a mounting and imaging protocol that minimized photo-bleaching, light scattering, and refractive index mismatch (see Materials and methods). We wrote a Matlab algorithm to identify all 8000 + nuclei centerpoints in each embryo in three dimensions, to remove populations unresponsive to P-Smad5 such as yolk syncytial nuclei and dividing cells (see Materials and methods), and to extract the P-Smad5 intensities associated with each nucleus (Figure 2A’–E’). Embryos were aligned by coherent point drift to a reference embryo to create ensembles of embryos suitable for statistical analysis (Myronenko et al., 2010a) (see Materials and methods). We used a band of cells around the margin of the embryo (Figure 2F’) to plot profiles from the dorsal-most to the ventral-most points to compare P-Smad5 gradient profiles between stages and between wild-type and mutant embryos (Figure 2F).
Figure 2
Dynamics of the WT P-Smad5 gradient across head and trunk patterning.
(A–E) Animal views of maximum projections (MP) of P-Smad5 stained individual embryos. (A’–E’) Animal views of nuclear intensities of all nuclei from the embryos shown above. (F) Average marginal intensities for 4.7–6.7 hpf (4.7: N = 3, 5.3: N = 4, 5.7: N = 13, 6.3: N = 11, 6.7: N = 4). Error bars indicate standard deviation. (G) Slope of the P-Smad5 gradients shown in panel F. Dotted line separates high slope (>0.5 a.u./deg) regions from low slope regions. (H) Change in cell number versus (vs.) developmental time of embryos fixed from a single cross and nuclei stained with Sytox Orange. (I) Cell number varies between different crosses of WT fish fixed at 5.7 hpf. (H,I) Gray dots are individual embryo cell counts. Red lines show the mean number of cells at a given time point, red boxes show 95% confidence interval, blue boxes show one standard deviation.
Our quantitative analysis revealed that the BMP gradient during DV patterning is quite dynamic. BMP signaling patterns prospective head and rostral trunk DV axial tissues during late blastula to mid-gastrula stages at ~5 to 7 hr post fertilization (hpf) in the zebrafish (Hashiguchi and Mullins, 2013; Kwon et al., 2010; Tuazon and Mullins, 2015; Tucker et al., 2008). We quantified the BMP signaling gradient at 30 min intervals across this period. We found that the ventral-most 30˚ undergoes about a 2-fold intensification from 4.7 to 6.7 hpf (Figure 2F). This is accompanied by a 3 to 5 fold increase in the slope of the gradient in ventrolateral regions of the embryo (0–75 degrees) over this 2 hr period (Figure 2G). Moreover, the lateral region encompassing the high slope (>0.5 A.U./degree) expands from a size of 20˚ to 75˚, meaning that by 6.7 hpf, nearly half the embryo falls within this high slope region. This contrasts with Drosophila DV patterning, where an initial broad, low-slope distribution of P-Mad is refined into a steep peak of BMP signaling covering only the dorsal-most 8% of the embryo (11 cell lengths) (Sutherland et al., 2003; Wang and Ferguson, 2005). This intensification of P-Mad is very rapid in Drosophila DV patterning, where P-Mad increases about three fold in the 30 min (min) between stages 5 and 6 (Ross et al., 2001; Sutherland et al., 2003; Wang and Ferguson, 2005), a process that we found is much slower in the zebrafish embryo: a 2-fold increase over a 2 hr period.
We then sought to determine if changes in cell number could account for the observed changes in the P-Smad5 gradient. We counted the number of cells in each embryo from a single timecourse and observed an approximately 70% increase in cell number from 4.7 to 6.7 hpf (Figure 2H). Cell nuclei do not change significantly in size during this time (Keller et al., 2008). The increase in cell number occurs throughout the embryo and is not restricted to a particular DV region, while the change in gradient amplitude (Figure 2F) and slope is restricted to the ventral third of the embryo (Figure 2G). Thus an increase in cell number does not account, via an unknown mechanism, for the increase in amplitude or slope. Additional support that cell number has little effect on gradient shape stems from the observation that the absolute number of cells at a given time point can vary by as much as 20% between different embryos within the same cross or between crosses (Figure 2I) with no detectable change in gradient shape or phenotype.
Mathematical modeling and computational screen of BMP gradient formation
We then performed a computational screen of mathematical models to investigate which gradient-forming mechanisms (Figure 1C–F) fit the WT P-Smad5 gradient profiles (Figure 2). To do so, we first needed to determine the expression domains of bmp, chordin, noggin, and tolloid to use for the mathematical model. We based the domain sizes on our own measurements (Figure 3), as well as published in situ hybridizations for bmp (Fürthauer et al., 2004; Ramel and Hill, 2013), chordin (Miller-Bertoglio et al., 1997), tolloid (Connors et al., 1999), and noggin (Dal-Pra et al., 2006). In animal-pole views of wholemount in situ hybridizations, we measured the chordin and noggin expression domain sizes to be 75 and 40 degrees in width, respectively (Figure 3A–C). We estimated the size of the bmp expression domain via wholemount in situ hybridizations (Figure 3D), however, bmp2b expression appeared graded and not as easily measured as chordin and noggin. We instead measured the bmp2b expression gradient at 5.7 hpf via fluorescent in situ hybridization on cross-sections of the DV margin (Figure 3F–H). We quantified the relative intensity of the fluorescent bmp2b in situ (Figure 3I black line) and used it to estimate the BMP production domain in our model (Figure 3I, blue line).
Figure 3
Measuring the bmp2b, chordin, and noggin expression domains.
Animal pole views of wholemount in situ hybridizations of the expression of (A) chd (N = 25), and (B) nog (N = 8) in WT embryos. (C) Measured domain size of chd and nog domains via wholemount in situ hybridization in WT and chd mutant embryos. (D–D’’’) bmp2b in chd ± embryos at 4.7 (N = 10), 5.3 (N = 15), 5.7 (N = 20), and 7 hpf (N = 16), and (E–E’’’) bmp2b expression in chd -/- embryos at 4.7 (N = 6), 5.3 (N = 16), 5.7 (N = 13), and 7 hpf (N = 12). (F–H) Fluorescent in situ hybridization (FISH) signal of bmp2b from a marginal slice at 5.7 hpf with a DAPI nuclear stain. Scale bars = 100 μm. (I) Quantification of FISH of bmp2b expression from ventral to dorsal (black line, N = 5) compared to the BMP production gradient used in the mathematical model (blue dotted line). Error bars indicate standard deviation.
We then developed a system of partial differential equations to model the interactions of BMP, Chordin, Noggin, and Tolloid (Supplementary Figure 4—figure supplement 4 and 1A). BMP, Chordin, Noggin, BMP-Chordin, and BMP-Noggin were modeled as diffusible species, while Tolloid was treated parametrically according to its domain of expression (Figure 4A). The zebrafish gastrula was reduced to a 1-dimensional (1D) half-circumference with a length of 700 µm. We selected a 1D description at the margin for simulation, since the inputs to the system based on gene expression are distributed symmetrically across the embryo, so a 1D model should largely reflect one in 3D, and the faster simulation time of a 1D model allowed us to explore far greater parameter space computationally, increasing our ability to discern biophysical constraints that distinguish between the mechanisms. Domains of production of BMP, Chordin, and Noggin were estimated as described (Figure 3, Figure 4D). The dissociation constants for BMP-Chordin and BMP-Noggin were set to 1 and 0.1 nM respectively, based on previously reported analysis (Piccolo et al., 1996; Troilo et al., 2014; Zimmerman et al., 1996). All remaining parameters (ie. the diffusion coefficients, production rates, decay rates, on and off binding rates) were varied over 4 orders of magnitude encompassing all biologically feasible values (Table 1).
Table 1
List of the parameter ranges used in the computational model-based screen.
Values range between the upper and lower bound. Note that the dissociation constant of BMP-Chd and BMP-Nog was held constant, but the on- and off- rates were allowed to vary.
| Parameter | Units | Symbol | Lower bound | Upper bound |
|---|---|---|---|---|
| BMP Production Rate | nM/s | ηB | 10−2 | 102 |
| BMP Decay Rate | 1/s | decB | 10−1 | 10−5 |
| BMP Diffusivity | µm2/s | DB | 10−2 | 102 |
| Chd Production Rate | nM/s | ηC | 10−2 | 102 |
| Chd Decay Rate | 1/s | decC | 10−1 | 10−5 |
| Chd Diffusivity | µm2/s | DC | 10−2 | 102 |
| Nog Production Rate | nM/s | ηN | 10−2 | 102 |
| Nog Decay Rate | 1/s | decN | 10−1 | 10−5 |
| Nog Diffusivity | µm2/s | DN | 10−2 | 102 |
| BMP-Nog Decay Rate | 1/s | decBN | 10−1 | 10−5 |
| BMP-Nog Diffusivity | µm2/s | DBN | 10−2 | 102 |
| BMP-Chd Decay Rate | 1/s | decBC | 10−1 | 10−5 |
| BMP-Chd Diffusivity | µm2/s | DBC | 10−2 | 102 |
| Chd Degradation by Tld | 1/s | λC | 100 | 10−4 |
| BMP-Chd Degradation by Tld | 1/s | λBC | 100 | 10−4 |
| Length of the Embryo | µm | - | 700 | 700 |
| Length of the Chd domain (from dorsal) | µm | - | 145 | 145 |
| Length of the Nog domain (from dorsal) | µm | - | 78 | 78 |
| Length of the Tolloid domain (from ventral) | µm | - | 400 | 400 |
| Dissociation Constant of BMP-Chd | nM | - | 1 | 1 |
| Dissociation Constant of BMP-Nog | nM | - | 0.1 | 0.1 |
| Time | min | t | 130 | 130 |
Figure 4 with 1 supplement see all
Creation of a mathematical model of BMP gradient formation.
(A) Depiction of the species and binding possibilities modeled. (B) BMP distributions of 10 individual model solutions (black dotted lines) plotted against the WT 5.7 hpf P-Smad5 gradient (blue line). Error bars indicate standard deviation. (C) Flow-chart of model mechanism classification. (C’) BMP mass balance from model labeled to indicate which terms contribute to the source-sink, counter-gradient, and transcriptional mechanisms at each point. (C”) Shuttling mechanism was defined by a 20% decrease at the ventral-most point in chd LOF compared to WT. (D) Expression domains of bmp (blue), tld (purple), chd (red), and nog (yellow) used in the model. (E) Expression domains of bmp (blue), tld (purple), chd (red), and nog (yellow) used for the alternative scenario where the bmp expression domain mirrors the measured P-Smad5 gradient. (D’, E’) Pie chart showing how many parameter combinations fit the WT data (blue) and how many failed to do so (grey). (D”, E”) Pie chart showing how many parameter combinations were classified to have a source-sink (blue), counter-gradient (red), transcriptional (orange), or shuttling (green) mechanism.
The equations were solved for the developmental window from ~3.5 hpf to ~5.7 hpf, since bmp and chordin are first expressed shortly after the mid-blastula-transition at 3 hpf (Koos and Ho, 1999; Leung et al., 2003; Shimizu et al., 2000; Solnica-Krezel et al., 2001). The equations were simulated 1,000,000 times, each time with a different combination of randomly selected parameters. Each parameter combination was then re-simulated without Chordin or Noggin to predict the BMP signaling gradient in a chordin or noggin loss-of-function (LOF) scenario. We then selected simulations that generated BMP profiles that fit our measured P-Smad5 gradient at 5.7 hpf with a normalized root mean squared deviation of 8% or less (Figure 4B). We also eliminated simulations that had significantly different BMP profiles when Noggin production was set to 0, since the loss of Noggin does not affect DV patterning in zebrafish or Xenopus (Dal-Pra et al., 2006; Khokha et al., 2005).
All simulation results that fit our data were classified into categories based on the biophysical process that dominated formation of the gradient shape: shuttling, source-sink, counter-gradient, or transcriptional (Figure 4C). We discerned between source-sink, counter-gradient, and transcriptional mechanisms by examining the balance of binding, diffusion, decay, and accumulation processes in the partial differential equation for the BMP species (Figure 4C–C’’). If 80% of the BMP ligand was degraded where it was produced or accumulated there, the simulation was classified as transcriptional (Figure 4C,C’). If the majority of BMP diffused away from its site of production rather than being bound by Chordin, the simulation was considered a source-sink mechanism (Figure 4C,C’). Conversely, if the majority of BMP was bound at its site of production by Chordin, the simulation was classified as a Chordin counter-gradient mechanism (Figure 4C,C’). We classified a simulation as shuttling if the ventral-most point in the predicted chordin-/- BMP profile was at least 20% lower than in WT, as shuttling by the antagonist leads to a net accumulation of ligand in the ventral-most region (Figure 4C,C’’). By comparison, shuttling in Drosophila accounts more significantly to the peak level, with a 50% decrease in the peak P-Mad level when the chordin homolog sog is deficient (Mizutani et al., 2005; Peluso et al., 2011; Sutherland et al., 2003). Shuttling of 20% is about one standard deviation greater than our measured embryo-to-embryo variability for peak P-Smad5 signaling levels (Figure 4B).
Multiple classes of mechanisms generated simulations fitting our WT data, including an antagonist counter-gradient, source-sink, and shuttling. Of the 1,000,000 randomly picked parameter combinations, 15,142 fit the experimentally measured WT signaling gradient (Figure 4D’). Among those that fit, 13,382 were classified as source-sink, 1710 as counter-gradient, and 50 as shuttling (Figure 4D”). Notably, no transcriptional simulations were found, because our measured bmp expression profile (Figure 3I) did not match our measured WT BMP signaling gradient (Figure 2), and therefore BMP must diffuse from its site of production or be bound by Chordin to fit our measured signaling gradient.
The mean of our bmp2b expression profile did not reflect the P-Smad5 gradient at 5.7 hpf (Figures 2F and 3I), however, the bmp2b expression profile displayed variability from embryo to embryo. To account for the possibility that the bmp expression profile reflects the WT P-Smad5 gradient, which fell within one standard deviation of the mean (Figures 2F and 3I), we simulated the mathematical model 250,000 additional times with a bmp expression input matching the WT P-Smad5 gradient (Figure 4E). We found that the transcriptional mechanism was the most abundant mechanism among the simulations, comprising 83,747 of the simulations (Figure 4E’,E’’). Surprisingly, no counter-gradient solutions were found when bmp expression matched P-Smad5. This is because Chordin binding to BMP would interfere with the shape of the BMP protein gradient, causing it to no longer match the measured BMP signaling gradient. Thus the transcriptional and counter-gradient mechanisms are incompatible with each other.
Constraints on biophysical parameters imposed by the WT gradient
Each mechanism required specific biophysical parameters to fit the experimentally measured DV signaling gradient. The source-sink mechanism required BMP to have a high diffusion rate, so BMP could diffuse to a dorsally-localized sink of antagonists (Figure 5A). The counter-gradient mechanism required Chordin to have a high diffusion rate, so Chordin could diffuse ventrally to generate an antagonist gradient (Figure 5B). When either BMP or Chordin had moderately high diffusion rates (e.g. near 1 um2/s), decay rates were low to allow each more time to diffuse dorsally or ventrally, respectively. When BMP and Chordin range are plotted on the same axis, the segregation of the source-sink and counter-gradient mechanisms based on range is readily apparent (Figure 5C). The shuttling mechanism required that BMP-Chordin have a high diffusion rate and low decay rate in order to diffuse ventrally where Tolloid cleaves Chordin, which then releases BMP (Figure 5D). Noggin diffusion and decay was not constrained in any of the three mechanisms (Figure 5E).
Figure 5
Biophysical values of individual simulations that fit the WT P-Smad5 gradient.
(A–L) Scatter plots comparing biophysical parameters of 1000 randomly selected solutions classified by mechanism that fit the WT data. Combinations that failed to fit the WT P-Smad5 gradient are small grey dots. We plot solutions as large circles colored according to their mechanism, which is based on definitions outlined in Figure 4C: counter-gradient (red), source-sink (blue), transcriptional (orange), or shuttling (green). We plotted additional shuttling solutions in order to better illustrate trends. (A–F) Simulations using domains displayed in Figure 4D. (G–L) Simulations using domains displayed in Figure 4E. (A,G) BMP diffusivity vs. BMP decay rate. (B,H) Chd diffusivity vs. Chd decay rate, which includesthe rate of Chd cleavage by Tld. (C,I) Range was estimated as sqrt(diffusivity/decay). (D,J) Diffusivity of BMP bound to Chd vs. decay rate of BMP bound to Chd. (E,K) Range of Nog protein. (F,L) Chd and BMP-Chd cleavage rate by Tld.
For the simulations performed where the bmp expression domain matched the signaling gradient, distinct biophysical requirements were also observed for each mechanism. The source-sink mechanism required BMP to have a high range, while the transcriptional mechanism required BMP to have a lower range (Figure 5I). The source-sink mechanism required either high diffusivity with predominantly high decay rates or low diffusivity with low decay rates, while the transcriptional mechanism filled a near complementary parameter space with progressively higher BMP diffusivity necessitating higher decay rates (Figure 5G), which allows the P-Smad5 gradient to match the bmp expression profile. Neither of these mechanisms constrained Chordin range (Figure 5H,I). Although only four shuttling mechanism solutions were identified, they required both BMP-Chordin range to be high and its decay rate to be low (Figure 5J). The shuttling solutions also required BMP range to be high (Figure 5G,I), because if Chordin moved ventrally to bind BMP in this scenario, it would interfere with the bmp expression gradient matching the BMP signaling gradient. In these simulations BMP moves dorsally to form the BMP-Chordin species that ultimately is shuttled back ventrally. The remaining parameters required for shuttling were similar for the two simulations (Figure 5E,F,K,L).
The chordin mutant gradient shows no evidence of shuttling
To determine whether Chordin shuttling of BMP ligand plays a functionally relevant role in generating the ventral P-Smad5 peak in zebrafish, as it does in Drosophila, we quantified the P-Smad5 gradient of chordin mutant embryos over a developmental time series (Figure 6A,B). If it does play a role, then we expect the ventral P-Smad5 peak to be reduced in chordin mutants compared to WT embryos. We found that the P-Smad5 gradient in chordin mutants showed a statistically significant increase in lateral regions of the embryo at the four time-points examined from 4.7 to 6.3 hpf (Figure 6C–F, Figure 6—figure supplement 1A). Importantly, no decrease in P-Smad5 was observed in the ventral region of chordin mutant embryos or in any region of the gradient. These results indicate that, unlike the Drosophila homolog Sog, Chordin plays no significant BMP shuttling role during zebrafish DV patterning. It is worth noting that in many simulations, small amounts of BMP ligand are shuttled short distances but do not impact the gradient significantly, and thus are not classified as shuttling. The P-Smad5 gradient in chordin mutants shows that this is minimal in zebrafish.
Figure 6 with 1 supplement see all
Effect of Chd on gradient shape and ligand shuttling.
(A,B) Animal views of average intensities from each time-point in (A) WT (4.7: N = 3, 5.3: N = 4, 5.7: N = 13, 6.3: N = 11) and (B) chd mutant (4.7: N = 3, 5.3: N = 5, 5.7: N = 11, 6.3: N = 9) embryos. (C–F) Average marginal intensities for WT (blue) and chd mutant (red) embryos from 5.7 to 6.3 hpf. (G) Average marginal intensities for WT (blue, N = 4) and bmp2/7 RNA injected embryos at 5.7 (grey, N = 4) and 6.3 hpf (black, N = 5). Error bars indicate standard deviation. (H,I) Fully ventralized (V5) embryos injected with 6 pg of bmp7a RNA and 12 pg of bmp2b RNA vs uninjected WT siblings. (J) WT (N = 9) vs chd+/- (N = 10) at 5.7 hpf.
Interestingly, the loss of chordin did not cause an increase in the ventral-most P-Smad5 peak level either (Figure 6C–F,S1A), suggesting that Chordin does not actively block BMP signaling there. However, Smad5 or another signal transducing component could be limiting in the ventral-most cells of WT embryos, rendering them unresponsive to further increases in free ligand. To investigate this possibility, we overexpressed Bmp2/7 ligand in WT embryos at a level that fully ventralizes (V5) them at 24 hpf (Figure 6H,I). We quantified P-Smad5 at 5.7 and 6.5 hpf and found that the gradient showed a significant increase in signaling embryo-wide over WT siblings, including in the ventral-most region (Figure 6G). These results indicate that BMP signaling in ventral regions is not saturated in WT embryos and that Chordin does not regulate the peak P-Smad5 levels by promoting or inhibiting signaling at these stages.
We next tested whether the P-Smad5 gradient is robust to heterozygosity of chordin. The Drosophila P-Mad gradient shows some small changes in sog heterozygotes (sog is the chordin homolog) compared to WT (Eldar et al., 2002; Umulis et al., 2010). In zebrafish chordin heterozygotes do not show any DV patterning phenotype at 24 hpf (Hammerschmidt et al., 1996). Consistent with this, we found that the chordin heterozygous and WT BMP signaling gradients were indistinguishable at 5.7 hpf (Figure 6J). Therefore, the BMP signaling gradient in zebrafish is robust to chordin heterozygosity and to a potential 50% decrease in Chordin levels, but not to the complete loss of Chordin.
Constraints on computational models imposed by the chordin mutant gradient
We then constrained the mathematical model-based computational screen with the WT, the chordin mutant and heterozygote P-Smad5 gradients at 5.7 hpf, assuming that heterozygotes have half the Chordin level of WT. We determined which mechanisms and simulations remained compatible with these new results (Figure 7A,B). We eliminated simulations that deviated by more than 8% from the chordin heterozygous or homozygous P-Smad5 gradients (Figure 7C–D). Many mathematical model simulations that fit the WT BMP signaling gradient no longer fit our chordin heterozygous and homozygous mutant data. Of the 15,142 simulations that fit the WT gradient alone, only 4059 fit the WT, chordin +/-, and chordin mutant gradients (Figure 7E–E’). Of those, all were either source-sink or counter-gradient mechanisms (Figure 7E’). All simulations classified as shuttling had chordin LOF BMP distributions that deviated from the measured chordin mutant P-Smad5 gradient by more than 17% (Figure 7C). Many, but not all, simulations classified as counter-gradient deviated by more than 8% in their BMP distributions from the measured chordin heterozygous P-Smad5 gradient (Figure 7D). Fitting the chd LOF and chd heterozygous data eliminated more counter-gradient than source-sink solutions, with 28% of the source-sink solutions remaining but only 15% of the counter-gradient solutions remaining.
Figure 7
Biophysical values of individual simulations that fit both the WT and chd LOF P-Smad5 gradients.
(A) BMP distributions of 5 individual modeling solutions (WT: black dotted lines, chd LOF: grey dotted lines) plotted against WT (blue line) and chd LOF (red line) 5.7 hpf P-Smad5 gradients. Error bars indicate standard deviation of experimental P-Smad5 intensity. (B) BMP distributions of 5 individual modeling solutions (WT: black dotted lines, chd +/-: grey dotted lines) plotted against WT (blue line) and chd +/- (green line) 5.7 hpf P-Smad5 gradients. Error bars indicate standard deviation of experimental P-Smad5 intensity. (C,D,F–J,L, M) 1000 randomly selected parameter combinations capable of fitting both the WT data, chd +/-, chd LOF data classified by mechanism. Larger circular points fit the WT P-Smad gradient and are colored based on their mechanism according to the definitions outlined in Figure 4C: counter-gradient (red), source-sink (blue), transcriptional (orange), or shuttling (green). Combinations that failed to fit the WT P-Smad5 gradient are small grey dots. (C–D) Normalized Root Mean Squared Deviation (NRMSD) between the measured P-Smad5 and the model BMP distributions. Black dotted lines mark the 8% threshold. (C) Comparing WT and chd LOF. (D) Comparing WT and chd +/-. (E) Parameter combinations that fit both the WT and chd LOF data (blue) and those that failed to do so (grey). (E’) parameter combinations were classified to have a source-sink (blue), counter-gradient (red), or shuttling (green) mechanism. (F–J) Simulation using the bmp expression domain displayed in Figure 4D. (L,M) Simulation using the bmp expression domain displayed in Figure 4E. (F,L) BMP diffusivity vs. BMP decay rate. Green dotted line marks the BMP diffusivity we measured using FRAP (4.4 µm2/s). (G,M) Range was estimated as sqrt(diffusivity/decay). (H) Chd diffusivity vs. Chd decay rate plus the rate of Chd cleavage by Tld. (I) Rate of Chd cleavage by Tld vs rate of BMP-Chd cleavage by Tld. (J) BMP-Chd diffusivity for source-sink and counter-gradient simulation solutions. (K) Pie chart showing the parameter combinations that fit the WT data (blue) or failed to do so (grey) for the alternative scenario where the bmp expression domain mirrors the measured P-Smad5 gradient (Figure 3I). (K’) Pie chart of the solutions that had a source-sink (blue), counter-gradient (red), transcriptional (orange), or shuttling (green) mechanism.
The remaining mechanisms required different and specific combinations of biophysical parameters. The source-sink solutions required a high BMP range of 60+ µm with a diffusivity above 1 µm2/s (Figure 7F,G). The counter-gradient mechanism required a lower BMP range, less than 60 µm, a high Chordin range above 40 µm with a diffusivity above 2 µm2/s (Figure 7G,H). Consistent with this, very high rates of Chordin cleavage by Tolloid restricted Chordin range and therefore was not compatible with the counter-gradient mechanism (Figure 7I). While BMP diffusivity needed to be high in source-sink solutions, BMP-Chd diffusivity was not restricted in either the source-sink or counter-gradient solutions (Figure 7J).
We then tested whether the transcriptional mechanism could also fit the chordin mutant data. We used the simulation in which the bmp expression profile matched the WT P-Smad5 gradient (Figure 4E). Of the 148,646 simulations that fit the WT data alone, 227 fit the WT, chordin heterozygous, and chordin mutant gradients, all of which were source-sink (Figure 7K,K’). These source-sink simulations required a high BMP diffusivity and range (Figure 7L,M), similar to what was observed in the simulation with the measured bmp expression profile (Figure 7F,G). However, while 83,747 simulations fit a transcriptional mechanism with the WT data alone, none fit the WT, chordin heterozygous, and chordin mutant gradients (Figure 7K,K’). This is because only a change in the bmp expression domain in the chordin mutant allows a transcriptional mechanism to fit both the WT and chordin mutant data.
The bmp expression domain is known to become responsive to BMP signaling, creating a positive feedback loop during gastrulation (Hammerschmidt et al., 1996; Nguyen et al., 1998; Schmid et al., 2000). Gastrulation begins in zebrafish at 6 hpf. While the initial bmp expression domains are established independently of BMP feedback, a BMP feedback loop becomes active with reported onset times ranging from ~5.5 to 6.5 hpf (Kishimoto et al., 1997; Miller-Bertoglio et al., 1999; Ramel and Hill, 2013; Schmid et al., 2000). To test whether the bmp expression domain changes in chordin mutants at 4.7 to 5.7 hpf, we compared the bmp2b domain size in sibling chordin-/- and chordin +/- embryos (Figure 3D,E). chordin heterozygotes display a WT phenotype (Hammerschmidt et al., 1996; Miller-Bertoglio et al., 1999) and we found that they also display a WT P-Smad5 gradient (Figure 6J). There was no discernable difference in the bmp2b domain size at 4.7, 5.3, or 5.7 hpf (Figure 3D,E), indicating that bmp transcriptional feedback is not active before 5.7 hpf. Similarly, the noggin expression domain size did not change in chordin mutants before 5.7 hpf (Figure 3C). Therefore, the increase in BMP signaling activity observed already at 4.7 hpf in chordin mutants (Figure 6A–C) precedes a change in bmp expression, showing that the transcriptional mechanism cannot account for the P-Smad5 gradient profiles prior to 5.7 hpf. We observed a change in bmp2b expression by 7 hpf, consistent with previous findings that bmp transcriptional feedback activates after gastrulation begins (Figure 3D’’’,E’’’) (Kishimoto et al., 1997; Miller-Bertoglio et al., 1999; Ramel and Hill, 2013; Schmid et al., 2000).
Fluorescence recovery after photobleaching to measure BMP diffusivity
The combination of WT and chordin mutant P-Smad5 gradients limits the number of computationally-derived model simulations and, importantly, reduces the number of mechanisms to two: source-sink and counter-gradient. We and others have largely purported the counter-gradient mechanism as acting in vertebrate DV patterning (Blitz et al., 2000; Hama and Weinstein, 2001; Little and Mullins, 2006; Sasai and De Robertis, 1997; Thomsen, 1997). To our surprise, the source-sink modeling simulations emerged more frequently within our computational screen than the counter-gradient simulations (Figure 7E’). To test the source-sink mechanism further, we investigated BMP ligand diffusivity, a biophysical parameter that must be high in this mechanism (Figure 7F). To test if BMP diffusivity excludes or supports the source-sink mechanism, we measured the effective diffusivity of the Bmp2b ligand using fluorescence recovery after photobleaching (FRAP).
We tagged Bmp2b by inserting the coding sequence of the fluorescent protein Venus between the pro- and mature coding domains of bmp2b. When mRNA of bmp2b-venus was injected into 1-cell stage zebrafish embryos, we detected both the pro- and mature domains of Bmp2b-Venus protein on a western blot (Figure 8A, black arrows, n = 3). We did not detect protein with a molecular weight equal to Venus alone in the embryos injected with bmp2b-venus, indicating that the Venus tag was not cleaved from Bmp2b during post-translational processing (Figure 8A, red arrow). The Venus tag did not interfere with Bmp2b activity, since the injected mRNA significantly ventralized WT embryos (Figure 8B, Row 1). To further assess the activity and range of the Bmp2b-Venus chimera, we tested if Bmp2b-Venus could rescue embryos deficient in Bmp2b. As previously reported for bmp2b RNA (Nguyen et al., 1998), we could rescue Bmp2b deficient embryos to a WT phenotype by injecting bmp2b-venus RNA (Figure 8B, Rows 5–7).
Figure 8
Measuring Bmp2b-Venus diffusivity via FRAP.
(A) Detection of Bmp2b-Venus and secreted Venus proteins by western blot. Embryos were injected with bmp2b-venus mRNA (250 pg) or secreted-Venus mRNA (200 pg) at the one-cell stage. Protein lysates were prepared at late blastula stage. In the Bmp2b-Venus overexpression sample, two major protein bands were detected by Venus antibody (black arrows). The larger molecular weight protein is the pro- and mature domains of Bmp2b with Venus protein (669 amino acids (AA),~74 KDa). The smaller protein is the mature domain of Bmp2b with Venus protein (376 AA,~41 KDa). The secreted Venus protein (248 AA,~27 KDa) is also detected in the secreted-Venus overexpression sample (red arrow). β-actin was used as a loading control. (B) 24 hpf phenotypes of embryos injected with the bmp2b-venus construct used for FRAP experiments, controls, and rescue. Dorsalization was classified as C5: Loss of all ventral structures; C4-C3: Loss of, or truncated tail; C2-C1: Loss of ventral tail fin. Ventralization is classified as V1: reduction is eye size; V2-V3: the eyes, notochord, and anterior brain are partially or completely absent; or V4-V5: complete loss of all dorsal structures. Fluorescent BMP-Venus (C) or Venus (D) recovery after photobleaching for 20 min. (E–G) Plots of fluorescent intensity recovery in the extracellular region. Bold lines are mean curves, thin lines are raw intensity data. (H) BMP diffusivity vs. BMP decay rate for simulations that fit WT, chd +/-, and chd -/- P-Smad5 profiles and were within 2 µm2/s of 4.4 µm2/s. Large blue circles are simulations classified as source-sink, red are counter-gradient, and small grey dots failed to fit the measured P-Smad5 profiles. (I) The mean BMP and Chd concentrations in all solutions that fit the WT, chd-/-, and chd +/- P-Smad5 data and within a diffusivity of 2.4 and 6.4 µm2/s that are also robust to uniform Chd production.
To perform the FRAP, we injected bmp2b-venus mRNA into a single blastomere at the 8 cell stage (Figure 8B, Row 3) and then photobleached a 160 µm cube of cells in 4.3 hpf embryos (Figure 8C). We then measured recovery of fluorescence over one hour. To ensure we only recorded extracellular Bmp2b-Venus, we photobleached a region away from the cells producing Bmp2b-Venus. The bleached region recovered fluorescence to its initial level in ~30 min (Figure 8C,E), corresponding to a measured Bmp2b-Venus effective diffusivity of 4.4 ± 0.4 µm2/s (SEM (standard error of the mean), n = 5). To ensure that we were only measuring the diffusivity of Bmp2b-Venus alone and not Bmp2b-Venus bound to Chordin, we repeated the FRAP experiment in Chordin deficient embryos (Figure 8B, Rows 3–4). Again, the bleached region recovered fluorescence to its initial level in ~30 min (Figure 8F), corresponding to a measured Bmp2b-Venus effective diffusivity of 4.0 ± 0.5 µm2/s (SEM, n = 5). To determine the extent to which Venus limited Bmp2b diffusion, we measured the diffusivity of Venus alone. The bleached region recovered fluorescence much more rapidly, and neared its initial level in 5 min (Figure 8D,G), corresponding to a measured Venus effective diffusivity of 16.3 ± 2.2 µm2/s (SEM, n = 5). This faster rate of Venus diffusion indicates that it per se does not reduce BMP2 diffusion.
A measured BMP diffusivity of ~4 µm2/s fits a large portion of the source-sink modeling simulations. In fact, 1421 source-sink simulations have diffusivities within 2 µm2/s of our measured diffusivity (Figure 8H). In contrast, only 31 counter-gradient simulations were within 2 µm2/s of our measured diffusivity, and these simulations have very high BMP decay rates (above 10−3/s) (Figure 8H). A decay rate of that magnitude would cause the half-life of BMP ligand in the embryo to be very short, <10 min for decay rates above 1 × 10−3/s. We can infer the BMP lifetime in the embryo based on our previous studies. We found previously that a pulse of injected BMP ligand protein in a bmp deficient embryo sustains P-Smad5 for at least 1.5 hr after injection (Little and Mullins, 2009). In another study, we found that P-Smad5 can persist for a maximum of 40–60 min in the absence of BMP signaling (Tucker et al., 2008). Thus, P-Smad5 persistence 1.5 hr after a BMP protein pulse indicates that the ligand remains for at least 30 to 50 min after injection, if P-Smad5 can persist for 40–60 min in the absence of signaling. A BMP ligand lifetime of 30 to 50 min is inconsistent with the low BMP half-lives required for the counter-gradient simulations, but consistent with the BMP half-lives of source-sink simulations (Figure 8H), suggesting that the source-sink mechanism much more likely establishes the BMP signaling gradient patterning the zebrafish DV axis.
Finally, we asked whether the source-sink and counter-gradient mechanisms are robust to uniform chordin expression. Uniformly expressing chordin by mRNA injection into one-cell stage chordin mutant embryos has been used to rescue chordin mutant embryos to adulthood (Fisher and Halpern, 1999), including in our study here. We simulated the system with ubiquitous Chordin production and determined that 426 of the 1452 solutions that fit our WT, chordin-/-, and chordin +/- data with a BMP diffusivity within 2 um2/s of our measured 4.4 um2/s retain a WT BMP gradient when Chordin is uniformly produced. We did not titrate the Chordin production rate, but had we done so, more simulations may have fit. Interestingly, the 426 remaining solutions are all source-sink mechanisms with gradients of Chordin that are high dorsally and low ventrally (Figure 8I).
Comparing zebrafish and Drosophila DV patterning
Our results show that the BMP signaling gradient patterning the zebrafish DV axis is markedly different from the one patterning the Drosophila DV axis. The zebrafish BMP signaling gradient is broad, reaching half of its maximum at ~40% of the total DV axis length (Figure 2F). In contrast, the Drosophila gradient is incredibly steep, reaching half of its peak at only ~10% of the total embryo DV axis length (Figure 9A) (Peluso et al., 2011; Sutherland et al., 2003). Similarly, the loss of the main BMP antagonist in either organism, Chordin or Sog, causes markedly different effects on the BMP signaling gradient (Figure 9A) (Mizutani et al., 2005; Peluso et al., 2011; Sutherland et al., 2003).
Figure 9
Comparing Zebrafish and Drosophila-like solutions.
(A) Depiction of the BMP gradients patterning the Drosophila and zebrafish DV axis. The Drosophila DV axis has been flipped to match the zebrafish. Solid lines are WT. Dotted lines are chd or sog LOF. (B) List of homologous proteins involved in DV patterning of zebrafish and Drosophila. (C–F) Solutions able to fit WT and chd LOF zebrafish data (blue) vs. solutions capable of fitting Drosophila-like WT and Drosophila-like sog LOF gradients (red). Parameter combinations that failed to fit either are represented as small grey dots. (C) Chd diffusivity vs. Chd decay rate, which includes the rate of Chd cleavage by Tld. (D) Diffusivity of BMP bound to Chd vs. decay rate of BMP bound to Chd. (E) Range was estimated as sqrt(diffusivity/decay). (F) Cleavage rate of Chd and BMP-Chd by Tolloid. (G) BMP diffusivity vs. BMP decay rate.
Zebrafish and Drosophila DV patterning differ in both length-scale and time-scale. The Drosophila embryo has a 250 µm half-circumference, while the zebrafish embryo has a 700 µm half-circumference. The zebrafish gradient is established gradually in ~2–3 hr and maintained for several hours (Ramel and Hill, 2013; Tucker et al., 2008), whereas the Drosophila BMP signaling gradient is established and patterns DV tissues in ~1 hr (Dorfman et al., 2001; Wang and Ferguson, 2005). Given these differences, we sought to determine if Drosophila-like shuttling simulations could exist with zebrafish time- and length-scales, and if so, how the biophysical parameters of the components would differ from those consistent with the WT and chordin mutant P-Smad5 gradients (Figure 7).
In the 1,000,000 random simulations tested, we found many simulations that could generate a steep gradient with extensive shuttling in zebrafish. Simulations were considered to be Drosophila-like if their WT gradient reached its half-maximum <10% of the total embryo circumference and the ventral peak of the chordin mutant was 50% lower than the ventral peak level of the WT curve (Mizutani et al., 2005). We found 251 simulations that fit the Drosophila-like WT and chordin mutant signaling gradients. We also excluded simulations with excessive BMP-Noggin interaction, as Drosophila does not possess noggin homologs (Figure 9B).
The Drosophila-like simulations required a very mobile Chordin and BMP-Chordin species. Drosophila-like simulations required Chordin to have a high range to move to encounter the BMP in the ventral region (Figure 9C,E). Similarly, BMP-Chordin needed to have a high range so it could be shuttled a sufficient distance towards the ventral-most region of the embryo (Figure 9D,E). The cleavage of free Chordin by Tolloid needed to be low to allow Chordin range to remain high (Figure 9F). Conversely, the cleavage of bound Chordin needed to be high to release the BMP from Chordin (Figure 9F). Chordin range must be high to allow the formation of a counter-gradient to block signaling in the lateral regions of the embryo (Figure 9E). Conversely, BMP range was relatively unrestricted in Drosophila-like simulations, as the shuttling mechanism relies more on BMP-Chordin mobility than BMP mobility (Figure 9G).
Discussion
Here we have quantified the BMP signaling gradient in WT and chordin zebrafish mutants by measuring with high precision the P-Smad5 immunofluorescence level in all ~8,000 + nuclei of an embryo, with high reproducibility within and between embryos at multiple developmental stages. We then used these data to inform a computational model-based screen of over 1,250,000 combinations of biophysical parameters of the major extracellular BMP modulators. We defined mathematical criteria to distinguish between four widely proposed mechanisms to set up the BMP signaling gradient. Our computational model-based screen excludes the shuttling and transcriptional mechanisms as possibilities for establishing our measured WT and chordin mutant P-Smad5 profiles, providing compelling evidence that the BMP signaling gradient patterning the zebrafish DV axis is established by either a counter-gradient or source-sink mechanism. We further determined that the effective diffusivity of the BMP ligand in the zebrafish embryo is relatively fast, consistent with 1421 source-sink simulations but only 31 counter-gradient ones (Figure 8H). Comparison of models that satisfy zebrafish or Drosophila-like gradient profiles suggests that the range of BMP-Chordin differs between zebrafish and Drosophila-like DV patterning mechanisms (Figure 9E), and the Drosophila-like patterning mechanism requires restricted Tolloid degradation rates for BMP-Chordin and Chordin, which were not observed for the zebrafish patterning mechanism (Figure 9F).
Fish vs. flies: A mechanism diverged
The shape of the BMP signaling gradient differs greatly between Drosophila and zebrafish DV patterning (Figure 9A). The parameters that drive the most significant difference between the Drosophila-like and zebrafish simulations are the mobility and processing rates of BMP-Chordin (Figures 7I and 9D,F). For shuttling to be possible in a system with a broad peak of signaling, as it is in zebrafish, the BMP-Chordin cleavage rate needed to be low enough to allow BMP-Chordin to move farther and distribute BMP over a larger region (Figure 5F,L). For shuttling to be possible in a system with a tight peak of signaling as it is in Drosophila, BMP-Chordin cleavage needed to be rapid to release BMP-Chordin over a smaller region (Figure 9F). We show that a shuttling mechanism is not functioning in zebrafish, indicating that this delicate balance of BMP-Chordin mobility and Tolloid cleavage has been lost or did not emerge in vertebrate DV patterning.
Sog and its vertebrate homolog Chordin differ in how they are processed by the metalloprotease Tolloid depending on whether it is bound to BMP ligand. Chordin can be cleaved by Tolloid whether bound to BMP ligand or not (Piccolo et al., 1997), while Sog is only cleaved when bound to BMP (Marqués et al., 1997). Interestingly, when Sog is mutated to allow it to be processed by Tolloid regardless of BMP binding, the shuttling of BMP-Sog complexes in flies is greatly reduced (Peluso et al., 2011), suggesting that this attribute of Sog is necessary for effective shuttling. However, surprisingly, we found numerous shuttling simulations in our zebrafish modeling-based screen with the opposite properties, high Chordin cleavage rates and low BMP-Chordin cleavage rates (Figure 5F,L). The requirement for preferential cleavage of BMP-Chordin by Tolloid only emerged when we screened for Drosophila-like simulations, in which shuttling was generating a tight peak of BMP signal (Figure 9A,F). This suggests that the preferential processing of BMP-Sog by Tolloid seen in Drosophila is not an inherent requirement of the shuttling mechanism, but may instead be a result of both the requirement to facilitate shuttling and to generate a steep gradient.
In Drosophila DV patterning, the BMP signaling gradient is so steep that its base falls well within the region of bmp expression, far from the sog/chordin expression domain (Figure 9A) (Francois et al., 1994; Holley et al., 1995). To suppress lateral BMP signaling and form the Drosophila-like simulations seen in our mathematical model-based screen (Figure 9A), Sog/Chordin needed to have a high range to diffuse far from its site of expression to inhibit BMP signaling over most of the bmp expression domain (Figure 9C,E). Therefore, the degradation of free Sog/Chordin by Tolloid needed to be low (Figure 9F). However, to generate a narrow peak of BMP signaling (Figure 9A), BMP-Chordin cleavage by Tolloid needed to be high (Figure 9F). Therefore, the requirement to preserve the range of action of free Chordin, combined with the requirement to rapidly cleave BMP-Chordin to generate a steep peak of signaling may explain why the preferential cleavage of BMP-Sog by Tolloid is needed for shuttling in Drosophila.
Comparing the source-sink and counter-gradient mechanisms
While the shuttling mechanism relies on the movement of the bound BMP-Chordin complex, the source-sink and counter-gradient mechanisms rely on the movement of unbound BMP and Chordin. The source-sink mechanism relies on BMP diffusing dorsally to bind Chordin. Conversely, the counter-gradient mechanism relies on Chordin diffusing ventrally to bind BMP. Consistent with this, we found that the majority of simulations consistent with the source-sink mechanism require a high BMP range and a high BMP diffusivity (above 1 µm2/s), while the counter-gradient mechanism requires a high Chordin range and high Chordin diffusivity (above 1 µm2/s) (Figures 5 and 7).
To illustrate the distinct manners by which a source-sink and counter-gradient mechanism would generate the zebrafish BMP signaling gradients observed, we graphically display in Figure 10 the relative contributions of BMP and Chordin diffusion and BMP-Chordin binding to gradient formation for the median simulations that fit our WT, chordin mutant, chordin +/-, and Bmp2b-Venus FRAP data. The primary differences between the source-sink and counter-gradient mechanisms manifest in the relative amount of Chordin protein that diffuses ventrally into the bmp expression domain and the primary role of Chordin in forming the BMP gradient. A counter-gradient mechanism leads to higher levels of Chordin that extend over a greater region of the ventral bmp expression domain compared to a source-sink mechanism (Figure 10A). Counter-gradient and source-sink mechanisms also differ significantly in where Chordin binds and inhibits BMP ligand activity along the DV axis, which is consistent with the distinct Chordin protein distributions (Figure 10A). In a counter-gradient mechanism Chordin binds BMP in a broader domain, extending over a much greater extent of the DV axis than in a source-sink mechanism (Figure 10B,C). In a source-sink mechanism, Chordin binds BMP largely in dorsal regions and extends little ventrally, effectively generating a driving force for the movement of BMP (Figure 10B,C). This dorsal sink of Chordin leads to a diffusive flux of BMP ligand down its concentration gradient (Figure 10C) that largely shapes the BMP signaling profile in a source-sink mechanism (Figure 10C).
Figure 10
How the source-sink and counter-gradient mechanisms shape the gradient.
(A) The mean BMP and Chd concentrations in all source-sink and counter-gradient solutions fitting WT, chd LOF, and chd heterozygous P-Smad5 data and within a diffusivity between 2.4 and 6.4 µm2/s. (B) The diffusive flux divided by the decay [(DBMP/decBMP)*(d[BMP]/dx)*(1/[BMP]max)] of BMP (blue) with units of 103*µm and rate of binding of BMP to Chd (kon*[BMP]*[Chd]) (red) with units of 3.6*10−2*sec−1 for representitive (B) Counter-Gradient and (C) Source-Sink solutions fitting WT, chd LOF, and chd heterozygous P-Smad5 data and within a diffusivity between 2.4 and 6.4 µm2/s (Figure 8).
Importantly, gradient formation by either of these mechanisms need not be exclusive, and instead characteristics of each can contribute to some extent in shaping sectors of the other’s gradient. In the source-sink simulation in Figure 10A,C, Chordin forms a small counter-gradient partially contributing to the gradient shape in this region. Thus, in some source-sink simulations, Chordin shaped the gradient by simultaneously binding ligand to block signaling in lateral positions and by establishing a sink that serves as a driving force for the diffusive flux of ligand from ventral positions dorsally towards the regions of higher Chordin. Though each solution is classified by the dominant mechanism, many share some aspects of both the source-sink and counter-gradient mechanisms.
BMP diffuses relatively freely
We measured BMP diffusivity for the first time in vertebrates. Using FRAP, we show that BMP can diffuse relatively freely with a diffusivity of 4.4 ± 0.4 µm2/s (Figure 8E), about 4-fold less than secreted Venus diffusion in the zebrafish blastula (≈16 µm2/s) (Figure 8G). Our measured BMP diffusivity is comparable to the diffusivity of Squint (Ndr1, D = 3.2 µm2/s), another TGF-β ligand in the zebrafish blastula that acts as a long-range mesoderm inducer (Müller et al., 2012). This high BMP diffusivity is consistent with previous BMP heterodimer protein injections into the extracellular space of BMP-deficient embryos, which could extend throughout the embryo and restore the WT P-Smad5 gradient within 1.5 hr, suggesting the BMP ligand can move rapidly (Little and Mullins, 2009).
Our computational screen found hundreds of simulations with a BMP diffusivity near 4.4 µm2/s. The vast majority of those simulations were classified as having a source-sink mechanism. The remaining few are classified as having a counter-gradient mechanism. This paucity of counter-gradient simulations is a reflection of the fine-tuning needed for this mechanism to work as compared to the source-sink mechanism. The counter-gradient mechanism requires a specific balance of Chordin diffusivity and decay as well as Tolloid degradation, while the source-sink mechanism does not (Figure 7H,I). Together, this suggests that the source-sink mechanism is more robust to changes in biophysical parameters than the counter-gradient mechanism, but does not rule out the counter-gradient mechanism entirely.
BMP transcriptional feedback: A symptom not a cause
The recent observation that bmp2b and bmp7a are expressed in a graded manner in WT embryos has lead to the hypothesis that the BMP signaling gradient may largely reflect the bmp expression gradient (Ramel and Hill, 2013). We sought to determine if the BMP signaling gradient in zebrafish is predominantly established by matching the bmp expression gradient. Three results do not support this hypothesis. First, the relative shape of the bmp2b expression domain (Figure 3) did not precisely match the P-Smad5 gradient (Figure 2). Second, we mathematically defined a gradient established by a transcriptional mechanism as one where 80% of the BMP accumulates or is degraded where it is produced, as opposed to binding antagonists or diffusing away. When we performed a mathematical model-based computational screen using a bmp expression profile that reflected the P-Smad5 WT gradient, no tested parameter combination fit both our measured WT and chordin LOF P-Smad5 gradients. For the transcriptional mechanism to work, the bmp expression domain would have to change in the chordin mutant condition to fit the mutant gradient profile.
Finally, we showed that feedback by BMP signaling on bmp expression does not begin until after 5.7 hpf. This is likely because the initial bmp expression domain is established by maternal factors and repression from Bozozok, a transcription factor activated by maternal Wnt signaling (Koos and Ho, 1999; Langdon and Mullins, 2011; Leung et al., 2003; Solnica-Krezel et al., 2001). FGF (Fibroblast Growth Factor) and Nodal signaling also repress bmp expression dorsally (Fürthauer et al., 2004; Kuo et al., 2013; Maegawa et al., 2006; Shimizu et al., 2000; Varga et al., 2007). Importantly, BMP signaling does not play a role in the initial establishment of the bmp and chordin expression domains at 4 hpf, as both are unchanged prior to ~6 hpf in BMP pathway mutants (Kishimoto et al., 1997; Miller-Bertoglio et al., 1997; Schmid et al., 2000). At or after ~6 hpf, the bmp expression domain begins to respond to changes in BMP signaling levels, as bmp expression decreases in BMP pathway mutants and increases in BMP antagonist mutants (Kishimoto et al., 1997; Miller-Bertoglio et al., 1999; Nguyen et al., 1998; Ramel and Hill, 2013; Schmid et al., 2000). We show that the bmp2b expression domain does not begin to shift in response to the loss of chordin until after 5.7 hpf (Figure 3), while the P-Smad5 gradient shifts as early as 4.7 hpf (Figure 6).
While we show that the transcriptional mechanism cannot entirely account for the measured P-Smad5 WT, chordin +/-, and chordin -/- profiles, the shape of the transcriptional gradient does still contribute to gradient shape. Four-fold more individual solutions fit the WT, chd +/-, and chd -/- P-Smad5 profiles when the bmp expression input matched our measured bmp2b level (0.4%) than when bmp expression was set to match the P-Smad5 profile (0.1%) (Figure 7E,K). Counter-gradient solutions only appeared when the bmp expression input was broader than the P-Smad5 profile (Figure 7E’,K’). More source-sink solutions were present when the bmp expression input was broader because BMP did not need to travel as far or fast (Figure 7F,G vs. Figure 7L,M) to reach the dorsal sink. Together, these effects show that the shape of the bmp expression domain contributes to BMP gradient formation even if it does not entirely account for it, similar to that found for the bicoid transcript and protein gradients that pattern the Drosophila AP axis (Little et al., 2011).
Integrated approach reveals source-sink mechanism
Although we and most others in the vertebrate field have contended that a Chordin (or BMP antagonist) counter-gradient drives formation of the BMP activity gradient in DV patterning, our studies here, intriguingly, suggest that an alternate source-sink mechanism may prevail. While the source-sink gradient mechanism is also modulated by Chordin, Chordin instead acts in a distinct manner as a sink, binding BMP ligand predominantly in dorsal regions, thus allowing a BMP diffusive gradient to form throughout most of the ventral half of the embryo. Key to deriving this alternate model was the integrated approach used that combined quantitative experimental analysis with computational modeling. Importantly, a role for Chordin in establishing a sink that drives gradient formation would not have been revealed to us had we not performed the mathematical model-based computational screen. By narrowing the compatible simulations successively with the WT, chordin +/-, and then chordin mutant P-Smad5 profiles, 15-fold more source-sink simulations than counter-gradient simulations perdured (Figure 7E’). Furthermore, the computational modeling also illuminated the BMP diffusivity parameter as one to further test the source-sink mechanism. Significantly, our measured Bmp2 diffusivity further supports the source-sink mechanism of gradient formation. Thus the seamless integration of quantitative experimental analysis with mathematical model-based computational screens has proved to be a highly successful approach to elucidating mechanisms of BMP gradient formation. Future studies will be required to definitively determine the mechanism and further test the source-sink and counter gradient models of BMP gradient formation.
Materials and methods
Zebrafish lines
Request a detailed protocolWT TU (RRID:ZIRC_ZL57) zebrafish were used as the wild-type strain. Adult chdtt250/tt250 (RRID:ZFIN_ZDB-GENO-060811-12) or chdtt250/Tu zebrafish were used to produce chd homozygous and heterozygous embryos. Adult chdtt250/tt250 homozygotes were generated by injecting chordin mRNA into 1-cell stage chd-/- embryos to rescue the ventralization phenotype and then were raised to adulthood. Genotyping of adults and embryos was performed using KASPar genotyping (Smith and Maughan, 2015) with primers designed and generated by KBioscience of the following sequences, where [G/A] is the WT/mutant nucleotide:
CTCCTTCGGTGGCCGCTTTTACTCTCTGGAAGACACGTGGCATCCAGATCTCGGAGAGCCGTTTGGTGTGATGCACTGCGTTATGTGTCATTGTGAGCCG[G/A]TGAGTTGTGCACAGTTCAGTTTGAAATCCATATTGAATCTGAATTGACTTCTGCTGCTGAGTTGCAACATTCACACCATATCTAAATTGAATTCATATT
mRNA injection for BMP overexpression
Request a detailed protocolWT zebrafish embryos were injected with mRNA at the 1-cell stage. bmp2b and bmp7a RNA were made using the SP6 MMessage Machine kit (Life Technologies AM1340). bmp2b or bmp7a cDNA in a pBluescript II KS-construct was linearized with NotI. To overexpress BMP, 6 pg of bmp7a RNA and 12 pg of bmp2b RNA were injected. Resulting embryos had a V5 fully ventralized phenotype at 24 hpf (Figure 6H,I).
In situ hybridization and domain size analysis
Request a detailed protocolWhole-mount in situ hybridizations were performed using RNA DIG probes as described using, chordin (Miller-Bertoglio et al., 1997) and noggin1 (Dal-Pra et al., 2006). RNA probes were generated using the Roche DIG RNA labeling kit (11277073910). Embryos were cleared in glycerol, and photographed using a Leica IC80 HD. Images were processed using ImageJ and MATLAB. In situs were stained with Anti-DIG-Alkaline Phosphatase (Roche 11093274910) and developed using BM Purple (Roche Life Sciences).
The sizes of the chordin and noggin expression domains were determined by image processing with MATLAB. Centerpoints of animal views of each embryo were determined by thresholding. The boundaries of the noggin and chordin expression domains relative to the center of the animal view were determined by a second threshold. These points were connected by line segments and the angle was measured (Figure 3A,B).
Fluorescent bmp2b in situ hybridization and image analysis
Request a detailed protocolFluorescence in situ hybridization (FISH) was performed on fixed cryosections using a RNAscope Fluorescent Multiplex Kit (Advanced Cell Diagnostics (ACD)). Embryos were fixed with 4% paraformaldehyde in PBS at 4°C overnight. Embryos equilibrate in 30% sucrose until they sink and incubated in fresh 30% sucrose for 3 days at 4°C. Cryosections (20 μm) at the marginal region were collected on slides, followed by air drying for 30 min at −20°C. In situ hybridization was performed according to the manufacturer’s instructions (ACD). A custom C2 probe was designed for bmp2b (#456471-C2). bmp2b probe was used at 1:10 dilution. Sections were stained for DAPI and images were acquired at 63 × oil objective using a Zeiss 800 upright confocal.
Relative intensity quantification of mRNA levels was performed on maximum intensity projections of 20 μm sections. Marginal cells were grouped into 10 degree intervals along the marginal circumference. Average intensity was quantified in each section using the MATLAB image analysis toolbox. Averaged bmp2b mRNA levels in 2.5 hpf embryos were used to measure background. We found equivalent intensity levels and distributions in the 2.5 hpf embyros and the dorsal bmp2b signal in Wt 5.7 hpf embryo suggesting limited to zero bmp2b expression in the dorsal region. For each cross-section, the right and left side of the distributions were averaged into a single ventral to dorsal profile.
Immunostaining
Request a detailed protocolEmbryos were fixed overnight in 4% paraformaldehyde at 4°C, blocked in NCS-PBST (10% fetal bovine serum, 1% DMSO, 0.1% Tween 20 in PBS), and probed overnight with a 1:100 dilution of anti-phosphoSmad1/5/8 antibody (Cell Signaling Technology, #9511, discontinued), followed by a 1:500 dilution of goat anti-rabbit Alexa Fluor 647-conjugated antibody (Thermo Fisher Scientific, Rockford, IL; Cat# A-21244, RRID:AB_2535812). Embryos were mounted in BABB (benzyl alcohol (Sigma B-1042) and benzyl benzoate (Sigma B-6630), 1:2 ratio) and scanned using a Zeiss LSM 710 confocal microscope with a LD LCI Plan-Achromat 25x/0.8 Imm Corr DIC M27 multi-immersion lens. The oil-immersion setting was used to reduce Mie scattering distortion, spherical aberrations, and chromatic aberrations by minimizing refractive index (R.I.) mismatch between the lens oil (R.I. = 1.518), the coverslip, BABB (R.I.≈1.56), and the light scattering particles in the embryo (R.I.≈1.56). Fluorophore bleaching was greatly reduced by precise embryo orientation, reducing sample thickness, and by high scan speeds using a Zeiss LSM 710 confocal microscope. Nuclei were visualized with Sytox Orange (Molecular Probes) or Sytox Green (Molecular Probes).
Summary of imaging and processing
Request a detailed protocolImmunostained P-Smad5 embryos were processed and imaged as described above. We observed minimal photo-bleaching and spherical aberration (Figure 11A–C). We wrote a Matlab algorithm capable of identifying all 8000 + nuclei centerpoints in each embryo in 3D, removing populations unresponsive to P-Smad5 (such as yolk syncytial nuclei and dividing cells), and extracting the P-Smad5 intensities associated with each nucleus (Figure 11D–L). The resulting individual digital embryos (Figure 2A’–E’) from each condition were averaged together to generate large datasets from which embryo-wide P-Smad5 levels could be quantified in WT and mutant conditions (Figure 6A–B).
Figure 11
Quantifying nuclear P-Smad5 intensities embryo-wide.
(A) Marginal P-Smad5 intensity from a chd LOF embryo imaged twice. (B) Average P-Smad5 intensity drop-off from photo-bleaching of all nuclei in embryos imaged twice (N = 5). (C) There is minimal intensity drop-off due to spherical aberration, as shown by the average intensity of the nuclear DNA stain (Sytox Orange) versus distance from the coverslip (4.7: N = 3, 5.3: N = 4, 5.7: N = 13, 6.3: N = 11, 6.7: N = 4). (D) Maximum projection of an animal view of a single embryo. (E) Nuclei centerpoints (red dots) identified from the sytox nuclear stain (blue). (F) Measured centerpoint nuclear intensities displayed as a heatmap. (G) P-Smad5 is absent in dividing cells (red stain, yellow arrows). Dividing cells have bright condensed chromatin (green stain, yellow arrows). (H) Bright condensed chromatin was used to identify dividing cells. Cells with a 40% elevated DNA stain over the mean (red line) were eliminated from the analysis. (I) Lateral view of a single embryo. (J) Sparse Yolk Syncytial Layer nuclei below the margin are eliminated. (K) Single lateral slice depicting the elimination of remaining yolk syncytial layer nuclei and enveloping layer nuclei by subtracting the outer 15% of all nuclei (filled in circles) to leave only deep cell nuclei (open circles). (L) Lateral view of embryo after outer 15% has been eliminated.
Image processing
Request a detailed protocolNuclear intensities of P-Smad5 were extracted from the stacks of images generated using Matlab algorithms (source code in Supplementary files 1).
The centerpoints of all the nuclei were located using the Sytox DNA stain. The ‘.lsm’ files were converted to ‘.tif’ files using ImageJ, and then imported into Matlab as 1024 × 1024 X Z multidimensional arrays. XY pixels were 0.55 um, Z pixels were 2.3 um. The images were then smoothed using a 9 × 9×3 kernel (most nuclei are 10 × 10×4 pixels large). Local minima and maxima were removed using the ‘imhmax’ and ‘imhmin’ functions. The remaining maxima were found using the ‘imregionalmax’ function on the entire 1024 × 1024 X Z array. Maxima closer together than six pixels were assumed to be in the same nucleus, and were combined. The remaining maxima were assumed to be the centerpoints of the nuclei (Figure 11E).
These centerpoints were used to extract P-Smad5 intensities on the P-Smad5 channel. P-Smad5 distribution in each nucleus was approximately uniform, so a small sphere within each nucleus was averaged to attain the P-Smad5 intensity. On the P-Smad5 channel, pixels within a spherical 6 × 6 × 3 kernel of each maxima were averaged.
Cell types unresponsive to Bmp signaling were removed. P-Smad5 appears to be uniformly distributed throughout the cytoplasm during cell division, making measurement impractical (Figure 11G). In dividing cells, chromatin condenses making DNA stains such as Sytox concentrated and bright. Cells with a bright DNA fluorescence staining above 140% of the mean DNA fluorescence were considered dividing and eliminated from the analysis. Extra-embryonic cells such as the Enveloping Layer (EVL) and the Yolk Syncytial Layer (YSL) did not appear to respond to BMP ligand the same way as the deep cells. These cell types are not patterned along the DV axis by BMP, and were eliminated from our analysis. To do so, sparse EVL cells located below the vegetal margin were eliminated by hand (Figure 11J). Next, the inner and outer layer of approximately 15% of the total cells was eliminated (Figure 11K,L). The remaining cells were assumed to be non-dividing deep cells.
Embryos of similar stages were then aligned and conformed to a template embryo of the same stage using Coherent Point Drift (CPD) (Myronenko and Song, 2010b). Embryos were aligned in the AP direction by fitting them to a sphere, finding a plane that spanned the marginal region, and rotating until that plane was aligned with the XY axis. Embryos were aligned in the DV direction using the embryonic shield as a morphological marker. Before 6 hpf, when the shield is not present, the embryos were aligned in the DV direction by fitting a polynomial regression to the P-Smad5 gradient around the margin and rotating until the max peak was ventral. Next, embryos were all aligned to a template using an affine CPD (Myronenko and Song, 2010b). This corrected for any distortions in embryo shape that may have occurred during fixation and staining.
Embryos from the same set were subjected to no normalization. Embryos stained and imaged on different days with different settings were normalized by multiplying the entire set by a single scalar value. To determine this normalization scalar, control WT embryos were always imaged in conjunction with each experimental condition. The scalar normalization value was determined by minimizing the sum of the error between the control WT embryos imaged on different days.
Generating and comparing P-Smad5 profiles around the margin of the embryo
Request a detailed protocolTo generate marginal profiles, a 40 um thick band of cells around the vegetal margin was chosen for each embryo. Cells within that band were grouped into 10 degree intervals and averaged together to form 36 individual points. The left and right side of the gradient were averaged together into a single ventral to dorsal profile. For 3-D embryo-wide averages, all nuclei were projected onto a sphere fitting the embryo. The sphere was then divided into 4800 approximately equilateral triangles. All nuclei falling within each triangle were averaged together. Slopes were obtained by fitting a lowess fit to the averaged 3-D data’s spherical coordinates phi and theta using the ‘fit’ function in Matlab. To determine if the marginal gradients of WT and chordin mutant embryos were significantly different, two-tailed T-Tests were performed with a rejection of the null hypothesis at the 5% significance level (Figure 6, Figure 6—figure supplement 1A).We observed a difference in WT vs. chordin mutant embryos that was much larger than our observed embryo-to-embryo variability (Figure 6). Our t-tests confirmed that our sample sizes are sufficient to discern differences between the WT and chordin mutant embryos (Figure 6, Figure 6—figure supplement 1A).
Mathematical model-based computational screen method
Request a detailed protocolFor each set of parameters defined in the parameter vector, we solved the five non-linear partial differential equations (PDEs) for BMP ligand, Chordin, Noggin, and the complexes of BMP-Chordin, BMP-Noggin in MATLAB (Figure 4, Figure 4—figure supplement 1A). Equations were solved on the half-circumference, with ‘symmetry’ boundary conditions imposed on the first and last node-point in the spatial discretization. The half-circumference was discretized into 36 node points with equidistant spacing and the second order spatial derivative is discretized via the finite difference method. The production regions of BMP ligand, Chordin, Noggin, are specified along the nodes by mapping the spatial position to subsequent node position (Figures 3 and 4D,E, Table 1). Tld is treated parametrically as a function of position according to its domain of expression (Figure 4D,E,Table 1). Time-stepping of the simulation is handled by the adaptive solver ode15s with a relative tolerance set to 1e-9. The model is solved for the developmental window that spans from 3.5 to 5.7 hpf and all measurements of model error are calculated at 5.7 hpf.
For each model iteration, parameters were selected from a uniform distribution in log space that covered four orders of magnitude within the physiological range for each parameter. Each parameter was selected independently of the remaining parameters. Adaptive and subsampling methods that increase parameter selection in regions with high variance in model output were not used in our parameter selection for a number of reasons. Firstly, a parameter matrix is produced that is then subdivided across distributed computers to solve the PDEs to produce a stored file of model solutions and parameter vectors associated with the stored solution. For each parameter vector, the PDE system is converted into a set of 180 ordinary differential equations after discretization and dynamically solved for WT, chd−/−, nog−/−, and chd+/- conditions using the implicit solver ode15s. Ode15s is well suited to problems with numerical stiffness that arise during numerical screens with random parameters. Thus, for an ensemble of 1 million parameter vectors, the system of differential equations are solved four times to simulate the wt and mutant conditions, increasing the total model evaluations to 4 million. Following calculation of the model solutions for each parameter vector, post processing, sorting, and calculation of model fitness against the data is handled by a separate program that operates on the stored solutions. The separation of model evaluation from model analysis allows for much greater flexibility, the total number of simulation results, and an ability to add additional simulation results to existing simulations that have previously been computed. If a job is interrupted, the index and solutions up to the parameter index are stored and simulation jobs can be restarted at the last iteration point.
Random sampling is not susceptible to the ‘curse of dimensionality’ that is common to regular grid spacing methods allowing for more efficient estimation of the NRMSD landscape in Figure 7(C,D), and random points fill in the parameter space more evenly than regular grid spacing that forces evaluations of the model within the same parameter hyperplane (Caflisch et al., 1998). Figure 7C,D demonstrates the dense coverage and sampling to estimate the NRMSD values and dense coverage in regions with acceptable NRMSD values. Random sampling of the biophysical parameters is common to these types of problems (Eldar et al., 2002; von Dassow et al., 2000).
For each parameter vector, the model is initially solved against WT conditions and subsequent simulations for mutant conditions are carried out for the same parameter vector by setting the corresponding production rates for the mutant to zero and re-simulating. Error between the model results and the fluorescent data are calculated via a two-step process. First, the amplitude of the P-Smad5 fluorescent-intensity data and model peak levels for free BMP are normalized as commonly done when calculating a residual with fluorescent intensity data (Hengenius et al., 2014; Pargett and Umulis, 2013). This approximation is valid considering that (1) BMP ligands are not saturating receptors and (2) Smad5 activity is not saturated (Figure 6G). The scaling parameter determined for model-fitness against P-Smad5 was then applied to the remaining model results to capture any changes in BMP levels in the mutant simulations. Residuals were calculated for WT and mutant conditions independently and simulations were scored for passing the WT and mutant conditions independently as opposed to using an aggregate residual. Simulations were classified as transcriptional, source-sink, counter-gradient, or shuttling.
Constructs, mRNA, and morpholinos (MOs) for FRAP
Request a detailed protocolSequence encoding fluorescent protein Venus was amplified from pBSK12-her1:Ub2-Venus (a gift from Sharon L. Amacher, Ohio State University, OH) (Delaune et al., 2012). This sequence was inserted between the pro- and mature domains of Bmp2b, two amino acids downstream of the proprotein convertase (PC) cleavage site (REKR) with a GSTGTTGGG linker separating the prodomain and the fluorescent protein and a GS linker (GGGGSGGGGS) separating the fluorescent protein from the mature domain. This fusion construct was modified from pCS2(+)-HA-Bmp2b (Little and Mullins, 2009). Sequences encoding Venus protein were also fused to the pro-domain of Bmp2b two amino acids downstream of the proprotein convertase (PC) cleavage site (REKR) with a GSTGTTGGG linker, to generate the secreted-Venus plasmid that lacks the mature Bmp2b ligand domain.
Capped mRNA was synthesized using the mMessage mMachine Kit (Ambion) with SP6 RNA polymerase according to the manufacturer's protocol. Vectors were linearized by digestion with NotI. mRNA encoding the Bmp2b-Venus or secreted-Venus, was injected into one- or eight-cell stage embryos. For rescue experiments, 1 ng of bmp2b MO2 (GTCTGCGTTCCCGTCGTCTCCTAAG) was injected along with 9 pg of bmp2b-venus RNA from a different needle (Imai and Talbot, 2001). To perform FRAP in the absence of Chordin, we injected embryos at the 1 cell stage with 1 ng of chordin MO1 (ATCCACAGCAGCCCCTCCATCATCC) (Nasevicius and Ekker, 2000). Next, bmp2b-venus RNA was injected into a single blastomere at the 8-cell stage. Associated phenotypes are shown in Figure 8B.
Western blot
Request a detailed protocolZebrafish embryos were lysed in Pierce RIPA buffer (89900, Thermo Scientific) supplemented with Halt Protease Inhibitor Cocktail (1862209, Thermo Scientific) and Phosphatase inhibitor Cocktail (1862495, Thermo Scientific). Protein samples mixed with Laemmli sample buffer (Bio-rad) were denatured by incubation for 5 min at 98o, and resolved by SDS-PAGE using Mini-PROTEAN TGX Gels (10%, Bio-rad) and transferred to PVDF membranes (Bio-rad). The membranes were blocked with 5% non-fat milk (Bio-rad) in PBST for 1 hr at room temperature, and incubated with primary antibodies in 2% BSA (Sigma) in PBST at 4°C overnight. After that, the membranes were incubated with HRP-coupled secondary antibodies for 1 hr at room temperature. Chemiluminescence was detected using Clarity Western ECL Substrate (Bio-rad) to obtain the image. Using stripping buffer (46430, Thermo Scientific), the membranes were reused to detect β-Actin as loading controls.
Fluorescence Recovery After Photobleaching (FRAP) mRNA encoding the Bmp2b-Venus fusion protein (50 pg) was injected at the one-cell stage to test the activity of the mRNA in a ventralization assay. Embryos used for FRAP were injected in one cell at the 8-cell stage. Embryos were mounted in 1% low melting point agarose (Sigma) in glass bottom microwell dishes (MatTek Corporation). FRAP experiments were performed using a LSM 800 confocal microscope (Zeiss) with a W Plan-Apochromat 20×/1.0 objective (D = 0.17 M27 75 mm). Photobleaching in a square region (160.4 μm × 160.4 μm) was performed through the depth of the blastoderm with 100% laser power for ~10 min. In the secreted Venus FRAP (Figure 8D), some recovery occurred at the t = 0 time point at the periphery of the bleached region. Recovery of fluorescence was monitored every 10 s in the same imaging plane.
Processing of FRAP data
Request a detailed protocolFrom the 8-cell stage injected embryos, regions lacking Bmp2b-Venus producing cells (visualized by high intensity signal throughout the cytoplasm) were identified. Cells displayed characteristic higher intensity signals in the intercellular space and no signal was detected intracellularly in the non-producing cells. Images are taken before the FRAP experiment commences, and saved every 10 s during acquisition. All files are exported in lossless TIFF format for subsequent quantification in MATLAB. To measure the recovery, all TIFF files are imported into MATLAB and the FRAP region is identified for subsequent measurements. Internal FRAP region is scaled from 8-bit [0 255] to [0, 1] and an extracellular mask is generated by removing background with a minimum threshold level set at 1% of the image maximum value. This excludes the intracellular compartments from biasing the average intensity calculations for the extracellular recovery. With background removed, recovery is calculated as the average fluorescence intensity of the extracellular fluorescence within the masked region.
Calculation of diffusion coefficients from FRAP data
Request a detailed protocolThe FRAP region is modeled using a finite-difference equation for diffusion in the FRAP region. Diffusion was estimated by measuring model recovery starting from zero initial conditions and constant concentration boundary conditions. Masked region dimensions for measurement were mapped directly to node-points and distances in the finite difference model to compare and optimize recovery in the mask region. A steepest-descent optimizer with multiple starts was used to estimate the diffusion coefficient for each FRAP experiment. We did not explicitly model the occlusion and tortuosity of the diffusion process caused by the arrangement of cells, nor did it account explicitly for binding and unbinding to HSPGs and other immobile binding components. The impact of binding and the role of occlusions in the diffusion path are very well known (Cussler, 2009; Garnett, 1904). Therefore our measured diffusion coefficients are the effective diffusion coefficients in zebrafish and not an intrinsic measurement of the diffusion coefficient in a free environment.
Experimental replicate information
All listed replicates in the figure legends are biological replicates.
Whole-mount in situ hybridization
Request a detailed protocolFor chordin and noggin domain size, all biological replicates are reported in Figure 3C, while representative images are shown in Figure 3A–B. All biological replicates were from a single cross stained with different probes. All embryos were reported with no elimination of outliers. For bmp2b expression profile, representative images are displayed for each condition and developmental stage. Biological replicates were taken from two experiments.
Immunohistochemistry
Request a detailed protocolAll biological replicates are indicated and averaged in Figure 2 and Figure 6 and associated legends. Biological replicates from the 4.7 and 5.3 hpf timepoints for the WT and chordin mutant embryos in Figure 2 and Figure 6 were from a single experiment. The 5.7 and 6.3 hpf timepoints were collected from 3 and 4 experiments respectively. The 6.7 hpf (WT-only) timepoint was collected from two experiments. In all cases, only embryos that were clearly damaged during the process were eliminated from the analysis. Mean intensities at each point were an average of a subset of the 8000 + nuclei present in each embryo. The XYZ coordinates and P-Smad5 values of all nuclei in all embryos used is attached as a supplemental file.
FRAP
Request a detailed protocolAll biological replicates are shown in Figure 8E–G. No outliers were eliminated.
Supplemental data files
Request a detailed protocolThe ‘Image Analysis’ folder in Supplementary file 1 contains all files used to identify nuclei, align, and extract P-Smad5 intensities from the ‘.lsm’ files generated by confocal imaging. A Read-Me is included with details on how to use the scripts. A conceptual description can be found above in the materials and methods section. All raw XYZ coordinates and P-Smad5 intensities are included in the ‘P-Smad Intensities’ folder. Each embryo is represented as a five by n array, where columns 1–3 are XYZ coordinates, column four is nuclear stain intensity, and column five is P-Smad5 intensity.
References
-
Generation and interpretation of FGF morphogen gradients in vertebratesCurrent Opinion in Genetics & Development 23:415–422.https://doi.org/10.1016/j.gde.2013.03.002
-
BMP signalling: agony and antagony in the familyTrends in Cell Biology 25:249–264.https://doi.org/10.1016/j.tcb.2014.12.004
-
Computer-aided design of thrombin inhibitorsNews in Physiological Sciences 13:182–189.
-
Morphogen interpretation: the transcriptional logic of neural tube patterningCurrent Opinion in Genetics & Development 23:423–428.https://doi.org/10.1016/j.gde.2013.04.003
-
The role of tolloid/mini fin in dorsoventral pattern formation of the zebrafish embryoDevelopment 126:3119–3130.
-
Diffusion: Mass Transfer in Fluid SystemsDiffusion: Mass Transfer in Fluid Systems, 3rd edition.
-
Noggin1 and Follistatin-like2 function redundantly to Chordin to antagonize BMP activityDevelopmental Biology 298:514–526.https://doi.org/10.1016/j.ydbio.2006.07.002
-
Biphasic activation of the BMP pathway patterns the Drosophila embryonic dorsal regionDevelopment 128:965–972.
-
Patterning the zebrafish axial skeleton requires early chordin functionNature Genetics 23:442–446.https://doi.org/10.1038/70557
-
Colours in metal glasses and in metallic filmsPhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 203:385–420.https://doi.org/10.1098/rsta.1904.0024
-
dino and mercedes, two genes regulating dorsal development in the zebrafish embryoDevelopment 123:95–102.
-
Making models match measurements: model optimization for morphogen patterning networksSeminars in Cell & Developmental Biology 35:109–123.https://doi.org/10.1016/j.semcdb.2014.06.017
-
The molecular nature of zebrafish swirl: BMP2 function is essential during early dorsoventral patterningDevelopment 124:4457–4466.
-
Maternal and zygotic control of zebrafish dorsoventral axial patterningAnnual Review of Genetics 45:357–377.https://doi.org/10.1146/annurev-genet-110410-132517
-
Extracellular modulation of BMP activity in patterning the dorsoventral axisBirth Defects Research Part C: Embryo Today: Reviews 78:224–242.https://doi.org/10.1002/bdrc.20079
-
Maternal and zygotic activity of the zebrafish ogon locus antagonizes BMP signalingDevelopmental Biology 214:72–86.https://doi.org/10.1006/dbio.1999.9384
-
Differential regulation of chordin expression domains in mutant zebrafishDevelopmental Biology 192:537–550.https://doi.org/10.1006/dbio.1997.8788
-
Formation of the BMP activity gradient in the Drosophila embryoDevelopmental Cell 8:915–924.https://doi.org/10.1016/j.devcel.2005.04.009
-
Point set registration: coherent point driftIEEE Transactions on Pattern Analysis and Machine Intelligence 32:2262-75.https://doi.org/10.1109/TPAMI.2010.46
-
Non-rigid point set registration: coherent point driftIEEE Transactions on Pattern Analysis and Machine Intelligence 32:2262–2275.https://doi.org/10.1109/TPAMI.2010.46
-
Effective targeted gene 'knockdown' in zebrafishNature Genetics 26:216–220.https://doi.org/10.1038/79951
-
Shaping BMP morphogen gradients through enzyme-substrate interactionsDevelopmental Cell 21:375–383.https://doi.org/10.1016/j.devcel.2011.06.025
-
Morphogen gradients: from generation to interpretationAnnual Review of Cell and Developmental Biology 27:377–407.https://doi.org/10.1146/annurev-cellbio-092910-154148
-
Temporal dynamics, spatial range, and transcriptional interpretation of the Dorsal morphogen gradientCurrent Opinion in Genetics & Development 22:542–546.https://doi.org/10.1016/j.gde.2012.08.005
-
Evo-devo aspects of classical and molecular data in a historical perspectiveJournal of Experimental Zoology 302:69–91.https://doi.org/10.1002/jez.b.20003
-
Gradients in the brain: the control of the development of form and function in the cerebral cortexCold Spring Harbor Perspectives in Biology 1:a002519.https://doi.org/10.1101/cshperspect.a002519
-
Ectodermal patterning in vertebrate embryosDevelopmental Biology 182:5–20.https://doi.org/10.1006/dbio.1996.8445
-
Dynamics and precision in retinoic acid morphogen gradientsCurrent Opinion in Genetics & Development 22:562–569.https://doi.org/10.1016/j.gde.2012.11.012
-
Equivalent genetic roles for bmp7/snailhouse and bmp2b/swirl in dorsoventral pattern formationDevelopment 127:957–967.
-
Creating gradients by morphogen shuttlingTrends in Genetics 29:339–347.https://doi.org/10.1016/j.tig.2013.01.001
-
SNP genotyping using KASPar assaysMethods in Molecular Biology 1245:243–256.https://doi.org/10.1007/978-1-4939-1966-6_18
-
The role of the homeodomain protein bozozok in zebrafish axis formationThe International Journal of Developmental Biology 45:299–310.
-
Temporally coordinated signals progressively pattern the anteroposterior and dorsoventral body axesSeminars in Cell & Developmental Biology 42:118–133.https://doi.org/10.1016/j.semcdb.2015.06.003
-
Computational analysis of BMP gradients in dorsal-ventral patterning of the zebrafish embryoJournal of Theoretical Biology 248:579–589.https://doi.org/10.1016/j.jtbi.2007.05.026
Decision letter
-
Robb KrumlaufReviewing Editor; Stowers Institute for Medical Research, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Systems biology derived source-sink mechanism of BMP gradient formation" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Robb Krumlauf as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The authors use quantitative imaging and modeling to come to the conclusion that a source-sink mechanism underlies the BMP/Chordin patterning system. This is novel as formation of a BMP gradient by a source-sink mechanism is distinct from previously reported models such as transcriptional, shuttling, or counter-gradient mechanisms. The measurement of BMP2 diffusion in vivo is novel. The topic is not only of broad interest but also quite controversial, and this study brings some much needed clarity to the field. The main weakness is the limited attempt to minimize the parameter space used for the modeling. Some of the measurements are difficult but there are several key measurements that need to be made. The computational and experimental analyses are incomplete, and some details of simulations have not been provided. For their model to be made more convincing a variety of concerns/comments should be fully addressed.
Major specific points:
1) The source-sink model is puzzling in light of previous reports showing that zebrafish and Xenopus chordin mRNA injections at the one-cell stage can fully rescue zebrafish chordino mutants (Fisher and Halpern, 1999; Schulte-Merker et al., 1997). In these embryos, injected chordin mRNA is present throughout the embryo and yet, the rescued mutants develop normally to adulthood (and presumably form a normal BMP gradient). The mRNA rescue experiments also call into question the notion that the transcriptional pattern is critical for BMP gradient formation. The authors should discuss the various models in light of these reports.
2) Relevant to the above, did the authors measure the effective diffusion of BMP2 by FRAP in chordin mutant embryos? This will provide experimental evidence to test the models.
3) Cell numbers increase substantially between 4 and 7 hpf, when the authors report a steep increase in the BMP gradient. How does the change in cell number in the zebrafish gastrula (between 4-7 hpf) compare to that during fly DV patterning? To what extent does this account for the increase in gradient slope in zebrafish?
4) The authors state "The equations were simulated 1,000,000 times".
Table 1 lists 15 free parameters; therefore each parameter can, in principle, only take ~2.5 different values when combined randomly (2.52^15 is ~10^6). However, the plots in Figure 4 suggest that many more values have been chosen for individual parameters. The authors should clarify how the parameters were varied in the simulations.
5) The authors state "To our surprise, there are greater than 50 times the number of source-sink to counter-gradient modeling solutions (Figure 6B'), suggesting that the source-sink mechanism predominates." However, the number of different parameter combinations that can explain the data cannot be taken as a measure of support for a specific modelling solution. Different parameter combinations can result in similar solutions simply because the model is not identifiable or because it is difficult to find a unique solution when trying to fit the data directly.
6) If the authors already have or are making a BMP-destabilized GFP/degradFP fusion, precise measurement of BMP2 decay rates would be useful to distinguish between the models. If these are available it would be useful but not essential for the paper.
7) Why were the equations solved for the developmental window from 3.5 to 5.7 hpf and not from 4.7 to 6.7 hpf, when the pSMAD intensities were measured?
8) The authors presumably solved the equations in 1D in order to sample parameter space in a reasonable amount of time. I think it would add to the paper if the authors showed that the same parameters found to satisfy the 1D case would also set up the 3D P-Smad gradient when solving the equations in the hemispherical geometry of the zebrafish embryo as in Zhang, Lander, and Nie, J. Theor. Biol., 2007.
9) Can the P-Smad5 gradient of the chordin heterozygotes be included? Why are they "not shown"? Does the model fit the chordin heterozygote P-Smad5 gradient?
10) The pSmad5 measurements are impressive but the models ultimately rest on the BMP transcript and (active and inactive) protein expression profiles. I understand that there are no good antibodies to detect BMP but minimally, the authors need to use now standard quantitative fluorescent in situ hybridization approaches to measure BMP transcript distribution.
11) The BMP FRAP experiment goes a long way to reduce parameter space but why not also measure BMP diffusion in the presence of Chordin? This experiment seems trivial but would greatly help to test the shuttling model.
12) The characterization of the bmp2-venus chimera is too superficial. Does it have the same range and activity as wild-type bmp2?
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Systems biology derived source-sink mechanism of BMP gradient formation" for further consideration at eLife. Your revised article has been favorably evaluated by Marianne Bronner as Senior editor, a Reviewing editor, and three reviewers. =
The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below. In the absence of measurements of BMP diffusion in the presence of chordin or Chordin protein expression, the alternate model Zinski et al. suggest for BMP gradient formation cannot be substantiated. To support their clams, it is crucial that the authors:i) determine how the BMP gradient forms in the presence of uniform chordin ii) show expression of tagged BMP / Chordin proteins at early gastrula stages in mutants injected with RNA. iii) clarify how random sampling was done.
Key comments:
1) In response to how BMP diffuses in the absence and presence of chordin (previous review comments 2 and 11), in the revised Zinski m/s, the authors examined BMP diffusion without chordin. However, they did not examine BMP diffusion with uniform chordin. This is a crucial test of the model which is missing.
2) In their explanation of how their data and preferred model can explain the previous reports of rescue of chordino mutants with uniform chordin injections, Zinski et al. contend that chordin RNA injections likely do not generate uniform Chordin protein expression in embryos. But no evidence is provided to support this view – either experimentally or in their simulations. (The triple morphant is not relevant to the comment).
3) Regarding how parameters were varied in their simulations, the authors state that they prefer random sampling as opposed to a regular grid.
Did the authors refine sampling depending on the outcomes? If solutions did not change in a region in parameter space, did they increase the sampling density?
https://doi.org/10.7554/eLife.22199.019Author response
Major specific points:
1) The source-sink model is puzzling in light of previous reports showing that zebrafish and Xenopus chordin mRNA injections at the one-cell stage can fully rescue zebrafish chordino mutants (Fisher and Halpern, 1999; Schulte-Merker et al., 1997). In these embryos, injected chordin mRNA is present throughout the embryo and yet, the rescued mutants develop normally to adulthood (and presumably form a normal BMP gradient). The mRNA rescue experiments also call into question the notion that the transcriptional pattern is critical for BMP gradient formation. The authors should discuss the various models in light of these reports.
The ubiquitous expression of chordin RNA, which can rescue chordin mutant embryos, likely does not generate a ubiquitous distribution of Chordin protein due to the presence of Tolloid and Bmp1a, two ventrally expressed metalloproteases that degrade Chordin protein. The rescue of embryos lacking Chordin, Tolloid, and Bmp1a using chordin RNA injection at the one-cell stage has not been attempted as it is a difficult experiment to execute, but we hypothesize that it would not be effective. We agree that the chordin expression domain is likely not critical to BMP gradient formation.
2) Relevant to the above, did the authors measure the effective diffusion of BMP2 by FRAP in chordin mutant embryos? This will provide experimental evidence to test the models.
We have now repeated our FRAP experiments in Chordin deficient embryos and have added it to the main text of the paper and Figure 8F. The effective diffusivity of Bmp2b-Venus in chd deficient embryos matched the effective diffusivity we had measured in WT embryos. This new experiment validates our previous assumption that the majority of BMP-Venus observed in our WT FRAP experiments was not bound to Chordin, as expected due to the ventralized embryonic phenotype at 24 hpf (Figure 8B row 2). However, performing the FRAP experiment in chd deficient embryos assures that 100% of Bmp2b-Venus observed is not bound to Chordin, removing the potential impact of a subpopulation of Chordin-bound Bmp2b-Venus on the FRAP measurements in WT embryos.
3) Cell numbers increase substantially between 4 and 7 hpf, when the authors report a steep increase in the BMP gradient. How does the change in cell number in the zebrafish gastrula (between 4-7 hpf) compare to that during fly DV patterning? To what extent does this account for the increase in gradient slope in zebrafish?
Fly DV patterning takes place during stage 5. Mitosis is paused during stage 5 until gastrulation begins at stage 6. Cells do not divide at this time, but nuclear localization of P-Mad increases from nearly none to a steep gradient in less than 1 hour from early stage 5 to early stage 6 when gastrulation begins. Thus, in this context while there are no changes in nuclei density, the gradient of PMad starts broad and refines to a higher amplitude and sharper profile [1]. However, in the context of patterning the embryo termini in Drosophila via Torso RTK signaling, the nuclei are present in a syncytium in the absence of cell membranes, where the nuclear density plays a role in shaping the gradient of phosphorylated ERK at the termini [2]. The nuclei density doubles with each successive division providing an increase in nuclear trapping of dpERK at the most terminal positions, resulting in a steeper gradient with higher peak amplitude and reduced gradient extent or length over time. The zebrafish embryo is not a synctium, so this same nuclear trapping mechanism could not occur. Moreover, the gradient DV length is constant over this time period, also inconsistent with this type of mechanism (Figure 2).
From our cell counts in embryos from a single timecourse (Figure 2H), we observe an approximately 70% increase in cell number from 4.7 to 6.7 hours post fertilization (hpf) and a 100% increase in nuclear P-Smad5 amount. Cell nuclei do not change significantly in size during this time [3]. The increase in cell number occurs throughout the embryo and is not restricted to a particular DV region. So half of the 70% increase in cell number occurs in lateral and dorsal regions, where the slope does not change over this time period (Figure 2G). The increase in slope observed during the 2 hour period occurs over the ventral half of the embryo, which accounts for 35% of the increase in cell number, while the gradient peak doubles in amplitude. Thus an increase in cell number, via an unknown mechanism, could not account for the increase in slope. Additional support that cell number has little effect on gradient shape, we observed that the absolute number of cells at a given time point can vary by as much as 20% between different embryos within the same cross or between crosses (Figure 2I) with no detectable change in gradient shape or phenotype.
4) The authors state "The equations were simulated 1,000,000 times".
Table 1 lists 15 free parameters; therefore each parameter can, in principle, only take ~2.5 different values when combined randomly (2.52^15 is ~10^6). However, the plots in Figure 4 suggest that many more values have been chosen for individual parameters. The authors should clarify how the parameters were varied in the simulations.
For each model iteration, parameters were selected from a uniform distribution in log space that covered four orders of magnitude within the physiological range for each parameter. Each parameter was selected independently of the remaining parameters. We prefer the coverage offered by this method, as opposed to a regular grid wherein each parameter has a set of discrete values and coverage at those values is complete. Random sampling is not susceptible to the “curse of dimensionality” that is common to regular grid spacing methods allowing for more efficient estimation of the NRMSD landscape in Figure 7(C,D), and random points fill in the parameter space more evenly than regular grid spacing that forces evaluations of the model within the same parameter hyperplane [4]. Random sampling of the biophysical parameters is common to these types of problems [5, 6].
5) The authors state "To our surprise, there are greater than 50 times the number of source-sink to counter-gradient modeling solutions (Figure 6B'), suggesting that the source-sink mechanism predominates." However, the number of different parameter combinations that can explain the data cannot be taken as a measure of support for a specific modelling solution. Different parameter combinations can result in similar solutions simply because the model is not identifiable or because it is difficult to find a unique solution when trying to fit the data directly.
We agree with this point and have changed the sentence so it now reads: "To our surprise, the source-sink modeling solutions emerged more frequently within our computational screen than the counter-gradient solutions (Figure 7E'). "
6) If the authors already have or are making a BMP-destabilized GFP/degradFP fusion, precise measurement of BMP2 decay rates would be useful to distinguish between the models. If these are available it would be useful but not essential for the paper.
We agree it would be useful, but these constructs are not currently available.
7) Why were the equations solved for the developmental window from 3.5 to 5.7 hpf and not from 4.7 to 6.7 hpf, when the pSMAD intensities were measured?
While dorsal-ventral cell fate for the head and trunk is largely specified from 4 to 7 hours post fertilization [7-10], bmp and chordin are first expressed after the mid-blastula transition (MBT) at 3 hpf [11-14]. We started our simulation at 3.5 hpf to account for the time needed for these proteins to be translated, folded, and secreted. We have now added that rationale to the manuscript.
8) The authors presumably solved the equations in 1D in order to sample parameter space in a reasonable amount of time. I think it would add to the paper if the authors showed that the same parameters found to satisfy the 1D case would also set up the 3D P-Smad gradient when solving the equations in the hemispherical geometry of the zebrafish embryo as in Zhang, Lander, and Nie, J. Theor. Biol., 2007.
We agree that this would add to the paper and the development of data-driven three dimensional models continues to be an effort of our work. We have made progress in developing data-driven 3D models and this is being worked on but not achievable within the timeframe of the revision. The primary scope of the modeling effort herein focused on broad coverage of potential mechanisms to provide an unbiased view of data-consistent mechanisms. At this time, the screen and scope of our study is only achievable with models with one space dimension. Moreover, the inputs to the system based on gene expression are distributed symmetrically across the embryo, so a 1D model should largely reflect one in 3D. We added the rationale for performing a 1D model now to the manuscript.
9) Can the P-Smad5 gradient of the chordin heterozygotes be included? Why are they "not shown"? Does the model fit the chordin heterozygote P-Smad5 gradient?
Thank you for the suggestion. We have now added a comparison of WT to chordin heterozygotes to the manuscript (Figure 6J). There is no significant change in the P-Smad5 gradient or embryonic phenotype between WT and chordin heterozygous embryos (Figure 6J). We also now required mathematical model solutions to fit our chordin heterozygous data, when the Chordin production rate was set to 50%. Forcing the model to fit our chordin heterozygous data decreased the overall number of solutions, but did not eliminate all solutions of either the source-sink or counter-gradient mechanism. The mathematical model in Figure 7–8 and associated text have been updated to include chordin heterozygous data.
10) The pSmad5 measurements are impressive but the models ultimately rest on the BMP transcript and (active and inactive) protein expression profiles. I understand that there are no good antibodies to detect BMP but minimally, the authors need to use now standard quantitative fluorescent in situ hybridization approaches to measure BMP transcript distribution.
bmp2b expression is now quantified by RNAScope fluorescent in situ hybridization and included in Figure 3. We have adjusted our BMP production gradient in our mathematical model to fit the measured bmp2b expression gradient and simulated the mathematical model an additional 1,000,000 times. The adjusted BMP production gradient increased the overall number of solutions, but did not modify our conclusions about the mechanisms. All modeling figures and data displays (Figure 4,5,7,9) have been updated incorporating the new data.
11) The BMP FRAP experiment goes a long way to reduce parameter space but why not also measure BMP diffusion in the presence of Chordin? This experiment seems trivial but would greatly help to test the shuttling model.
While we agree that measuring the diffusivity of BMP bound to Chordin would be very useful for further testing the shuttling model, the experiment is not feasible without extensive overexpression or the development of additional tools. To see Bmp2b-Venus, we must moderately overexpress it. To ensure all Bmp2b-Venus is bound to Chordin, we would need to heavily overexpress Chordin. Chordin has been shown to interact with Tsg [15], Bmper [16-18], Ont1 [19], and HSPGs [20]. By overexpressing both BMP and Chordin, these interacting factors, which may limit BMP-Chordin diffusion could be diluted out, causing FRAP to overestimate the diffusivity of the BMP-Chd complex. Since we do not have time to develop more sophisticated tools to test this, we think it is best not to include a double overexpression experiment.
12) The characterization of the bmp2-venus chimera is too superficial. Does it have the same range and activity as wild-type bmp2?
To assess whether Bmp2b-Venus is properly processed and remains attached to the Venus tag, we performed a Western Blot for Venus on embryos injected with venus or bm2b-venus RNA. We have added this experiment to the paper (Figure 8A).
To further assess the activity and range of the Bmp2b-Venus chimera, we rescued embryos lacking Bmp2b with bmp2b-venus RNA. It has been previously reported that embryos lacking Bmp2b can be rescued using bmp2b RNA [21]. We were able to rescue embryos to a WT phenotype by injecting bmp2b-venus RNA (Figure 8B, rows 5-7). This result has been added to the text (Figure 8). However, a direct comparison of bmp2b RNA activity and bmp2b-Venus activity is not possible, as RNA activity can vary from synthesis reaction to synthesis reaction. In addition to this we have now measured secreted Venus diffusion rates with our approach and added additional measurements for Bmp2b-Venus diffusion in chordin MO-injected embryos (Figure 8).
At this time we are unable to compare the range of Bmp2b to that of Bmp2b-Venus. Range has two components: diffusivity and decay rate. Measuring diffusivity via FRAP requires the protein in question to be tagged fluorescently. While decay rates of Bmp2b and Bmp2b-Venus could be measured and compared by injecting Bmp2b and Bmp2b-Venus protein and performing a western blot time-series, previous estimates of Tgf-β decay rates have relied exclusively on tagged proteins [22].
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Key comments:
1) In response to how BMP diffuses in the absence and presence of chordin (previous review comments 2 and 11), in the revised Zinski m/s, the authors examined BMP diffusion without chordin. However, they did not examine BMP diffusion with uniform chordin. This is a crucial test of the model which is missing.
We presume that by requesting that we measure BMP diffusion with uniform chordin that the reviewers are asking us to measure bound BMP-Chd diffusion. We do assert that the Shuttling mechanism requires high diffusivity of BMP-Chordin (Figure 5D x-axis), however, our chordin-/- mutant data already exclude the Shuttling mechanism (Figure 6). However, when discussing the Source-Sink versus the Counter-Gradient mechanism in Figure 7, we omitted an important plot of BMP-Chd diffusivity in our latest submission, which addresses the reviewers’ point. We thank the reviewers for making us aware of this. A plot of BMP-Chd diffusivity shows that possible values of BMP-Chd diffusivity over the 4 orders of magnitude tested in our simulations does not distinguish between the Counter-Gradient and Source-Sink mechanisms. Both of these mechanisms are compatible with BMP-Chd diffusivities over 4 orders of magnitude from nearly immobile (10-2 µm2/s) to free diffusion (102 µm2/s). Therefore, measuring BMP-Chd diffusivity would not help discern between the two remaining models, Counter-Gradient and Source-Sink. We added this BMP-Chd diffusivity plot now to the manuscript to clarify this important point (Figure 7J).
2) In their explanation of how their data and preferred model can explain the previous reports of rescue of chordino mutants with uniform chordin injections, Zinski et al. contend that chordin RNA injections likely do not generate uniform Chordin protein expression in embryos. But no evidence is provided to support this view – either experimentally or in their simulations. (The triple morphant is not relevant to the comment).
In comment 1 of the previous review, the reviewers asked us to discuss the rescue of chordino mutants with uniform chordin RNA injections. We address this point now by simulating the system with ubiquitous Chordin production. We found that 426 of the 1452 solutions that fit our WT, chordin-/-, and chordin+/- data and BMP diffusivity (within 2 um2/s of our measured 4.4 um2/s) retain a WT BMP gradient when Chordinis uniformly produced. The 426 remaining solutions are all source-sink mechanisms with steep gradients of Chordin that are high dorsally and low ventrally. The results of the simulations show that the BMP gradient is generated by Tolloid degrading Chordin ventrally. We include these simulation results that show how uniform chordin RNA can rescue a chordin mutant in the Results section of our paper, as a new panel Figure 8I. We thank the reviewers for the suggestion, as it provides further support for the source-sink mechanism.
3) Regarding how parameters were varied in their simulations, the authors state that they prefer random sampling as opposed to a regular grid.
Did the authors refine sampling depending on the outcomes? If solutions did not change in a region in parameter space, did they increase the sampling density?
For each model iteration, parameters were selected from a uniform distribution in log space that covered four orders of magnitude within the physiological range for each parameter. Each parameter was selected independently of the remaining parameters. Adaptive and subsampling methods that increase parameter selection in regions with high variance in model output were not used in our parameter selection for a number of reasons. Firstly, a parameter matrix is produced that is then subdivided across distributed computers to solve the pdes to produce a stored file of model solutions and parameter vectors associated with the stored solution. For each parameter vector, the PDE system is converted into a set of 180 ordinary differential equations after discretization and dynamically solved for wt, chd-/-, nog-/-, and chd+/- conditions using the implicit solver ode15s. Ode15s is well suited to problems with numerical stiffness that arise during numerical screens with random parameters. Thus, for an ensemble of 1 million parameter vectors solved, the system of differential equations are solved 4 times to simulate the wt and mutant conditions, increasing the total model evaluations to 4 million. Following calculation of the model solutions for each parameter vector, post processing, sorting, and calculation of model fitness against the data is handled by a separate program that operates on the stored solutions. The separation of model evaluation from model analysis allows for much greater flexibility, the total number of simulation results, and an ability to add additional simulation results to existing simulations that have previously been computed. If a job is interrupted, the index and solutions up to the parameter index are stored and simulation jobs can be restarted at the last iteration point.
Random sampling is not susceptible to the “curse of dimensionality” that is common to regular grid spacing methods allowing for more efficient estimation of the NRMSD landscape in Figure 7(C,D), and random points fill in the parameter space more evenly than regular grid spacing that forces evaluations of the model within the same parameter hyperplane (Caflisch et al., 1998). Figure 7C, D demonstrates the dense coverage and sampling to estimate the NRMSD values and dense coverage in regions with acceptable NRMSD values. Random sampling of the biophysical parameters is common to these types of problems (Eldar et al., 2002; von Dassow et al., 2000).
We have added this explanation to the Materials and methods section as further clarification.
1. Wang, Y. and E.L. Ferguson, Spatial bistability of Dpp–receptor interactions during Drosophila dorsal–ventral patterning. Nature, 2005. 434(7030): p. 225-9.
2. Coppey, M., et al., Nuclear trapping shapes the terminal gradient in the Drosophila embryo. Curr Biol, 2008. 18(12): p. 915-9.
3. Keller, P.J., et al., Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy. Science, 2008. 322(5904): p. 1065-9.
4. Caflisch, A., R. Walchli, and C. Ehrhardt, Computer-Aided Design of Thrombin Inhibitors. News Physiol Sci, 1998. 13: p. 182-189.
5. Eldar, A., et al., Robustness of the BMP morphogen gradient in Drosophila embryonic patterning. Nature, 2002. 419(6904): p. 300-4.
6. von Dassow, G., et al., The segment polarity network is a robust developmental module. Nature, 2000. 406.
7. Tucker, J.A., K.A. Mintzer, and M.C. Mullins, The BMP signaling gradient patterns dorsoventral tissues in a temporally progressive manner along the anteroposterior axis. Dev Cell, 2008. 14(1): p. 108-19.
8. Hashiguchi, M. and M.C. Mullins, Anteroposterior and dorsoventral patterning are coordinated by an identical patterning clock. Development, 2013. 140(9): p. 1970-80.
9. Tuazon, F.B. and M.C. Mullins, Temporally coordinated signals progressively pattern the anteroposterior and dorsoventral body axes. Semin Cell Dev Biol, 2015. 42: p. 118-33.
10. Kwon, H.J., et al., Identification of early requirements for preplacodal ectoderm and sensory organ development. PLoS Genet, 2010. 6(9): p. e1001133.
11. Leung, T., bozozok directly represses bmp2b transcription and mediates the earliest dorsoventral asymmetry of bmp2b expression in zebrafish. Development, 2003. 130(16): p. 3639-3649.
12. Solnica-Krezel, L. and W. Driever, The role of the homeodomain protein Bozozok in zebrafish axis formation. Int J Dev Biol, 2001. 45(1): p. 299-310.
13. Koos, D. and R. Ho, The nieuwkoid dharma homeobox gene is essential for bmp2b repression in the zebrafish pregastrula. Developmental Biology, 1999. 215(2): p. 190–207.
14. Shimizu, T., et al., Cooperative roles of Bozozok/Dhama and Nodal-Related protein in the formation of the dorsal organizer in zebrafish. Mech Dev, 2000. 91(1-2): p. 293-303.
15. Troilo, H., et al., Structural characterization of twisted gastrulation provides insights into opposing functions on the BMP signalling pathway. Matrix Biol, 2016. 55: p. 49-62.
16. Rentzsch, F., et al., Crossveinless 2 is an essential positive feedback regulator of Bmp signaling during zebrafish gastrulation. Development, 2006. 133(5): p. 801-11.
17. Zhang, J.L., et al., Binding between Crossveinless-2 and Chordin von Willebrand factor type C domains promotes BMP signaling by blocking Chordin activity. PLoS One, 2010. 5(9): p. e12846.
18. Ambrosio, A.L., et al., Crossveinless-2 Is a BMP feedback inhibitor that binds Chordin/BMP to regulate Xenopus embryonic patterning. Dev Cell, 2008. 15(2): p. 248-60.
19. Inomata, H., T. Haraguchi, and Y. Sasai, Robust stability of the embryonic axial pattern requires a secreted scaffold for chordin degradation. Cell, 2008. 134(5): p. 854-65.
20. Jasuja, R., et al., Cell-surface heparan sulfate proteoglycans potentiate chordin antagonism of bone morphogenetic protein signaling and are necessary for cellular uptake of chordin. J Biol Chem, 2004. 279(49): p. 51289-97.
21. Nguyen, V.H., et al., Ventral and lateral regions of the zebrafish gastrula, including the neural crest progenitors, are established by a bmp2b/swirl pathway of genes. Dev Biol, 1998. 199(1): p. 93-110.
22. Muller, P., et al., Differential diffusivity of Nodal and Lefty underlies a reaction-diffusion patterning system. Science, 2012. 336(6082): p. 721-4.
https://doi.org/10.7554/eLife.22199.020Article and author information
Author details
Funding
National Institute of General Medical Sciences (R01GM056326)
- Joseph Zinski
- Mary C Mullins
Eunice Kennedy Shriver National Institute of Child Health and Human Development (R01HD073156,T32 HD08318)
- Joseph Zinski
- Ye Bu
- Xu Wang
- Wei Dou
- David Umulis
- Mary C Mullins
National Science Foundation
- Joseph Zinski
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank the UPenn Microscopy Core, especially Andrea Stout and Xinyu Zhao; undergraduate students Samantha Warrick and Dan Zhao; funding from NIH grants T32 HD08318, R01GM056326, R01HD073156, and an NSF Graduate Fellowship; editing help from Francesca Tuazon; and Caroline Hill, Sharon Amacher, Bernard Thisse, and Christine Thisse for reagents.
Ethics
Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to approved institutional animal care and use committee (IACUC) protocols (#803105, #804931) of the University of Pennsylvania Perelman School of Medicine.
Reviewing Editor
- Robb Krumlauf, Stowers Institute for Medical Research, United States
Publication history
- Received: October 12, 2016
- Accepted: August 8, 2017
- Accepted Manuscript published: August 9, 2017 (version 1)
- Version of Record published: September 8, 2017 (version 2)
Copyright
© 2017, Zinski 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,652
- Page views
-
- 642
- Downloads
-
- 58
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
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)
Systems biology derived source-sink mechanism of BMP gradient formation
eLife 6:e22199.
https://doi.org/10.7554/eLife.22199
Further reading
-
- Computational and Systems Biology
- Neuroscience
Inhibition is crucial for brain function, regulating network activity by balancing excitation and implementing gain control. Recent evidence suggests that beyond simply inhibiting excitatory activity, inhibitory neurons can also shape circuit function through disinhibition. While disinhibitory circuit motifs have been implicated in cognitive processes including learning, attentional selection, and input gating, the role of disinhibition is largely unexplored in the study of decision-making. Here, we show that disinhibition provides a simple circuit motif for fast, dynamic control of network state and function. This dynamic control allows a disinhibition-based decision model to reproduce both value normalization and winner-take-all dynamics, the two central features of neurobiological decision-making captured in separate existing models with distinct circuit motifs. In addition, the disinhibition model exhibits flexible attractor dynamics consistent with different forms of persistent activity seen in working memory. Fitting the model to empirical data shows it captures well both the neurophysiological dynamics of value coding and psychometric choice behavior. Furthermore, the biological basis of disinhibition provides a simple mechanism for flexible top-down control of the network states, enabling the circuit to capture diverse task-dependent neural dynamics. These results suggest a biologically plausible unifying mechanism for decision-making and emphasize the importance of local disinhibition in neural processing.
-
- Computational and Systems Biology
Different strains of a microorganism growing in the same environment display a wide variety of growth rates and growth yields. We developed a coarse-grained model to test the hypothesis that different resource allocation strategies, corresponding to different compositions of the proteome, can account for the observed rate-yield variability. The model predictions were verified by means of a database of hundreds of published rate-yield and uptake-secretion phenotypes of Escherichia coli strains grown in standard laboratory conditions. We found a very good quantitative agreement between the range of predicted and observed growth rates, growth yields, and glucose uptake and acetate secretion rates. These results support the hypothesis that resource allocation is a major explanatory factor of the observed variability of growth rates and growth yields across different bacterial strains. An interesting prediction of our model, supported by the experimental data, is that high growth rates are not necessarily accompanied by low growth yields. The resource allocation strategies enabling high-rate, high-yield growth of E. coli lead to a higher saturation of enzymes and ribosomes, and thus to a more efficient utilization of proteomic resources. Our model thus contributes to a fundamental understanding of the quantitative relationship between rate and yield in E. coli and other microorganisms. It may also be useful for the rapid screening of strains in metabolic engineering and synthetic biology.
