Response to Nodal morphogen gradient is determined by the kinetics of target gene induction
Abstract
Morphogen gradients expose cells to different signal concentrations and induce target genes with different ranges of expression. To determine how the Nodal morphogen gradient induces distinct gene expression patterns during zebrafish embryogenesis, we measured the activation dynamics of the signal transducer Smad2 and the expression kinetics of long- and short-range target genes. We found that threshold models based on ligand concentration are insufficient to predict the response of target genes. Instead, morphogen interpretation is shaped by the kinetics of target gene induction: the higher the rate of transcription and the earlier the onset of induction, the greater the spatial range of expression. Thus, the timing and magnitude of target gene expression can be used to modulate the range of expression and diversify the response to morphogen gradients.
https://doi.org/10.7554/eLife.05042.001eLife digest
How a cell can tell where it is in a developing embryo has fascinated scientists for decades. The pioneering computer scientist and mathematical biologist Alan Turing was the first person to coin the term ‘morphogen’ to describe a protein that provides information about locations in the body. A morphogen is released from a group of cells (called the ‘source’) and as it moves away its activity (called the ‘signal’) declines gradually. Cells sense this signal gradient and use it to detect their position with respect to the source. Nodal is an important morphogen and is required to establish the correct identity of cells in the embryo; for example, it helps determine which cells should become a brain or heart or gut cell and so on.
The zebrafish is a widely used model to study animal development, in part because its embryos are transparent; this allows cells and proteins to be easily observed under a microscope. When Nodal acts on cells, another protein called Smad2 becomes activated, moves into the cell's nucleus, and then binds to specific genes. This triggers the expression of these genes, which are first copied into mRNA molecules via a process known as transcription and are then translated into proteins. The protein products of these targeted genes control cell identity and movement.
Several models have been proposed to explain how different concentrations of Nodal switch on the expression of different target genes; that is to say, to explain how a cell interprets the Nodal gradient. Dubrulle et al. have now measured factors that underlie how this gradient is interpreted. Individual cells in zebrafish embryos were tracked under a microscope, and Smad2 activation and gene expression were assessed. Dubrulle et al. found that, in contradiction to previous models, the amount of Nodal present on its own was insufficient to predict the target gene response. Instead, their analysis suggests that the size of each target gene's response depends on its rate of transcription and how quickly it is first expressed in response to Nodal.
These findings of Dubrulle et al. suggest that timing and transcription rate are important in determining the appropriate response to Nodal. Further work will be now needed to find out whether similar mechanisms regulate other processes that rely on the activity of morphogens.
https://doi.org/10.7554/eLife.05042.002Introduction
The Nodal signaling pathway plays essential roles in animal development. Nodal signaling induces and patterns mesendoderm and establishes left-right asymmetry (Conlon et al., 1994; Shen, 2007; Grande and Patel, 2009; Schier, 2009; Duboc et al., 2010; Shiratori and Hamada, 2014). The Nodal signaling pathway regulates dozens of genes, ranging from transcription factors to cytoskeletal components, in order to pattern embryonic tissues (Bennett et al., 2007; Liu et al., 2011; Fodor et al., 2013). In embryonic stem cells, Nodal signaling is required for self-renewal as well as specification of endoderm and mesoderm (James et al., 2005; Vallier et al., 2005; Schier, 2009; Oshimori and Fuchs, 2012; Chen et al., 2013). Nodal signals can form concentration gradients and can act as morphogens (Chen and Schier, 2001; Williams et al., 2004; Müller et al., 2012, 2013; Xu et al., 2014a). It is unclear, however, how different Nodal concentrations induce different target genes and give rise to different cell types.
The classic morphogen threshold model postulates that Nodal signals are secreted from a source and form a concentration gradient that induces different fates in the target tissue according to local ligand concentration (Ashe and Briscoe, 2006; Barkai and Shilo, 2009; Rogers and Schier, 2011). According to this model, high-threshold genes require high levels of Nodal signaling and thus are expressed close to the source (short-range genes), whereas low-threshold genes require lower levels of Nodal and are expressed at a greater distance from the source (long-range genes). Studies of mesendoderm patterning by Nodal in fish and frog have provided five lines of evidence that support the concentration threshold model. First, Nodal signals are produced locally starting at mid-blastula stages, and by the beginning of gastrulation, cells overlapping or close to the Nodal source express endodermal markers, while cells farther away express mesodermal genes (Feldman et al., 1998; Sampath et al., 1998; Gritsman et al., 2000; Chen and Schier, 2001; Harvey and Smith, 2009). Second, a gradient of activated Smad2, the principal transducer of the pathway, peaks at the Nodal source (Faure et al., 2000; Yeo and Whitman, 2001; Harvey and Smith, 2009) with high levels of activated Smad2 in endodermal progenitors and lower levels in mesodermal precursors. Third, reduction of Nodal signaling during blastula stages leads to the absence of endodermal fates but leaves most mesodermal fates intact (Schier et al., 1997; Feldman et al., 1998, 2000; Gritsman et al., 2000; Dougan et al., 2003; Hagos and Dougan, 2007). Fourth, ubiquitous low concentrations of Nodal induce mesodermal markers, whereas high Nodal concentrations induce endodermal markers (Gritsman et al., 2000; Thisse et al., 2000; Dougan et al., 2003). Fifth, an ectopic source of Nodal can induce short- and long-range expression of endodermal and mesodermal markers, respectively (Thisse et al., 2000; Chen and Schier, 2001; Williams et al., 2004; Müller et al., 2012; Xu et al., 2014a). These observations suggest that different concentration thresholds induce different gene expression patterns.
In addition to the contribution of Nodal concentration to target gene induction, the timing of signaling affects Nodal interpretation. For example, the Nodal gradient is not static as signaling activity increases in range and amplitude between the initiation of Nodal expression and the onset of zebrafish gastrulation 2 hr later (Harvey and Smith, 2009; Müller et al., 2012). Moreover, delayed activation or premature inhibition of Nodal activity affects mesendoderm patterning (Gritsman et al., 2000; Dougan et al., 2003; Hagos and Dougan, 2007). Two models have addressed how duration of exposure and changes in concentration contribute to Nodal interpretation. In the snapshot model, cells rapidly adapt their output to the increasing concentration of Nodal, regardless of the duration and history of exposure (Rogers and Schier, 2011). Indeed, increases in activated Smad2 levels are accompanied by an expansion of target gene expression domains (Harvey and Smith, 2009). In this model, the only role of time is to allow the gradient to expand and reach the thresholds that trigger the expression of short- and long-range genes (Harvey and Smith, 2009). The alternative ‘cumulative dose’ or ‘integration’ model postulates that the duration of Nodal signaling plays a critical role in Nodal interpretation. Cells adopt progressively more marginal fates with increasing duration of exposure to Nodal (Gritsman et al., 2000; Hagos and Dougan, 2007). In this model, induction of long-range genes only requires Nodal for short periods of time, whereas activation of short-range genes depends on an extended period of exposure to high Nodal levels. It has therefore been suggested that the total cumulative dose of Nodal signaling determines the cell fate but it is unclear at which level in the pathway a cumulative dose would be measured (Hagos and Dougan, 2007).
Studies of TGFβ signaling in other contexts have suggested additional mechanisms for the time-dependent interpretation of Nodal signaling. In Xenopus, analysis of signaling by Activin, a TGFβ signal related to Nodal, has suggested a ratchet model: the response to the signal is maintained once the ligand has bound the receptor. Indeed, a short pulse of Activin is sufficient to induce and maintain target gene expression several hours after the pulse (Gurdon et al., 1995, 1998; Dyson and Gurdon, 1998; Bourillot et al., 2002). This molecular ‘memory’ has been shown to rely on the persistence of active receptor-ligand complexes (Jullien and Gurdon, 2005) and allows changes in signaling output only in response to increasing Activin concentrations but not to decreasing concentrations.
Cell culture studies have suggested that time-dependent Nodal interpretation is dictated by the dynamics of the signaling pathway (Inman et al., 2002; Xu et al., 2002; Nicolás et al., 2004; Schmierer and Hill, 2005; Guzman-Ayala et al., 2009). TGFβ signaling pathways operate through distinct steps: ligand binding to its receptor, phosphorylation and nuclear accumulation of Smad2, and induction of target gene expression (ten Dijke and Hill, 2004; Massagué, 2012). Several studies have revealed parameters that affect the levels of activated Smad2. For example, cultured human keratinocytes take approximately 60 min of ligand exposure to generate the maximum level of activated Smad2 (Inman et al., 2002). Other studies have shown that the rates of Smad2 phosphorylation and nucleo-cytoplasmic transport affect signaling output (Clarke et al., 2006; Zi and Klipp, 2007; Schmierer et al., 2008; Vizán et al., 2013) or that the speed and frequency of TGFβ ligand presentation influences target gene response (Sorre et al., 2014). These cell culture studies highlight the potential roles of signaling dynamics in target gene induction but it is unclear how these dynamics affect the response to Nodal in vivo.
To distinguish between the numerous proposed mechanisms for Nodal morphogen interpretation, we studied the temporal and spatial dynamics of Smad2 activation and target gene induction in the early zebrafish embryo. We find that not only Nodal concentration and time of exposure but also the kinetics of target gene induction are key determinants of the response to Nodal morphogens. In particular, our study indicates that a target gene's transcription rate and onset of activation are major determinants of expression range, revealing previously unrecognized layers in the interpretation of morphogen gradients.
Results
Smad2 is essential for Nodal signaling
Smad2 activation has been used as a read-out for Nodal signaling, but it has been unclear whether this transcriptional regulator is the main transducer of Nodal signaling in zebrafish (Dick et al., 2000; Jia et al., 2008). To test the role of Smad2 in Nodal signal transduction, we used TILLING (Wienholds et al., 2003) to recover a non-sense mutation in smad2 and generated embryos lacking maternal and zygotic Smad2 (MZsmad2; Figure 1A,B). Endoderm and head and trunk mesoderm are absent in MZsmad2 embryos, a phenotype very similar or identical to Nodal loss-of-function mutants (Figure 1C–E) (Feldman et al., 1998; Gritsman et al., 1999). MZsmad2 mutants could be rescued by ubiquitous expression of wild-type Smad2 and GFP-Smad2 and the larval lethality of Zsmad2 mutants could be rescued to adulthood by a GFP-Smad2 transgene (Figure 1F–H; Table 1). Moreover, neither Nodal nor Activin displayed any activity in MZsmad2 mutants (Figure 1I,J). These results demonstrate that Smad2 is an essential transducer of Nodal signaling during mesendoderm specification.
Spatio-temporal map of Smad2 activity and target gene expression
The nuclear accumulation of GFP-Smad2 is a well-established reporter of TGFβ signaling (Nicolás et al., 2004; Xu and Massagué, 2004; Harvey and Smith, 2009). This approach has been applied in embryos to visualize how the activated Smad2 gradient evolves over time (Harvey and Smith, 2009), but it has not yet been determined how Smad2 activity changes in individual cells and how cell movements might influence gradient interpretation (Xiong et al., 2013) (Figure 2A). We therefore generated stable transgenic lines in which both GFP-Smad2 and histone H2B-RFP were ubiquitously expressed (Figure 2B,C), and tracked GFP-Smad2 nuclear accumulation at the single cell level over time and space.
-
Figure 2—source data 1
- https://doi.org/10.7554/eLife.05042.006
To enable accurate quantification, we determined how Smad2 phosphorylation, GFP-Smad2 phosphorylation and GFP-Smad2 nucleo-cytoplasmic (NC) ratio increased with increasing Nodal signaling (Figure 2—figure supplement 1). These calibrations established that the GFP-Smad2 NC ratio could serve as a read-out for pathway activity and confirmed the graded nuclear accumulation of Smad2-GFP along the vegetal–animal axis (Harvey and Smith, 2009) (Figure 2D, Figure 2—figure supplement 1).
To follow the trajectory of each cell, we tracked individual blastomeres over time (Figure 2D–F, Figure 2—figure supplement 2, Figure 2—source data 1), determined their GFP-Smad2 NC ratio, and measured their distance from the margin. The resulting spatio-temporal map of Smad2 activity revealed that (1) the position of cells relative to the margin did not change extensively until the onset of gastrulation (Figure 2E); (2) cells close to the margin tended to activate Smad2 early and reached the highest levels of activated Smad2; (3) cells located farther away from the margin tended to activate Smad2 with a delay and the levels of activated Smad2 remained low (Figure 2F,G). Thus, during the 1.5 hr from mid- to late-blastula stage a low-amplitude short-range gradient of activated Smad2 is transformed into a high-amplitude long-range gradient.
To determine how the expression range of Nodal target genes correlates with Smad2 activity, we analyzed the expression of the long-range and short-range genes ntl and gsc, respectively (Figure 2H, Figure 2—figure supplements 3, 4). ntl was first faintly detected in a few cells on the presumptive dorsal side of the embryo at the mid-blastula stage. Subsequently, its expression domain intensified and progressively extended animally until the onset of gastrulation. By contrast, gsc expression initiated ∼30 min later and remained confined to a narrow domain on the dorsal side (Figure 2H). Comparing the spatio-temporal maps of Nodal target gene expression and Smad2 activity confirmed that the long-range target gene ntl was induced at both high and low levels of activated Smad2, whereas the expression of the short-range gene gsc correlated with high Smad2 levels and sustained Smad2 activity.
Testing the threshold and ratchet models
The spatiotemporal maps of Smad2 activity and target gene expression are consistent with previous proposals postulating that signaling thresholds determine target gene induction (Harvey and Smith, 2009). To directly test the threshold model of Nodal signaling, we wished to determine whether high Smad2 activity fully predicts the activation of both short- and long-range Nodal target genes. Using transplantation assays, we exposed GFP-Smad2 cells to high Nodal levels for different periods of time and analyzed the relationships between activated Smad2 levels and target gene expression (Figure 3A). GFP-Smad2 NC ratios were similar in cells exposed to Nodal for either 1 or 2 hr (Figure 3B,E). However, while the long-range gene ntl was expressed both after one or 2 hr of exposure to Nodal (Figure 3D,G), the expression of the short-range gene gsc was only detected after 2 hr of exposure (Figure 3C,F). These results are inconsistent with the strictest forms of the threshold model—the level of Smad2 activity at a given time predicts target gene expression—and reveal that the duration of signaling influences morphogen interpretation primarily at the level of target gene induction.
The spatiotemporal maps of Smad2 activity and target gene expression support one prediction of the ratchet model—cells respond to increases in ligand concentration. To test the other tenet of the ratchet model—cells remember the highest ligand concentration they have been exposed to—we determined whether response is refractory to decreasing Nodal levels. We transplanted GFP-Smad2 cells from the blastula margin (where Nodal concentration is high) to the animal pole (where Nodal concentration is low). Inconsistent with the ratchet model, Smad2 activity progressively decreased and reached basal levels after ∼60 min (Figure 3—figure supplement 2A). Similarly, the expression of the long-range gene ntl disappeared over time (Figure 3—figure supplement 2B). Thus, pathway activity and target gene expression cannot be maintained for extended periods after transient exposure to Nodal.
A kinetic model for Nodal morphogen gradient interpretation
Since the threshold and ratchet models do not fully account for Nodal morphogen interpretation, we sought an alternative model based on the biochemistry and biophysics of signaling. The changes in Smad2 activity and gene expression suggested that the kinetics of signal transduction and gene induction might be major factors in Nodal morphogen interpretation. To determine how time and concentration might translate into pathway activity and target gene response, we developed a mathematical description of the kinetics of Nodal signaling (Chen et al., 2010) (Source code 1). To reduce the complexity of the system and the numbers of free parameters, we focused on the three main steps in the pathway (Figure 4A): (1) the diffusion of Nodal from a local source, (2) the Nodal-dependent phosphorylation of Smad2 (pSmad2), and (3) the pSmad2-dependent transcription of target genes. Three coupled differential equations were formulated to implement the kinetic model. All equations were based on standard reaction-diffusion models and mass-action kinetics.
Equation 1 describes the change of Nodal (N) levels over time. Nodal is produced from a source, diffuses and is degraded. Nodal levels at a distance from the source increase with increases in Nodal production (P) and diffusion (D) and with decreases in clearance (k1).
Equation 2 describes the change in activated (phosphorylated) Smad2 (Sp) levels over time. Smad2 activation is proportional to Nodal and non-activated (non-phosphorylated) Smad2 (S) concentrations. Thus, when Nodal concentration increases, activated Smad2 levels increase. Smad2 is deactivated (de-phosphorylated) at rate k3.
Equation 3 describes the induction of Nodal target genes (RNAtarget) over time. For each target gene, levels of expression and dynamics of induction are defined by its maximal transcription rate (α), degradation rate of its RNA (β), and the affinity of pSmad2 for its promoter/enhancer (Kd). The expression of a given target gene increases as α increases, Kd decreases, or the degradation rate decreases. As the concentration of pSmad2 increases, target gene transcription increases. The Hill coefficient n defines the cooperativity that modulates the sensitivity of the response.
Constraining the kinetic model through in vivo measurements
To test the effectiveness of the kinetic model in explaining and predicting Nodal gradient interpretation, we wished to run simulations with a realistic set of parameters. The effective diffusion coefficients and clearance rates of Nodals have been experimentally determined (Müller et al., 2012), but other parameters of the system have not been measured. Exploring the contribution of each of these parameters in regulating target gene expression revealed that multiple parameter combinations could simulate the expression patterns observed in vivo (Figure 4B and Figure 4—figure supplement 1). In particular, the range of expression is affected most dramatically by changes in transcription rate, Kd or Hill coefficient. We therefore decided to constrain the parameter space by performing a detailed quantification of pSmad2 levels and target gene expression at different Nodal concentrations and durations of exposure. To precisely control the levels and timing of ligand input, we injected different amounts of recombinant mouse Nodal protein into the extracellular space of blastula embryos that lacked endogenous Nodal ligands (Figure 5A). Nodal-injected embryos were collected at different time points and pSmad2 levels were determined by Western blotting. Target gene expression levels were measured by NanoString analysis using a custom-designed codeset (Figure 5—source data 1). This technique combines fluorescently barcoded probes with microimaging to detect and count hundreds of transcripts simultaneously in a single hybridization reaction and without amplification. It thus avoids the primer-specific amplification biases of qRT-PCR experiments and allows the direct measurement and comparison of transcript levels (Geiss et al., 2008; Su et al., 2009; Strobl-Mazzulla et al., 2010).
-
Figure 5—source data 1
- https://doi.org/10.7554/eLife.05042.017
-
Figure 5—source data 2
- https://doi.org/10.7554/eLife.05042.018
-
Figure 5—source data 3
- https://doi.org/10.7554/eLife.05042.019
-
Figure 5—source data 4
- https://doi.org/10.7554/eLife.05042.020
-
Figure 5—source data 5
- https://doi.org/10.7554/eLife.05042.021
-
Figure 5—source data 6
- https://doi.org/10.7554/eLife.05042.022
-
Figure 5—source data 7
- https://doi.org/10.7554/eLife.05042.023
Phospho-Smad2
As predicted by the kinetic model and in agreement with previous studies (Zi et al., 2011), increasing Nodal levels induced higher levels of pSmad2 until a plateau was reached. High Nodal levels lead to a very rapid activation of Smad2, reaching a quasi-steady state after 1 hr. In contrast, low Nodal levels induced lower levels of Smad2 activation and quasi-steady state Smad2 activation was only reached 2 hr after ligand exposure (Figure 5B).
We subjected the kinetic model to a fitting procedure to identify values that would best reflect the experimental data (see ‘Materials and methods’). For Smad2 activation, we found phosphorylation rates (range from 2.3 × 10−6 to 4.0 × 10−6 nM−1s−1) and turnover rates (range from 0.9 × 10−4 to 2.7 × 10−4 s−1) similar to previous studies performed in cell culture and Xenopus animal cap cells (Bourillot et al., 2002; Schmierer and Hill, 2005; Schmierer et al., 2008).
Target gene expression
To measure target gene expression, we first identified genes in our NanoString codeset that were directly regulated by Nodal signaling using three criteria: (1) increased expression upon increase of Nodal levels (Figure 5C), (2) Nodal-mediated gene induction in the presence of translation-blocking cycloheximide (Figure 5—figure supplement 1), and (3) binding by Smad2 in the vicinity of transcription start sites (often in conjunction with the co-regulator FoxH1, Figure 5—figure supplement 1, Figure 5—source data 2–5) (Liu et al., 2011; Yoon et al., 2011). This analysis identified 47 direct targets of Nodal signaling.
NanoString analysis allowed precise comparisons of transcript levels in response to different levels and duration of Nodal exposure (Geiss et al., 2008; Strobl-Mazzulla et al., 2010; Nam and Davidson, 2012). Target genes had specific response profiles (Figure 5C, Figure 5—source data 6). For example, ntl, a typical long-range target, was induced at low Nodal concentrations and its expression reached high NanoString counts (∼2500) at high Nodal concentrations. The induction of ntl expression was rapid: ntl mRNA accumulated within 30 min after injection of intermediate levels of Nodal and continuously increased over time. By contrast, gsc, a short-range Nodal target, required higher concentrations of Nodal to be detected above background levels, and its NanoString counts at high Nodal concentrations were 25 times lower than those of ntl. The induction of gsc was slow: mRNA accumulation was only detected after 60 min. These measurements reveal striking differences in the transcriptional magnitude and timing of Nodal-induced gene expression.
To examine whether the kinetic model could capture the behavior of individual target genes, we screened for gene-specific parameter combinations that satisfied the constraints imposed by the NanoString measurements (see ‘Materials and methods’, Source code 2). Parameter value search was limited to defined intervals: 10−3 to 101 counts per second for the transcription rate (Miller et al., 2011; Tu et al., 2014), 0.1 to 100 nM for Kd (Xi et al., 2011; Geertz et al., 2012; Gentsch et al., 2013), 1 × 10−5 to 1 × 10−3 s−1 for the degradation rate (Rabani et al., 2011) and Hill coefficients from 1 to 4. Reflecting the different transcript levels measured by NanoString, transcription rates varied widely between genes, ranging from 0.0016 to 9.3 counts/s (mean 0.54 ± 1.66). In contrast, Kds only ranged from 0.73 to 42 nM, with more than 75% of Kds between 5 and 10 nM (Figure 5—source data 7). For example, we found Kds of Smad2 for ntl and gsc of 4.9 nM and 5.4 nM, respectively, while the transcription rates were 0.67 counts/s for ntl and 0.032 counts/s for gsc. These differences explain why only prolonged exposure to Nodal induced gsc in the test of threshold model (Figure 3): only very low gsc RNA counts (∼50) can be detected 1 hr after Nodal exposure. In contrast, ntl, which was rapidly induced in the test of the threshold model, was induced at high levels (∼1000 counts) after 1 hr. The finding that a realistic set of parameter combinations satisfied the constraints imposed by the NanoString measurements suggests that the kinetic model provides a suitable description of Nodal morphogen interpretation.
The kinetic model predicts the range of target gene expression
To test whether the kinetic model can recapitulate and predict the temporal and spatial pattern of Nodal target gene expression, we ran simulations in a one-dimensional column of cells spanning the vegetal–animal axis. We let Nodal be produced and diffuse from a point source and used the parameters identified in the previous section to simulate the spatial Smad2 activation and transcriptional response over time. Using these conditions, the spatiotemporal pattern of activated Smad2 correlated well with the endogenous pattern of Smad2 activation (Figure 5B): a short-range low-amplitude gradient was transformed over time into a long-range high-amplitude gradient, as observed in vivo for the GFP-Smad2 activity gradient (Figure 2F,G) (Harvey and Smith, 2009).
The simulated spatiotemporal patterns of gene expression also fit well with the in vivo data. For example, in our simulations, ntl expression began in cells close to the margin 45 min after ligand production started, and the range of ntl continuously increased and reached cells located more than 100 µm away from the margin after 3 hr (Figures 2H, 6A). By contrast, gsc expression was delayed and its range of expression was confined to cells close to the source (Figures 2H, 6B).
To test the predictive power of the kinetic model, we determined the expression patterns of genes that had not been analyzed in detail with respect to their range. As predicted by the simulations, foxa3 mRNA rapidly accumulated at high levels up to four cell tiers (∼60 µm) from the source and then extended up to 80–100 µm by the onset of gastrulation (Figure 6C). efnb2a mRNA also readily accumulated but was expressed in a narrower domain, as predicted from the simulations (Figure 6D). These results reveal the power of the kinetic model in recapitulating and predicting the response of target genes to Nodal morphogen signaling.
Transcription rate predicts range of target gene expression
Since the kinetic model predicted target gene expression, we wished to determine which parameters were the major contributors to the range of gene expression. In the simulations described above (Figure 4B and Figure 4—figure supplement 1), genes whose Kd is low and maximal transcription rate is high are expressed at high levels and at long range. In contrast, the degradation rate influences the range of expression only when mRNA half-lives are very short (Figure 4—figure supplement 1). In agreement with the simulations, we found that genes that are highly induced by Nodal generally display a long range of expression (Figure 7A,B). Strikingly, the maximal transcription rate, not the Kd or the degradation rate, was the best predictor of gene expression range (Figure 7C,D). For example, while the degradation rate, Kd and Hill coefficient for the long-range gene foxa3 and the short-range gene gsc are very similar (Figure 5—source data 7), their maximal transcription rate, and therefore their maximal level of expression, differ by a factor of 20. These results raise the possibility that the maximal transcription rate is a major contributor to target gene expression range: the higher the maximal rate of transcription, the longer the range. Moreover, multiple hypotheses analysis indicates that a model in which the maximal transcription rate is gene-specific and Kd is identical for all the genes performs better than a model where Kd is gene-specific and the maximal transcription rate is constant (see Nodal signaling modeling section in ‘Materials and methods’). Although the Kd may affect target gene response, these analyses indicate that the maximal transcription rate is a key parameter in determining the range of expression.
Delayed response to Nodal restricts the range of target gene expression
To analyze additional predictions generated by the kinetic model, we asked whether a delay in gene induction might affect target gene response. To simulate this scenario, we extended the kinetic model with a co-factor that is produced later than and independently of Nodal and acts together with Smad2 to activate gene transcription. Simulations revealed not only the expected delay but also a reduced range of target gene induction: a long-range gene could be transformed into a short-range gene by introducing a delay in gene induction (Figure 8A).
To determine whether such delayed genes might exist in vivo, we screened our NanoString data for Nodal targets whose induction upon Nodal exposure was delayed (Figure 8B, Figure 8—figure supplement 1). We discovered a small set of genes that were induced slowly after Nodal exposure at 3.5 hpf but more rapidly after exposure at 4.5 hpf (Figure 8B, Figure 8—figure supplement 1). For example, when Nodal was injected at 3.5 hpf, efnb2b was only induced after approximately 2 hr. By contrast, when Nodal was injected at 4.5 hpf, the delay in efnb2b induction was reduced by more than 30 min (Figure 8B). In contrast, most other genes responded rapidly to Nodal exposure at either time point (Figure 5—figure supplement 1). This result revealed that the delay was gene-specific and did not reflect a general lack of competence to respond to Nodal signaling or activate gene expression. Similar to the canonical target genes, genes with delayed induction contained Smad2/FoxH1 binding sites and showed a clear response upon injection of Nodal (Figure 8—figure supplement 1) but their induction was abolished in the presence of cycloheximide (Figure 8—figure supplement 1). These Nodal target genes are therefore likely to be regulated not only by Nodal signaling but additional factors. Strikingly, and as predicted by the delayed induction model, efnb2b expression in the embryo was detected only late (6 hpf) and at a short range (5 cell tiers) (Figure 8C). Similarly, the Nodal target gene bra (Martin and Kimelman, 2008) could only be induced shortly before gastrulation and, as predicted by the model, was expressed at low levels and at a short range (Figure 8—figure supplement 2). These results reveal that a delay in transcriptional response can be used to limit the range of morphogen-induced gene expression.
Discussion
Numerous models have been proposed to explain how morphogen gradients are interpreted to generate diverse gene expression patterns. To interrogate these models, we have taken a quantitative approach to measure the parameters that underlie gradient formation (Müller et al., 2012) and interpretation (this study). This approach reveals that the kinetics of target gene induction is a major determinant of morphogen interpretation and suggests that a kinetic model of morphogen interpretation is better suited for the Nodal morphogen system than the prominent threshold and ratchet models.
The kinetic model recapitulates the dynamics of Smad2 activation and reveals how distinct gene expression patterns can be generated: (1) the Nodal morphogen gradient forms and extends through diffusion; (2) rapid phosphorylation generates a corresponding gradient of activated Smad2; target genes are induced based on (3) their affinity for activated Smad2, (4) their maximal transcription rate, and (5) their competence to respond to activated Smad2. Thus, a target gene can be induced rapidly and at a long range by high transcription rate, high Smad2 affinity and early onset of induction. Conversely, low affinity for Smad2, low transcription rate or late onset of induction generate short-range gene expression patterns.
Our analysis identifies transcription rates and induction delays as two novel strategies to modulate morphogen interpretation. Previous models of morphogen interpretation have emphasized the importance of differential DNA (or chromatin) affinity: the higher the affinity for the transcription regulator, the longer the range of target gene expression. Our results do not contradict such models but reveal that in a rapidly developing system, the intrinsic rate of transcription of a target gene can be a major determinant of gene expression range: high affinity binding sites cannot overcome the limits imposed on gene expression range by low levels of intrinsic transcription. Similarly, delays in transcriptional onset can turn high affinity target genes into short-range genes.
The Nodal morphogen system stands in contrast to two other well-studied morphogen systems, Sonic Hedgehog (Shh) and Bicoid. The Shh gradient patterns the dorsoventral axis of the mammalian neural tube over several days (Briscoe and Ericson, 1999; Nishi et al., 2009; Cohen et al., 2013). Although different concentrations of Shh elicit different transcriptional responses, feedback and cross-regulatory interactions are key determinants of patterning. For example, cells are progressively desensitized to Shh activity by upregulating patched1 (Dessaud et al., 2007), and downstream targets regulate each other to generate discrete domains of expression (Balaskas et al., 2012). Thus, in contrast to the Nodal system which rapidly establishes target gene expression patterns, the Shh system makes extensive use of feedback inhibition and cross-regulation. At the other extreme, the Bicoid morphogen has already formed a quasi steady-state gradient before its target genes can be activated during zygotic genome activation (Driever and Nüsslein-Volhard, 1988; Gregor et al., 2007; Porcher et al., 2010). Bicoid concentration and affinity to regulatory chromatin elements are important (but not the sole [Ochoa-Espinosa et al., 2009; Chen et al., 2012]) determinants of target gene expression along the anterior-posterior axis of the Drosophila embryo (Driever et al., 1989; Burz et al., 1998). Thus, in contrast to the Nodal morphogen gradient, which evolves and is interpreted continuously under pre-steady state conditions, the Bicoid morphogen system makes only limited use of temporal strategies to modulate target gene response.
The influence of transcription rates and delays in morphogen interpretation raises the question how these processes might be regulated at the molecular level. Transcriptional delay might be achieved by a co-activator for target gene induction. Alternatively, a repressor might have to be eliminated for a target gene to become competent to respond. Transcription rates might be influenced by local chromatin structure, promoter strength, and by co-activators that boost or repressors that dampen the levels of target gene expression (Li et al., 2007; Lupien et al., 2008; Hager et al., 2009; Kanodia et al., 2012; Peterson et al., 2012; Coulon et al., 2013; Oosterveen et al., 2013; Foo et al., 2014; Xu et al., 2014b). In either case, our study suggests that the intensity and onset of target gene transcription can be major determinants in shaping morphogen gradient interpretation. Similar mechanisms might modulate other rapid and dynamic pattern formation processes (Bolouri and Davidson, 2003; Lewis, 2003; de-Leon and Davidson, 2010; Oates et al., 2012).
Materials and methods
Fish strains and transgenics
Fish were raised and maintained under standard conditions. Wild-type embryos were collected from TLAB in-crosses. MZoeptz57 embryos were obtained as previously described (Zhang et al., 1998; Gritsman et al., 1999). Mutations in the smad2 gene (ENSDARG0000006389, zv9) were screened for in the sperm of ENU-treated males by TILLING with primers encompassing exons 9 and 10.
Outer primers
Request a detailed protocolOSm2F1: 5′-CAATGGAGATAAGCCTGTGGC.
OSm2R4: 5′-TCTGCAAATGTTTTAAGCACTATTTCAG.
Inner primers
Request a detailed protocolISm2F2: 5′-CTGTAATTTTAATATATTCATTTTTGCTGGC.
ISm2R3: 5′-TGCATAAATCTAATTGGCATTTTTAGATAAACC.
A stop mutation was selected in exon 9 (arg335 > Stop) and heterozygous fish were produced by in vitro fertilization. After at least three outcrosses of smad2vu99/+ fish with a TLAB strain, smad2vu99/smad2vu99 germline carrier fish were obtained by the germline transplantation technique as described (Ciruna et al., 2002). The carrier fish were in-crossed to produce MZsmad2 embryos. gfp-smad2 and h2b-rfp transgenic strains were generated by the meganuclease I-SceI technique as described (Thermes et al., 2002). egfp was introduced in frame upstream of the zebrafish smad2 coding sequence. Both gfp-smad2 and h2b-rfp were inserted downstream of a 5.3 kb fragment of the zebrafish β-actin promoter (gift of F Maderspacher) and the cassettes were subcloned into the I-SceI vector. GFP-Smad2 and H2B-RFP transgenic fish were intercrossed to generate double trans-heterozygotes.
Embryo manipulations
mRNA injections
Request a detailed protocolConstructs were cloned in pCS2+ and mRNA for smad2, gfp-smad2, squint and activin were synthesized using the mMessage mMachine kit (Ambion, Grand Island, NY). cyclops and squint morpholinos (MOs) (Gene Tools, Philomath, OR) were previously described (Feldman and Stemple, 2001; Karlen and Rebagliati, 2001). Dechorionated embryos were injected at the one-cell stage with 0.5–1 nl of mRNA at the appropriate concentration.
Transplantations
Request a detailed protocolCell transplantations were performed by mouth pipetting. 20–50 cells were transplanted from a donor (gfp-smad2/h2b-rfp or wild-type) to a host embryo (wild-type or squint injected MZoep). To test Nodal signaling maintenance after input removal, marginal cells of gfp-smad2/h2b-rfp or wild-type donor embryos at 30% epiboly were transplanted into the animal pole of stage-matched wild-type host embryos. Transplanted embryos were further incubated and processed for either confocal microscopy or for in situ hybridization. To analyze the relationships between activated Smad2 levels and Nodal target gene expression, animal pole cells of gfp-smad2/h2b-rfp embryos at 3.5 hpf or 4.5 hpf were transplanted into the animal pole of squint injected stage-matched MZoep embryos. Under these conditions, only donor cells in transplanted embryos can respond to Nodal. Transplanted embryos were incubated for 1 or 2 hr and processed for confocal microscopy and for in situ hybridization.
Nodal induction dynamics
Request a detailed protocolWild-type embryos were first injected at the one-cell stage with 1 nl of cyclops and squint MOs mixture (at 0.2 mM and 4 µg/µl, respectively) to inhibit endogenous Nodal signals. Morphants were further injected with 0.5–1.5 nl of recombinant mouse Nodal protein (rmNodal, R&D Systems, Minneapolis, MN) at different concentrations in the extracellular space at 4 hpf (sphere stage). To test transcriptional competence, rmNODAL was injected at different stages ranging from 3.25 to 5.25 hpf. Cycloheximide (Sigma, Saint Louis, MO) was applied to dechorionated embryos in embryo medium at 50 µg/µl at 3 hpf. 30 min later embryos were injected with 100 nM Nodal protein and were incubated for 90 min before being processed for total RNA extraction and NanoString analysis.
Embryo samples processing
In situ hybridization
Request a detailed protocolIn situ hybridization on whole mount embryos were carried out using standard protocols. ntl, gsc, foxa3, efnb2a, efnb2b and bra probes were previously described (Schulte-Merker et al., 1992, 1994; Bennett et al., 2007; Martin and Kimelman, 2008). A ntl-gsc fusion probe was constructed by amplifying the full length gsc mRNA using the primers cg [ggatcc] ATGCCCGCTGGGATGTTTAGTATC and ataagaat [gcggccgc] TTAGATATTACTTTAATATTTGTTCCTGTTTTCAGGC and cloning into a plasmid containing full length ntl mRNA using BamHI and NotI enzymes. This construct was linearized with NotI and transcribed using T3. The mRNA was injected into embryos collected from a TLAB incross at the 1-cell stage at four different doses (2 pg, 8 pg, 25 pg, and 100 pg). Embryos were cultured to the 128–256-cell stage and fixed in 4% formaldehyde overnight. They were dehydrated, rehydrated, and stained using standard methods. The same concentrated probe stocks were used that had been used in Figure 2, which were freshly diluted 1:100 in hybridization medium. Following staining, embryos were cleared in benzyl benzoate/benzyl alcohol (2:1 vol/vol) Five random embryos were chosen for each probe–dose combination and imaged from the animal cap in one session with no changes to the settings of the microscope. These images were blinded (so that the file names no longer reflected the treatment), converted to grayscale, and the region of the animal cap that represented a single layer of stained cells was selected in ImageJ with the Lasso tool. The mean brightness of this region and a similarly sized background region were calculated. Each image file was quantitated three separate times and averaged, and the background brightness was subtracted. Brightness measurements were inverted (so that more staining would be a higher number), and their mean and standard deviation was plotted.
Western blotting
Request a detailed protocolWestern blots were performed using standard procedures and the signal was detected by chemiluminescence (ECL plus, Amersham, Piscataway, NJ). Phosphorylated Smad2 was probed using a rabbit anti-phospho-Smad2 (Ser465/467) antibody (1:2000 dilution, Cell Signaling Technology, Danvers, MA, #3104) and total Smad2 was probed using a rabbit anti-Smad2/3 antibody (1:2000 dilution, Cell Signaling Technology, #3102). Signals were quantified in ImageJ with the Gel Analyzer function and the ratio between phospho-Smad2 and total Smad2 signals was calculated.
NanoString
Request a detailed protocolTotal RNA from 5 to 10 embryos for each data point was extracted using the RNAeasy mini kit (Qiagen, The Netherlands) and 100 ng of input RNA was processed through the nCounter assay using standard protocols (NanoString Technologies, Seattle, WA) (Kulkarni, 2011). Samples were first normalized to positive controls included in the codeset. The codeset content was further normalized to 11 reference genes to correct for difference in sample input between assays, according to manufacturer's guidelines.
qPCR
Request a detailed protocolTotal RNA from 5 to 10 embryos was extracted using the RNAeasy mini kit (Qiagen) and 100 ng of input RNA was used to synthesize cDNA with the iScript kit (Bio-Rad, Hercules, CA). qPCR reactions were performed in duplicates using the Go Taq qPCR kit (Amersham) on a MX3000P qPCR instrument (Agilent Technologies, Santa Clara, CA). Relative expression of a given gene was calculated by the ∆Ct procedure using eef1a1l1 as a reference. Primers used for qPCR analysis are as followed:
bra
F: 5′ CTGTAGGGAACTCCTCTCAGT
R: 5′ AAGCAGCTGTGTCGTATAAAG
eef1a1l1
F: 5′ AGAAGGAAGCCGCTGAGATGG
R: 5′ TCCGTTCTTGGAGATACCAGCC
flh
F: 5′ GGCGGAGATGAGAGAACGAAC
R: 5′ GATAGCAGAACACGGGATAGC
gsc
F: 5′ GAGACGACACCGAACCATTT
R: 5′ CCTCTGACGACGACCTTTTC
Time lapse imaging, image processing, analysis and cell tracking
Request a detailed protocolLive embryos were embedded in 0.8% low melting point agarose on a glass bottom culture dish (MatTek, Ashland, MA), with the marginal region facing the objective. The dish was filled with fish water (Instant Ocean sea salt [0.6 g/l] in RO water, 0.01 mg/l methylene blue) to prevent dehydration. Images were acquired on a PASCAL confocal microscope (Zeiss, Germany) using a 25× objective (LCI Plan-Neofluar/0.8) equipped with a heated stage set at 28°C. Samples were simultaneously excited with an argon laser at 488 nm and a Helium laser at 546 nm. Four confocal planes were imaged at 3 µm intervals (512 × 512 size, 12-bit depth, line averaging eight times) every 3 min for a period of 3 hr. Embryonic position of the recorded field was assessed morphologically at the end of the imaging session. Image stacks were processed using custom-made Matlab scripts to measure centroid localization of nucleus, nucleo-cytoplasmic ratio of GFP-Smad2 intensity, distance from the margin, and cell tracks (see Supplementary file 1). Channels of stacked confocal images were split in ImageJ and saved as grayscale TIFF image sequences (8-bit for H2B-RFP, 16-bit for GFP-Smad2). H2B-RFP images were further converted to binary images, by applying a threshold using Otsu's method. Objects smaller than 20 pixels were then removed, and the resulting images were segmented using the Moore-Neighbor tracing algorithm modified by Jacob's stopping criteria. The centroid location and area of each nucleus were then extracted. The binary image of the H2B-RFP was used as a mask on the corresponding GFP-Smad2 image to extract nuclear only- and cytoplasmic only GFP-Smad2 signals. The ratio between the mean nuclear GFP intensity and the mean cytoplasmic GFP intensity was used to define Smad2 activity at the single cell level. In MZoep mutants (Nodal insensitive), the mean NC ratio value is 1.19 ± 0.07. Cell tracking was perfomed using the nearest-neighbor strategy based on the centroid position of each nucleus at different time frames. The nearest centroid of the next frame was selected as being part of the cell track if it was less than 10 pixels apart. This process was reiterated through all the frames to generate cell tracks. Based on visual checks of the resulting tracks, ∼90% of the tracks are estimated to be accurate. The distance between each centroid and the margin was measured at each time point. The position of the margin was defined using a user interface: the maximal projection of the H2B-RFP channel was displayed and six reference points were manually selected along the yolk-blastoderm boundary. The whole margin position was then extrapolated by fitting a polynomial curve. The fitted function was used to determine the distance of each centroid from the margin.
Smad2/FoxH1 chromatin immuno-precipitation
Request a detailed protocolEmbryos for Smad2 and FoxH1 ChIP were collected at dome stage after 5 pg squint mRNA injection or after treatment with the Nodal signaling inhibitor SB505124 (Sigma S4696) at 20 µM final. For FoxH1 ChIP, embryos were injected with 5 pg of FoxH1-flag mRNA at 1-cell stage, and anti-flag antibody was used for the pull down.
For each ChIP, 800 embryos were collected and fixed in 1.85% formaldehyde for 15 min at 20°C. Formaldehyde was quenched by adding glycine to a final concentration of 0.125 M. Embryos were rinsed three times in ice-cold PBS, and resuspended in cell lysis buffer (10 mM Tris-HCl pH7.5/10 mM NaCl/0.5% NP40) and lysed for 15 min on ice. Nuclei were collected by centrifugation, resuspended in nuclei lysis buffer (50 mM Tris-HCl pH 7.5/10 mM EDTA/1% SDS) and lysed for 10 min on ice. Samples were diluted three times in IP dilution buffer (16.7 mM Tris-HCl pH 7.5/167 mM NaCl/1.2 mM EDTA/0.01% SDS) and sonicated to obtain fragments of ∼500 bp. Triton X-100 was added to a final concentration of 0.75% and the lysate was incubated overnight while rotating at 4°C with 25 µl of protein G magnetic Dynabeads (Invitrogen) pre-bound to an excess amount of antibody. Antibodies used were anti-FLAG M1 (Sigma F3165), anti-Smad2/3 (Invitrogen, Grand Island, NY 51–1300). Bound complexes were washed six times with RIPA (50 mM HEPES pH7.6/1 mM EDTA/0.7% DOC/1% Igepal/0.5 M LiCl) and TBS and then eluted from the beads with elution buffer (50 mM NaHCO3/1% SDS). Crosslinks were reversed overnight at 65°C and DNA purified by the QIAquick PCR purification kit (Qiagen). Libraries were prepared according to the Illumina sequencing library preparation protocol and sequenced on an Illumina HiSeq 2000. ChIP-seq reads were mapped to the zebrafish genome (UCSC Zv9 assembly) and peaks were called using MACS (Zhang et al., 2008).
Nodal signaling modeling
The goal of modeling Nodal signaling is to predict the range of expression of Nodal target genes in the embryonic blastula and to analyze the key parameters regulating gene response.
Equations 1–3 were used to model the kinetics of Nodal signaling.
P: Production rate of Nodal from the source, where when x ≤ 25 µm and P = 0 when x > 25 µm.
DN: Diffusion coefficient of Nodal.
k1: clearance rate of Nodal.
k2: activation (phosphorylation) rate of Smad2.
k3: de-activation (de-phosphorylation) rate of Smad2.
α: maximal transcription rate of Nodal target gene.
β: degradation rate of Nodal target gene.
Kd: effective dissociation constant of activated Smad2 for target gene enhancer.
n: Hill coefficient.
We assume that the pool of total Smad2 remains constant such that Stotal = S + Sp.
We used two different scenarios to reflect the experimental set up: the ‘homogenous’ scenario, where ectopic Nodal ligand is injected uniformly into a Nodal depleted embryo, and the ‘spatial gradient’ scenario, where Nodal is produced locally on one side of a one-dimension column of cells.
Parameterization
Request a detailed protocolKnown parameters: DN: 1.5 µm2/s; k1: 1 × 10−4 s−1 (Müller et al., 2012); Estimated parameters: Smad2total: 25 nM/cell (Schmierer et al., 2008); Unknown parameters: γ, k2, k3, α, β, Kd, n.
Homogenous model
Request a detailed protocolIn this model, we assume that the exogenous Nodal concentration is uniformly distributed in the embryo and reaches steady-state shortly after injection. In this case,
Initial conditions:
Unknown parameters were then retrieved through minimization of the residual sum of squared errors for the fitted model using the Nelder-Mead simplex method (using a constrained version of the MATLAB function fminsearch) using three different initial guesses spanning the parameter space. The best set of parameters was selected according to the highest coefficient of determination. Parameter confidence intervals of 95% were computed from the residuals and the coefficient covariance.
Spatial gradient model
Request a detailed protocolIn this model, we consider the behavior of a column of cells spanning the vegetal–animal axis of the embryo (500 µm in length) during a 3-hr time span. Parameter values for Smad2 activation and gene induction were taken from the homogenous model. The unknown parameter left in this spatial model is γ, the maximal production of Nodal. As expected for a morphogen molecule, the system is very sensitive to the levels of Nodal. We thus manually set γ to a value where the simulated expression pattern of the well-characterized Nodal target ntl fits the in vivo distribution (γ = 0.03 nM/s). Each gene was then individually simulated using its specific parameters, and its range of expression was defined as the distance from the source where the expression drops below 100 counts.
Initial conditions:
Boundary conditions
Request a detailed protocolSimulations were solved numerically using the MATLAB pdepe function.
Delay model
Request a detailed protocolIn order to model the delay in gene induction, we introduced a co-factor Y required to activate target gene transcription in cooperation with Smad2. In this case,
and
with k4 = 4.6 × 10−5 nM/s and with initial conditions at T0 = 0, Y0 = 0. To abolish delay, we solved (3) with initial conditions where at T0, Y0 = Yfinal.
Model comparison and complexity
Request a detailed protocolMost models of morphogen signaling and interpretation use large numbers of parameters and can thus suffer from overfitting. We thus considered six different models with different numbers of parameters to describe Smad2-dependent transcription of Nodal target genes and compared the probability that these models can generate the data.
The same rate of Smad2 activation is shared among these models:
where Sp, N and S are phosphorylated Smad2, Nodal and non-phosphorylated Smad2 concentrations, respectively, with k1 = 3.1 × 10−6 nM−1s−1 and k2 = 1.8 × 10−4 s−1.
All models for RNA production have two terms: a pSmad2-dependent mRNA transcription rate, and a linear mRNA degradation rate.
In Model 1, we assume that the effective transcription rate is linearly proportional to pSmad2 concentration:
In Model 2, we assume that the transcription rate is regulated by a dissociation constant for pSmad2:
Model 3 is similar to Model 2, with the addition of a maximum transcription rate coefficient:
Model 4 is similar to Model 2, with the addition of a Hill coefficient and a fixed transcription rate coefficient.
A = 0.1 count/s, corresponding to the mode value of the maximal transcription rate coefficient distribution of our fully developed model (see Table 1).
Model 5 is similar to Model 4, except that in this case, the maximal transcription rate coefficient is let free while the dissociation constant is fixed:
C = 6.7 nM, corresponding to the mode value of the dissociation constant distribution of our fully developed model (see Table 1).
Finally, Model 6 is the fully developed model:
Assuming that our NanoString measurements {yi} are noisy with a standard deviation of {σi}, we can consider yi as a Gaussian random variable with a mean value f(ti;θ) of the underlying model containing a vector of parameters θ and a variance σi2 (Bialek, 2012). We thus have,
and
The probability of the data given the underlying model is
Given
Therefore minimizing χ2 by fitting the parameters θ increases the probability that the model could have produced the data. However, different classes of models with different numbers of parameters whose values are unknown have to be considered. To determine the probability of the data given a class of models with unknown K parameters, an integration over all the possible values of the parameters, weighted by some prior knowledge, has to be computed
where P(θ) is the probability of the a priori distribution of the parameters.
χ2 is proportional to N, and we can write
where
We use a saddle point approximation such that
where θ* is the value at which f(θ) is minimized, and H is the Hessian matrix of the second derivatives at this point. Taking the negative log probability of the data given the model class, we have
Since H is a K × K matrix, , and we finally have
The negative log probability measures the length of the shortest code for the data being generated given the class of models. This length depends on the sample size of the data, the number of parameters, the quality of the fit, and some prior on the parameters that we consider flat in our case. Therefore, the model giving the smallest value of the code is to be considered the best model explaining the data given the sample size. The NanoString data and associated noise (which is ∼10% of the count value based on the analysis of the positive spikes across a cartridge) are identical in all our models, so we are left to compare
Calculating the mean value across all genes, we found
CM1 = 601.1
CM2 = 767.6
CM3 = 681.0
CM4 = 264.3
CM5 = 134.5
CM6 = 204.1
CM5 < CM6 < CM4 < CM1 < CM3 < CM2.
Thus, among all the six different models we considered, model M5 and M6 are the most probable models given the data, highlighting the importance of the maximal transcription rate.
References
-
The interpretation of morphogen gradientsDevelopment 133:385–394.https://doi.org/10.1242/dev.02238
-
Robust generation and decoding of morphogen gradientsCold Spring Harbor Perspectives in Biology 1:a001990.https://doi.org/10.1101/cshperspect.a001990
-
Nodal signaling activates differentiation genes during zebrafish gastrulationDevelopmental Biology 304:525–540.https://doi.org/10.1016/j.ydbio.2007.01.012
-
Transcriptional regulatory cascades in development: initial rates, not steady state, determine network kineticsProceedings of the National Academy of Sciences of USA 100:9371–9376.https://doi.org/10.1073/pnas.1533293100
-
A changing morphogen gradient is interpreted by continuous transduction flowDevelopment 129:2167–2180.
-
The specification of neuronal identity by graded Sonic Hedgehog signallingSeminars in Cell & Developmental Biology 10:353–362.https://doi.org/10.1006/scdb.1999.0295
-
Classic and contemporary approaches to modeling biochemical reactionsGenes & Development 24:1861–1875.https://doi.org/10.1101/gad.1945410
-
Production of maternal-zygotic mutant zebrafish by germ-line replacementProceedings of the National Academy of Sciences of USA 99:14919–14924.https://doi.org/10.1073/pnas.222459999
-
Systems theory of Smad signallingSystems Biology 153:412–424.https://doi.org/10.1049/ip-syb:20050055
-
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
-
Eukaryotic transcriptional dynamics: from single molecules to cell populationsNature Reviews Genetics 14:572–584.https://doi.org/10.1038/nrg3484
-
A primary requirement for nodal in the formation and maintenance of the primitive streak in the mouseDevelopment 120:1919–1928.
-
Endogenous patterns of TGFbeta superfamily signaling during early Xenopus developmentDevelopment 127:2917–2931.
-
Massively parallel measurements of molecular interaction kinetics on a microfluidic platformProceedings of the National Academy of Sciences of USA 109:16540–16545.https://doi.org/10.1073/pnas.1206011109
-
Direct multiplexed measurement of gene expression with color-coded probe pairsNature biotechnology 26:317–325.https://doi.org/10.1038/nbt1385
-
Time-dependent patterning of the mesoderm and endoderm by Nodal signals in zebrafishBMC Developmental Biology 7:22.https://doi.org/10.1186/1471-213X-7-22
-
smad2 and smad3 are required for mesendoderm induction by transforming growth factor-beta/nodal signals in zebrafishThe Journal of Biological Chemistry 283:2418–2426.https://doi.org/10.1074/jbc.M707578200
-
Pattern formation by graded and uniform signals in the early Drosophila embryoBiophysical Journal 102:427–433.https://doi.org/10.1016/j.bpj.2011.12.042
-
Digital multiplexed gene expression analysis using the NanoString nCounter systemCurrent Protocols in Molecular Biology, Chapter 25, Unit25B.10, 10.1002/0471142727.mb25b10s94.
-
Information processing at the foxa node of the sea urchin endomesoderm specification networkProceedings of the National Academy of Sciences of USA 107:10103–10108.https://doi.org/10.1073/pnas.1004824107
-
Global identification of SMAD2 target genes reveals a role for multiple co-regulatory factors in zebrafish early gastrulasThe Journal of Biological Chemistry 286:28520–28532.https://doi.org/10.1074/jbc.M111.236307
-
TGFβ signalling in contextNature Reviews Molecular Cell Biology 13:616–630.https://doi.org/10.1038/nrm3434
-
Dynamic transcriptome analysis measures rates of mRNA synthesis and decay in yeastMolecular Systems Biology 7:458.https://doi.org/10.1038/msb.2010.112
-
Analysis of Smad nucleocytoplasmic shuttling in living cellsJournal of Cell Science 117:4113–4125.https://doi.org/10.1242/jcs.01289
-
Modeling the spatio-temporal network that drives patterning in the vertebrate central nervous systemBiochimica Et Biophysica Acta 1789:299–305.https://doi.org/10.1016/j.bbagrm.2009.01.002
-
Anterior-posterior positional information in the absence of a strong Bicoid gradientProceedings of the National Academy of Sciences of USA 106:3823–3828.https://doi.org/10.1073/pnas.0807878105
-
SoxB1-driven transcriptional network underlies neural-specific interpretation of morphogen signalsProceedings of the National Academy of Sciences of USA 110:7330–7335.https://doi.org/10.1073/pnas.1220010110
-
The Harmonies Played by TGF-? in stem cell biologyCell Stem Cell 11:751–764.https://doi.org/10.1016/j.stem.2012.11.001
-
Morphogen gradients: from generation to interpretationAnnual Review of Cell and Developmental Biology 27:377–407.https://doi.org/10.1146/annurev-cellbio-092910-154148
-
Spatial reconstruction of single-cell gene expressionNature Biotechnology, 10.1038/nbt.3192.
-
Nodal morphogensCold Spring Harbor Perspectives in Biology 1:a003459.https://doi.org/10.1101/cshperspect.a003459
-
The one-eyed pinhead gene functions in mesoderm and endoderm formation in zebrafish and interacts with no tailDevelopment 124:327–342.
-
Mathematical modeling identifies Smad nucleocytoplasmic shuttling as a dynamic signal-interpreting systemProceedings of the National Academy of Sciences of USA 105:6608–6613.https://doi.org/10.1073/pnas.0710134105
-
Expression of zebrafish goosecoid and no tail gene products in wild-type and mutant no tail embryosDevelopment 120:843–852.
-
The protein product of the zebrafish homologue of the mouse T gene is expressed in nuclei of the germ ring and the notochord of the early embryoDevelopment 116:1021–1032.
-
Nodal signaling: developmental roles and regulationDevelopment 134:1023–1034.https://doi.org/10.1242/dev.000166
-
TGFβ signaling in establishing left-right asymmetrySeminars in Cell & Developmental Biology 32:80–84.https://doi.org/10.1016/j.semcdb.2014.03.029
-
Histone demethylase JmjD2A regulates neural crest specificationDevelopmental Cell 19:460–468.https://doi.org/10.1016/j.devcel.2010.08.009
-
New insights into TGF-beta-Smad signallingTrends in Biochemical Sciences 29:265–273.https://doi.org/10.1016/j.tibs.2004.03.008
-
I-SceI meganuclease mediates highly efficient transgenesis in fishMechanisms of Development 118:91–98.https://doi.org/10.1016/S0925-4773(02)00218-6
-
Quantitative developmental transcriptomes of the sea urchin Strongylocentrotus purpuratusDevelopmental Biology 385:160–167.https://doi.org/10.1016/j.ydbio.2013.11.019
-
Activin/Nodal and FGF pathways cooperate to maintain pluripotency of human embryonic stem cellsJournal of Cell Science 118:4495–4509.https://doi.org/10.1242/jcs.02553
-
Efficient target-selected mutagenesis in zebrafishGenome Research 13:2700–2707.https://doi.org/10.1101/gr.1725103
-
Nucleocytoplasmic shuttling of signal transducersNature Reviews Molecular Cell Biology 5:209–219.https://doi.org/10.1038/nrm1331
-
HEB and E2A function as SMAD/FOXH1 cofactorsGenes & Development 25:1654–1661.https://doi.org/10.1101/gad.16800511
-
Model-based analysis of ChIP-Seq (MACS)Genome Biology 9:R137.https://doi.org/10.1186/gb-2008-9-9-r137
Article and author information
Author details
Funding
National Institutes of Health (NIH)
- Alexander F Schier
- Lilianna Solnica-Krezel
The funder had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank Bérengère Dubrulle-Bréon for suggestions on modeling, Sharad Ramanathan for help with cell tracking algorithms, modeling and statistical analysis, Will Talbot, Susan Mango, Sharad Ramanathan, Nathan Lord and members of the Schier lab for critical reading of the manuscript. This work was supported by NIH grants to LS-K and AFS.
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 (#25-08) of Harvard University.
Copyright
© 2015, Dubrulle 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
-
- 4,586
- views
-
- 907
- downloads
-
- 94
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Developmental Biology
Mammalian sperm delve into the female reproductive tract to fertilize the female gamete. The available information about how sperm regulate their motility during the final journey to the fertilization site is extremely limited. In this work, we investigated the structural and functional changes in the sperm flagellum after acrosomal exocytosis (AE) and during the interaction with the eggs. The evidence demonstrates that the double helix actin network surrounding the mitochondrial sheath of the midpiece undergoes structural changes prior to the motility cessation. This structural modification is accompanied by a decrease in diameter of the midpiece and is driven by intracellular calcium changes that occur concomitant with a reorganization of the actin helicoidal cortex. Midpiece contraction occurs in a subset of cells that undergo AE, and live-cell imaging during in vitro fertilization showed that the midpiece contraction is required for motility cessation after fusion is initiated. These findings provide the first evidence of the F-actin network’s role in regulating sperm motility, adapting its function to meet specific cellular requirements during fertilization, and highlighting the broader significance of understanding sperm motility.
-
- Developmental Biology
- Genetics and Genomics
O-GlcNAcylation is an essential intracellular protein modification mediated by O-GlcNAc transferase (OGT) and O-GlcNAcase (OGA). Recently, missense mutations in OGT have been linked to intellectual disability, indicating that this modification is important for the development and functioning of the nervous system. However, the processes that are most sensitive to perturbations in O-GlcNAcylation remain to be identified. Here, we uncover quantifiable phenotypes in the fruit fly Drosophila melanogaster carrying a patient-derived OGT mutation in the catalytic domain. Hypo-O-GlcNAcylation leads to defects in synaptogenesis and reduced sleep stability. Both these phenotypes can be partially rescued by genetically or chemically targeting OGA, suggesting that a balance of OGT/OGA activity is required for normal neuronal development and function.