Optogenetic investigation of BMP target gene expression diversity

Signaling molecules activate distinct patterns of gene expression to coordinate embryogenesis, but how spatiotemporal expression diversity is generated is an open question. In zebrafish, a BMP signaling gradient patterns the dorsal-ventral axis. We systematically identified target genes responding to BMP and found that they have diverse spatiotemporal expression patterns. Transcriptional responses to optogenetically delivered high- and low-amplitude BMP signaling pulses indicate that spatiotemporal expression is not fully defined by different BMP signaling activation thresholds. Additionally, we observed negligible correlations between spatiotemporal expression and transcription kinetics for the majority of analyzed genes in response to BMP signaling pulses. In contrast, spatial differences between BMP target genes largely collapsed when FGF and Nodal signaling were inhibited. Our results suggest that, similar to other patterning systems, combinatorial signaling is likely to be a major driver of spatial diversity in BMP-dependent gene expression in zebrafish.


Introduction
Embryogenesis is orchestrated by signaling pathways that activate spatiotemporally diverse patterns of gene expression. A prominent theory relating signaling to gene expression diversity is the gradient threshold model, in which a signaling gradient across a tissue defines unique spatial gene expression domains by activating target genes at different signaling thresholds ( Figure 1A; Sharpe, 2019;Briscoe and Small, 2015;Dubrulle et al., 2015;Rogers and Schier, 2011;Barkai and Shilo, 2009;Ashe and Briscoe, 2006). Gene expression patterns can also be influenced by signaling dynamics and expression kinetics (Sagner and Briscoe, 2017) as well as interactions with other signaling pathways (Briscoe and Small, 2015). However, in many patterning systems the factors leading to diverse developmental gene expression profiles are incompletely characterized. Here, we investigate how signaling levels, target gene expression kinetics, and combinatorial signaling contribute to gene expression diversity during dorsal-ventral patterning in zebrafish.
We focused on patterning mediated by BMP, a TGF-b superfamily member with important developmental roles across the animal kingdom (reviewed in Zinski et al., 2018). BMP ligands bind and assemble complexes of type I and II receptor serine/threonine kinases, resulting in the phosphorylation of the signal transducers Smad1/5/9 and activation of BMP target genes ( Figure 1B; Derynck and Budi, 2019). The regulation of BMP gradient formation during early development has been analyzed in a variety of organisms including Drosophila, Nematostella, and Xenopus (Genikhovich et al., 2015;Iber and Gaglia, 2007;Mizutani et al., 2005;Plouhinec et al., 2013) as well as zebrafish. During late blastula and early gastrulation stages in zebrafish embryos, graded transcription and subsequent diffusion of BMP ligands, together with dorsally secreted BMP inhibitors such as Chordin, generate a ventrally-peaking gradient of BMP signaling that patterns the dorsal-ventral axis ( Figure 1C; Pomreinke et al., 2017;Zinski et al., 2017). Loss of BMP signaling results in dorsalization, whereas excess BMP signaling produces ventralized embryos (Zinski et al., 2018). The degree of dorsalization or ventralization can be modulated by mutations in BMP pathway components with different strengths (Mintzer et al., 2001;Barth et al., 1999;Nguyen et al., 1998;Mullins et al., 1996) or by injecting different amounts of mRNA encoding pathway activators or inhibitors (Schumacher et al., 2011;Dick et al., 2000;Kishimoto et al., 1997;Neave et al., 1997). These observations have led to the suggestion that BMP functions as a morphogen to pattern the dorsal-ventral axis by activating different target genes at different signaling level thresholds ( Figure 1A; Zinski et al., 2018;Tuazon and Mullins, 2015;Schumacher et al., 2011;Barth et al., 1999;Nguyen et al., 1998;Neave et al., 1997;Mullins et al., 1996). However, overexpression and genetic manipulations can affect the duration of signal exposure, dysregulate other signaling pathways, and modify earlier aspects of development such as morphogenetic movements, complicating the interpretation of these experiments. Moreover, patterning of the dorsal-ventral axis by BMP and the germ layers by FGF and Nodal occurs simultaneously in zebrafish (Zinski et al., 2018), and although these pathways are known to interact, how FGF and Nodal influence the spatiotemporal expression of BMP target genes has not been systematically assessed.
To identify the factors that contribute to differences in BMP target gene expression and rule out factors that do not contribute, we first identified BMP targets in early zebrafish embryos and quantified their diverse spatial ( Figure 1) and temporal ( Figure 2) expression patterns. We then used an optogenetic approach to generate acute BMP signaling pulses ( Figure 3) and found that while most target genes can respond to early BMP signaling (Figure 4), differential transcription kinetics do not fully account for the observed expression differences ( Figure 5). Further, target gene responses to high-and low-amplitude signaling pulses suggest that not all spatiotemporal target gene expression differences are due to different signaling activation thresholds ( Figure 6). In contrast, inhibition of FGF and Nodal signaling homogenized the spatial expression patterns of BMP targets, suggesting that combinatorial regulation by BMP, FGF, and Nodal is a major driver of BMP target gene spatial diversity (Figure 7).
According to the gradient threshold model, target genes are activated by distinct signaling levels, leading to different spatial domains of target gene expression in the presence of a signaling gradient ( Figure 1A). To determine whether the BMP patterning system fits this paradigm, we first sought to characterize both BMP signaling distribution and spatial target gene expression. We assessed spatial BMP signaling from 2.75 to 7.25 hpf (256-cell stage -60% epiboly) using immunofluorescent Embryos were collected every 30 min from 2.75 to 7.25 hpf, and transcript levels were quantified using NanoString technology. Error bars represent standard error. Temporal profiles were fit with the cumulative distribution function of the normal distribution (gray lines), and activation time (arrowheads) was defined as the average time point at which the curves reached about two mean average deviations (i.e., 1:5 Á t ) from the inflection point n (excluding the maternally deposited genes id2a [Chong et al., 2005] and smad6a [White et al., 2017]). NanoString probes for two high-confidence activated BMP target genes (apoc1l and znfl2b) were not functional. (O) Average gene expression spatial range is plotted against average activation time. See the Figure 2-source data 1 file for source data. The online version of this article includes the following source data for figure 2: Source data 1. Source data for Figure 2. stainings to detect the BMP signal transducer pSmad1/5/9. We imaged embryos using selective plane illumination microscopy (SPIM) and quantified fluorescence along the dorsal-ventral axis (Figure 1-figure supplement 1E, Materials and methods). Similar to previous studies (Pomreinke et al., 2017;Zinski et al., 2017;Ramel and Hill, 2013;Tucker et al., 2008), we observed a ventrally-peaking BMP signaling gradient that increases in amplitude over time ( Figure 1D,E). We then used fluorescence in situ hybridization and SPIM to quantify the spatial expression profiles of BMP target genes along the dorsal-ventral axis at shield stage (~6.75 hpf) and found that target genes have different expression profiles along this axis ( Figure 1F-Y, Figure 1-figure supplement 1E, Materials and methods; some genes could not be quantified due to weak signal or inability to reliably identify the ventral side). The shape of the expression profiles can be well described by bell curves. We therefore used regression analysis with the Gaussian function and defined the range of each target gene as Using this definition, spatial gene expression broadness ranges from 40-100% dorsal-ventral embryo length ( Figure 1F-Y). Strikingly, pronounced differences along the orthogonal animal-vegetal axis were also evident: Genes were either uniformly expressed along this axis on the ventral side (sizzled, ved, apoc1l, and bambia), restricted to the margin (cdx4 and eve1), or excluded from the margin (foxi1, klf2b, gata2a, and tfap2c) ( Figure 1C,F-O'). Margin exclusion resulted in distinct dorsal-ventral profiles in which mRNA levels peak around 30% embryo length ( Figure 1R,T,W,X), compared to non-excluded genes that peaked more ventrally ( Figure 1P,Q,S,U,V,Y). Therefore, some of the spatial diversity in BMP target gene expression arises from differences along the animal-vegetal axis, orthogonal to the dorsal-ventral BMP signaling gradient.
The gradient threshold paradigm ( Figure 1A) implies that genes with broad ranges should be activated by lower signaling levels. Since signaling levels increase over time (Figure 1D,E;Pomreinke et al., 2017;Zinski et al., 2017;Ramel and Hill, 2013;Tucker et al., 2008), we sought to determine whether more broadly expressed targets were activated earlier. To assess temporal expression of BMP targets, we used NanoString molecular barcoding (Kulkarni, 2011) to measure transcript levels from 2.75 to 7.25 hpf (256-cell stage -60% epiboly) (Figure 2A-N). The shape of the temporal expression profiles can be well approximated by the modified cumulative distribution function of the normal distribution We used this function for regression analysis of the temporal expression profiles and defined activation times as the average time point at which the curves reached about two mean average deviations (i.e., 1:5 Á t ) from the inflection point n. BMP target gene activation times defined in this way ranged from 3.4 to 6.3 hpf ( Figure 2). The gradient threshold model predicts a monotonic decrease when comparing range and activation time. While this relationship is not observed for the entire dataset ( Figure 2O), there is a decreasing monotonic trend when foxi1, eve1, and cdx4 are excluded (note that in contrast to the other genes, the expression of eve1 and cdx4 was only quantified in the embryonic margin Materials and methods]). This suggests the possibility that subsets of BMP target genes may behave consistently with the gradient threshold model. We therefore sought to investigate the relationship between BMP signaling and target gene expression further using an optogenetic strategy.
Reversible optogenetic activation of BMP signaling in vivo using Opto-BMP Opto-BMP at the one-cell stage and their uninjected siblings were reared in the dark or exposed to blue light for 10 h starting 70-80 min post-fertilization. Ventralization phenotypes V1-V4 (indicating excess BMP signaling) were scored at 1 day post-fertilization. Number of embryos: uninjected dark = 59, Opto-BMP dark = 53, uninjected light = 55, Opto-BMP light = 60. (C-E) Uninjected and Opto-BMP-injected embryos were exposed to blue light (2300 lux) for 30 min starting at high stage (3.5 hpf) or shield stage (6.75 hpf) and fixed during and after exposure. pSmad1/5/9 immunofluorescence was quantified and plotted in (C) as Opto-BMP minus uninjected signal with piecewise linear interpolation between timepoints; Figure 3 continued on next page blue light-homodimerizable LOV domain (Rogers and Müller, 2020;Takahashi et al., 2007) and targeted the fusions to the membrane using a myristoylation motif ( Figure 3A), similar to previous approaches (Ramachandran et al., 2018;Vopalensky et al., 2018;Sako et al., 2016). Blue light (~450 nm) exposure should lead to dimerization of the LOV domains and interaction of the BMP kinase domains, activating BMP signaling ( Figure 3A). Injection of mRNA encoding Opto-BMP into zebrafish embryos at the one-cell stage resulted in strong ventralization in light-reared embryos, consistent with excess BMP signaling, whereas darkreared siblings were mostly aphenotypic ( Figure 3B and To facilitate optogenetic experiments, we developed a light exposure device by embedding blue LEDs into the lid of a standard six-well plate and controlling light intensity and dynamics with a single-board computer (Figure 3-figure supplement 1B-D, Materials and methods). Using the LED array, we exposed uninjected and Opto-BMP-injected embryos to blue light for 30 min during high (3.5-4 hpf) or shield (6.75-7.25 hpf) stages, fixed embryos during and after exposure, and quantified BMP signaling using pSmad1/5/9 immunofluorescence (Figure 3C-E and Figure 3-figure supplement 1I,J). At both stages, Opto-BMP embryos showed a dramatic increase in BMP signaling within 10 min of light exposure, and signaling levels returned to normal after light removal. These experiments demonstrate that Opto-BMP reversibly activates BMP signaling in zebrafish embryos in response to light.
Most BMP target genes are competent to respond to BMP at early stages BMP target genes are activated over a range of developmental stages, from 3.4 to 6.3 hpf ( Figure 2). Time-dependent differences in competence -a gene's ability to respond to signaling -may underlie the diversity in activation timing ( Figure 4A). To test this, we quantified BMP target gene expression in uninjected and Opto-BMP-injected embryos exposed to 30 min blue light during either high (3.5-4 hpf) or shield stage (6.75-7.25 hpf) (Figure 4 and In response to a strong BMP signaling pulse at high or shield stage ( Figure 3C-E and Figure 3figure supplement 1I,J), we observed corresponding significant pulses of BMP target gene expression for all genes except crapb2b and cdx4 ( Figure 4K,L). While cdx4 is not competent to respond to an early BMP signaling pulse and crabp2b did not clearly respond to either an early or late signaling pulse, all other tested high-confidence BMP target genes responded at high stage. Therefore, differences in competence to respond to BMP signaling at early stages do not explain the majority of diversity in activation timing.

Transcription kinetics in response to BMP do not fully explain spatiotemporal expression
Target gene transcription kinetics can play important roles in defining spatial expression domains. For example, it has been suggested that Nodal target genes with faster transcript accumulation rates have broader spatial expression domains (Dubrulle et al., 2015). To investigate how the transcription kinetics of BMP target genes may influence their spatiotemporal expression patterns, we assessed the dynamics of target gene responses ( Figure 4) to  . We reasoned that the early activation timing and broad spatial range of some BMP targets might be explained by more rapid transcription in response to BMP. In this paradigm, early BMP signaling activates expression of all target genes at the same time, but transcripts of more slowly transcribed genes only accumulate to detectable levels at later stages, causing them to appear to be 'late-activated' ( Figure 5A). Similarly, broader spatial ranges could be caused by faster accumulation of rapidly produced transcripts that would therefore be detectable farther from the ventral side than more slowly produced transcripts.
To determine whether higher transcript accumulation rates correlate with broader spatial ranges or earlier activation times, we first assessed maximum transcript counts in response to BMP signaling pulses at high or shield stage ( Figure 4). Assuming similar transcript degradation kinetics, transcripts with faster production rates should accumulate to higher levels in response to a BMP signaling pulse ( Figure 5A). However, we observed a weak negative correlation ( Figure 5B) or no correlation ( Figure 5D) between maximum transcript counts and activation time, and found similar results for range ( Figure 5C,E). This suggests that differences in transcript accumulation rates in response to BMP do not fully account for differences in activation timing and spatial broadness. We then used a second approach to assess transcript accumulation kinetics that does not require the assumption of similar transcript degradation rates ( Figure 5F-Z). We fitted the transcription data from the shield-stage BMP signaling pulse with a model involving the known pSmad1/5/9 input ( Figure 3C, Figure 5L, and Figure 3-figure supplement 1J) and parameters reflecting transcript induction (s) and decay (l) ( Figure 5M-Z, Materials and methods). Each of the three experimental repeats was fitted individually, and average s and l values were calculated for each gene. We found a weak negative correlation between s and activation time ( Figure 5F), and no correlation between s and range ( Figure 5G). We also observed a weak negative correlation between l and activation time ( Figure 5H), and no obvious correlation between l and spatial broadness ( Figure 5I). These results are consistent with the maximum transcript count analysis ( Figure 5B-E) and with an alternative fitting approach (Figure 5-figure supplement 1G-W, Materials and methods). In addition, we observed a strong positive correlation between maximum transcript count and s ( Figure 5J), and no correlation between maximum transcript count and l ( Figure 5K), suggesting that production dominates transcription kinetics, and supporting the use of maximum transcript count as a proxy for induction rate.
Together, our analyses indicate that differential transcription kinetics in response to BMP signaling play a minor role in generating the distinct spatiotemporal expression patterns of BMP target genes.

Differential activation thresholds do not fully explain spatiotemporal expression
In the gradient threshold paradigm, target genes are activated by distinct signaling thresholds that define gene expression ranges ( Figure 1A). This model therefore predicts that broadly expressed genes, but not narrowly expressed genes, should be activated by low levels of signaling ( Figure 6A).
To test this idea, we exposed uninjected and Opto-BMP-injected embryos to high-(3900 lux) or low-intensity (70 lux) blue light for 10 or 20 min at shield stage -resulting in high-or low-amplitude BMP signaling pulses, respectively ( Figure 6B) -and then quantified BMP target gene responses using NanoString technology. As expected, target activation was generally stronger following higher amplitude, longer duration pulses ( Figure 6C-F and Figure 6-figure supplement 1). However, after a 10 min low-amplitude exposure, the third most narrowly expressed gene, foxi1, was significantly activated, whereas the broader genes were not robustly induced ( Figure 6C). A longer 20 min low-amplitude pulse significantly activated both narrowly and broadly expressed genes ( Figure 6E). A 10 min low-amplitude pulse significantly activated two of the top 50% earliest expressed genes (foxi1 and smad7), whereas a 20 min low-amplitude pulse significantly activated both early and lateexpressed genes ( Figure 6D,F). High-amplitude pulses activated genes of all ranges and activation times ( Figure 6C-F).  Figure 2) or spatial range (C,E). Error bars represent standard error, gray lines represent linear fits, r s = Spearman correlation coefficient, r p = Pearson correlation coefficient. crabp2b is not included due to lack of significant induction. (F-L) All three target gene response repeats were fitted with a model of induction and decay (Materials and methods). The average induction constant (s) is plotted against activation time (F), spatial range (G), or maximum transcript count (J). The average decay constant (l) is plotted against activation time (H), range (I), or maximum transcript count (K). Error bars represent standard error, r s = Spearman correlation coefficient, r p = Pearson correlation coefficient. crabp2b is not included due to lack of significant induction. pSmad1/5/9 immunofluorescence ( Figure gata2a  klf2b  bambia  tfap2c  eve1  foxi1  sizzled  ved  cdx4  gata2a  klf2b  bambia  tfap2c  eve1  foxi1  sizzled  ved  cdx4   gata2a  klf2b  bambia  tfap2c  eve1  foxi1  sizzled  ved  cdx4  gata2a  klf2b  bambia  tfap2c  eve1  foxi1  sizzled  ved  cdx4   klf2b  foxi1  gata2a  bambia  smad7  crabp2b  tfap2c  ved  sizzled  bmp4  eve1  cdx4  klf2b  foxi1  gata2a  bambia  smad7  crabp2b  tfap2c  ved  sizzled  bmp4  eve1  cdx4   klf2b  foxi1  gata2a  bambia  smad7  crabp2b  tfap2c  ved  sizzled  bmp4  eve1  cdx4  klf2b  foxi1  gata2a  bambia  smad7  crabp2b  tfap2c ved sizzled bmp4 eve1 cdx4 Figure 6. Differential sensitivity to BMP does not fully explain target gene expression diversity. (A) The activation threshold model predicts that broadly expressed genes will be activated by lower amplitude signaling. (B) pSmad1/5/9 immunofluorescence in uninjected and Opto-BMP-injected embryos exposed to 10 (triangle, dashed line) or 20 (diamond, solid line) min of 70 (light pink) or 3900 (magenta) lux blue light starting at shield stage.
Immunofluorescence was quantified and plotted as Opto-BMP signal -uninjected with piecewise linear interpolation between timepoints; error bars represent standard error (see Materials and methods for statistical analysis). Embryos for the transcriptional response experiment were collected 30 (x), 40 (square), or 50 (circle) min after the start of light exposure. (C-F) Transcriptional responses in Opto-BMP embryos exposed to conditions shown in (B) were quantified using NanoString technology and are plotted against spatial range (C,E) or activation Figure 6 continued on next page Our experiments exposing embryos to different amplitude BMP signaling pulses therefore suggest that not all spatiotemporal target gene expression differences are due to different signaling activation thresholds, although a subset may be (see Discussion).

FGF and Nodal modify BMP signaling and target gene expression
We noted that BMP target genes have unique expression patterns along the animal-vegetal axis that contribute to differences in their dorsal-ventral expression profiles ( Figure 1F-Y). Specifically, 6 out of the 10 spatially quantified high-confidence BMP target genes are either restricted to (cdx4, eve1) or excluded from (foxi1, klf2b, gata2a, tfap2c) the margin. We wondered how regulation by additional signaling pathways active at the margin might contribute to these differences. We focused on the FGF and Nodal pathways, which regulate mesoderm and mesendoderm specification, respectively, and are known to influence BMP signaling (Figure 7; Rogers and Müller, 2019).
To assess the effects of FGF and Nodal signaling on BMP target gene expression, we inhibited these pathways using the small molecule inhibitors SU-5402 (Mohammadi et al., 1997) (Varga et al., 2007;Londin et al., 2005;Koshida et al., 2002). Additionally, FGF restricts the expression of bmp (Londin et al., 2005;Fürthauer et al., 2004;Fürthauer et al., 1997), and we detected increased bmp2b expression in FGF/Nodal-inhibited embryos ( We reasoned that the shift of margin-excluded genes into the margin could either be due to loss of FGF/Nodal activity, or due to enhanced BMP signaling at the margin ( Figure 7D   . FGF and Nodal contribute to the spatial diversity of BMP target genes. (A) pSmad1/5/9 immunofluorescence in untreated embryos at the indicated times (data also shown in Figure 1E). (B-E) Embryos were treated with 10 mM FGF inhibitor SU-5402 (B), 50 mM Nodal inhibitor SB-505124 (C), or both (D) starting at 2 hpf, or injected with 0.5 pg bmp2b mRNA at the one-cell stage (E). The BMP signaling gradient was quantified along the dorsal-ventral axis at shield stage using pSmad1/5/9 immunofluorescence; error bars represent standard error. Note that the embryos in panels B-E came from different experiments and were processed and imaged on different days, but untreated controls were always siblings of treated embryos and processed and imaged simultaneously. Although transcript levels of several BMP targets were higher in treated compared to untreated embryos at later stages, their activation times were still diverse, suggesting that inputs other than FGF and Nodal are responsible for differences in activation times.
Finally, we noticed that much of the spatial diversity in BMP target gene expression along the dorsal-ventral axis collapsed in embryos lacking both FGF and Nodal signaling ( Figure 7F,I). To quantify the decrease in spatial diversity, we calculated the spatial coefficient of variation in untreated and treated embryos (see Materials and methods). Strikingly, embryos lacking both FGF and Nodal had lower coefficients of variation at almost all positions along the dorsal-ventral axis compared to untreated embryos ( Figure 7N). Together, our results identify combinatorial FGF and Nodal signaling as a major driver of spatial diversity in BMP target gene expression.

Discussion
Minor roles for differential responses to BMP in generating spatiotemporal diversity Signaling gradients are frequently observed in developing tissues, including the embryonic axes of gastrulating zebrafish, the neural tube in mice, and the wing precursor in Drosophila (Briscoe and Small, 2015;Schier and Talbot, 2005). However, how gradients are interpreted by cells is complex to ascertain. The gradient threshold model proposes that gene-specific activation thresholds are responsible for differences in the spatial expression of target genes (Sharpe, 2019;Briscoe and Small, 2015;Dubrulle et al., 2015;Rogers and Schier, 2011;Barkai and Shilo, 2009;Ashe and Briscoe, 2006). Can gradients be reliably generated and signaling thresholds accurately interpreted with high sensitivity, or do gradients simply provide a 'rough framework' for patterning that is refined over time by other mechanisms such as target gene cross-talk (Briscoe and Small, 2015;Chen et al., 2012) or cell sorting (Akieda et al., 2019;Xiong et al., 2013)? In the former case, is such precision actually required for patterning?
In the context of zebrafish dorsal-ventral patterning, our data suggest minor roles for gene-specific activation thresholds in generating BMP target gene expression diversity. We did not find a clear monotonically decreasing relationship between activation time and gene expression range       ( Figure 2O), suggesting that more broadly expressed genes are not consistently more likely to be activated by the low levels of BMP present early ( Figure 1D-E). We were also unable to detect an unambiguous correlation between range and the levels of signaling required for activation ( Figure 6C,E). This suggests that not all BMP target expression boundaries are positioned by genespecific BMP signaling thresholds ( Figure 1A).
An alternative model proposes that diversity in spatiotemporal target gene expression is due to differences in expression kinetics. For example, it was shown that Nodal targets with higher transcript accumulation rates in response to Nodal signaling have broader spatial expression domains (Dubrulle et al., 2015). To determine whether the BMP patterning system might function similarly, we examined the transcriptional responses of BMP target genes (Figures 4 and 5, and Figure 5figure supplement 1) to optogenetically generated pulses of BMP signaling ( Figure 3C-D and Figure 3-figure supplement 1I,J). We did not detect a strong correlation between transcript induction rates and activation time or spatial range ( Figure 5 and Figure 5-figure supplement 1G-W). Therefore, differential transcription kinetics in response to BMP are unlikely to account for spatiotemporal expression diversity.
Our results do not rule out the possibility that a different subset of BMP target genes may behave more consistently with these models. We focused on a set of high-confidence BMP targets (Figure 1-figure supplement 1), but other known targets were excluded from our analyses (Supplementary file 1). For example, the BMP target gene tp63 (Bakkers et al., 2002) is not expressed at shield stage, and was therefore excluded since it was not downregulated by chordin overexpression in our RNA-sequencing experiment (Supplementary file 1). We note that a subset of three genes (sizzled, ved, and bambia) that are neither restricted to nor excluded from the margin do show a monotonically decreasing relationship between range and activation time ( Figure 2O) as well as activation dynamics that could be roughly commensurate with signaling input ( Figure 6C and Figure 6-figure supplement 1), consistent with the gradient threshold model. However, it remains to be determined to what extent this subset of genes (or others) quantitatively follows the input-output relationships predicted by the gradient threshold model.
Our results also do not rule out other mechanisms of BMP signaling interpretation. For example, the graded distribution of many genes ( Figure 1P-Y) could be consistent with a model in which gene expression is roughly proportional to the level of BMP signaling. In addition, BMP signaling duration may encode specific responses in vivo. Future work is needed to better define the relationship between BMP signaling levels and gene expression and to determine how BMP signaling dynamics are interpreted in embryos. Our study highlights the promise of optogenetic approaches in such investigations (Rogers and Müller, 2020). In contrast to pharmacological or genetic methods, optogenetic strategies can provide fast, tunable, and reversible spatiotemporal manipulation of signaling in vivo ( Figure 3, Figure 6B, and Figure 3-figure supplement 1E-H), allowing more thorough characterization of input/output relationships.
In addition, our observations indicate that BMP signaling precision may not be required for proper patterning, or that the system is robustly buffered. For example, most embryos experiencing transient activation of BMP signaling lack gross morphological defects ( Figure Figure 5-figure supplement 1). How patterning recovers from such insults will be an interesting avenue for future study. Together with previous work (reviewed in Zinski et al., 2018), several of our observations indicate that feedback is an important feature of the BMP patterning system: Five out of 16 high-confidence BMP target genes affect BMP signaling (Figure 1-figure supplement 1), and embryos can experience a dip in signaling levels after a signaling pulse ( Figure 3C and Figure 3-figure supplement 1J). Cell sorting strategies that sharpen gene expression boundaries may also contribute to the observed recovery from BMP signaling manipulation (Akieda et al., 2019;Xiong et al., 2013).

Margin restriction and competence of BMP target genes
One unresolved question from our study is the restriction of the BMP target genes eve1 and cdx4 to the margin ( Figure 1I,I',L,L' and Figure 7-figure supplement 1). Consistent with previous work (Swanhart et al., 2010;Ota et al., 2009;Bennett et al., 2007;Ho et al., 2006;Londin et al., 2005;Shimizu et al., 2005;Rentzsch et al., 2004), in the absence of FGF or Nodal, eve1 and cdx4 were still expressed at the ventral margin (Figure 7-figure supplement 1; we note, however, conflicting reports with dominant-negative FGF receptors [Ota et al., 2009;Kudoh et al., 2004;Griffin et al., 1995]). Inhibition by animal pole factors or a requirement for signaling pathways at the margin such as Wnt or retinoic acid might play a role in their margin restriction.
Both eve1 and cdx4 are also activated relatively late in development ( Figure 2F,J), and cdx4 is not competent to respond to an early BMP signaling pulse ( Figure 4K). FGF and Nodal have no obvious roles in regulating their activation timing or competence since their temporal expression was not significantly affected by loss of FGF/Nodal signaling (Figure 7-figure supplement 2Q,U). Understanding how the activation timing of all BMP target genes including eve1 and cdx4 is regulated is an important future goal.  (Varga et al., 2007;Maegawa et al., 2006;Londin et al., 2005;Fürthauer et al., 2004;Kudoh et al., 2004;Koshida et al., 2002;Fürthauer et al., 1997), as well as inactivating Smad1 (Sapkota et al., 2007;Pera et al., 2003;Kretzschmar et al., 1997). Loss of Nodal did not detectably alter BMP signaling at shield stage. This is surprising because early expression of fgf is thought to depend on Nodal (van Boxtel et al., 2015;Maegawa et al., 2006;Mathieu et al., 2004;Gritsman et al., 1999;Rodaway et al., 1999), although low levels of fgf3 appear to be present at late blastula stages in Nodal signaling mutants (Mathieu et al., 2004), and weak FGF activity is detectable in Nodal inhibitor-treated embryos (van Boxtel et al., 2015). Nodal can activate chordin expression independently of FGF (Varga et al., 2007), and chordin is detectable albeit reduced in Nodal signaling mutants (Gritsman et al., 1999), suggesting that the reduction in chordin caused by Nodal loss is not sufficient to affect BMP signaling during early gastrulation. Future work is needed to explain why FGF, but not Nodal loss enhances BMP signaling at early gastrulation, and why simultaneous loss increases not only the amplitude but the broadness of the BMP signaling gradient.
Inhibition of FGF, Nodal, or both together also had distinct effects on BMP target gene expression ( Figure 7F-I, Figure 7-figure supplements 1, 2 and 3). Although Nodal loss did not detectably alter the BMP signaling gradient ( Figure 7C and The margin exclusion of the BMP target genes foxi1, klf2b, gata2a, and tfap2c can be explained by FGF/Nodal-mediated inhibition ( Figure 7K). Loss of either FGF or Nodal signaling shifted the expression of margin-excluded genes toward the margin, although the shifts were most dramatic in the absence of both ( Figure 7F Figure 7F,I,K,N), creating distinct dorsal-ventral profiles for margin-excluded genes ( Figure 1R,T,W,X) compared to non-excluded genes ( Figure 1P,Q,S,U,V,Y).
Our results suggest that much of the spatial diversity in BMP target gene expression arises from combinatorial signaling. A similar strategy is thought to regulate Bicoid target genes during Drosophila embryogenesis: Gene expression boundary shifts in response to Bicoid manipulation are often inconsistent with the gradient threshold model (Chen et al., 2012;Ochoa-Espinosa et al., 2009), and activation thresholds do not appear to explain target gene expression profiles at the DNA level (Ochoa-Espinosa et al., 2005). Rather, Bicoid is thought to act within a system of repressive pathways that regulate Bicoid target gene expression (Chen et al., 2012). During zebrafish dorsal-ventral patterning, FGF and Nodal affect BMP target gene expression in two ways: by restricting BMP signaling ( Figure 7B

RNA-sequencing
Wild type TE zebrafish embryos were dechorionated with Pronase (Roche, 11459643001) and injected at the one-cell stage with 10 pg mRNA encoding zebrafish Bmp2b, 100 pg mRNA encoding zebrafish Chordin, or left uninjected (Pomreinke et al., 2017). When uninjected siblings reached shield stage (~6.75 hpf), embryos were snap-frozen in liquid nitrogen. 10 embryos were collected per sample, three samples per condition.
To prepare total RNA, the TRIzol reagent (Invitrogen, 15596026) manufacturer's protocol was followed until aqueous phase recovery, then 6.25 ml Co-Precipitant Pink (Bioline, BIO-37075) was added to 250 ml aqueous phase, followed by 375 ml 100% EtOH. After vortexing briefly, samples were transferred to RNeasy Mini kit (Qiagen, 74104) spin columns and centrifuged at 13600 rpm at 4˚C for 1 min. Flow-through was discarded and columns were washed twice with RPE buffer (Qiagen). RNA was eluted in 50 ml H 2 O. Total RNA concentration was measured using a NanoDrop spectrophotometer (ThermoScientific). 3-5 mg total RNA per sample were provided to LCG Genomics GmbH (Berlin, Germany) for sequencing and differential expression analysis. Sequences were aligned against the reference genome Danio rerio GRCz10 with STAR 2.4.1b, and differential expression analysis was carried out with edgeR 3.2.3, DESeq 1.12.0, and Cuff diff 2.1.1. The p-value threshold for differentially expressed genes was set to 0.05.
Note that endogenous bmp2b and chordin were not distinguishable from injected mRNAs in bmp2b-or chordin-injected embryos, respectively, and were therefore excluded from consideration as BMP target genes.

Opto-BMP constructs
Opto-BMP constructs are based on Opto-Acvr constructs (Sako et al., 2016). These pCS2+-based Opto-Acvr constructs encode proteins that are tethered to the plasma membrane by an N-terminal myristoylation motif. Next to the membrane is a Nodal receptor kinase domain, followed by the light-oxygen-voltage (VfLOV) domain Aureochrome1 from Vaucheria frigida (Takahashi et al., 2007), and finally a C-terminal HA tag. Using splicing by overlap extension (SOE) PCR (Horton et al., 2013), Nodal receptor kinase domains in Opto-Acvr were swapped with putative kinase domains from the type I zebrafish BMP receptors Alk3 (NM_131621, bp 691-1566) (Nikaido et al., 1999) and Alk8 (NM_131345, bp 622-1497) (Mintzer et al., 2001;Yelick et al., 1998), and the type II zebrafish receptors BMPR2a (NM_001039817, bp 571-3009) and BMPR2b (NM_001039807, bp 598-1536) (Monteiro et al., 2008). In all cases except for Opto-BMPR2a, all residues after the transmembrane domain until the end of the kinase domain were included. Opto-BMPR2a contains all residues after the transmembrane domain until the end of the protein; the kinase domain-only construct was inactive.
An equimolar combination of mRNA encoding Opto-Alk3 (5.2 pg), Opto-Alk8 (5.2 pg), and Opto-BMPR2a (8.9 pg) was found to optimally induce BMP signaling in the light but not in the dark ( Figure 3B, Figure 3-figure supplement 1A,K,L), and was used in all Opto-BMP experiments described here.

LED array
To facilitate optogenetic experiments requiring control of light intensity and exposure duration, an embedded system-based controller was developed (Figure 3-figure supplement 1B-D). To maximize the versatility of the setup for different applications, a single-board computer was deployed (Raspberry Pi 3 model B, running under a Linux kernel, version 4.9). The controller was programmed to generate signals that modulate the duration and intensity of light. The generated signal was further amplified to drive the load of the LED array. A two-stage Darlington amplifier was used (TIP122 complementary power NPN Darlington -STMicroelectronics) to raise the ceiling of the current of amplification. The Darlington pair was used in a common emitter configuration in order to achieve a large power gain. The loads were operated on a constant voltage source provided by a regulated power supply (Disrelec Group AG, RND 320-KD3000D). During initial trials, brief, weak signal spikes could be detected, and an RC filter was subsequently used across the load to dampen any sporadic light flashes. The LED array constituted the circuit load; these LEDs were glued into the plastic cover of 6-well plates (Greiner Bio-One, 657160) (Figure 3-figure supplement 1B). Blue Nichia (NSPB510AS) or Everlight (1363-2SUBC/C470/S400-A4) LEDs were used in the array. Both LEDs emitted maximal spectral intensity at 470 nm, with the Nichia LEDs having a broader radiation angle, tighter spectral distribution, and less variable performance. During experiments, the LED array was placed inside a temperature-controlled incubator (Thermo Scientific Heratherum IMC 18, 50125882) set to 28˚C. Dark fabric was taped to the interior of the incubator door to prevent outside light from entering.
The circuit schematic (Figure 3-figure supplement 1D) shows how the generated square wave was used to drive the LED array. One of the Raspberry Pi's GPIO pins was used as a pulse-width modulation (PWM) output to produce signal. The raspberry-gpio-python module (https://sourceforge.net/projects/raspberry-gpio-python) was used to interface the GPIO. A pulse program was written in Python, which allows for variable parameter settings: GPIO pin number, modulation frequency (10 kHz is the NPN Darlington amplifier linear limit), pulse duration, and duty cycle.

Cycloheximide experiment
For the cycloheximide (Sigma, C4859) experiment in Figure 1-figure supplement 1D, embryos from wild type TE incrosses were dechorionated using Pronase (Roche, 11459643001) and injected at the one-cell stage with 5.2 pg opto-Alk3 + 5.2 pg opto-Alk8 + 8.9 pg opto-BMPR2a mRNA (Figure 3-figure supplement 1A). Control siblings were left uninjected, and embryos were sorted into agarose-coated 6-well plates and incubated at 28˚C. 70-90 min post-fertilization at the 4-16 cell stage, unfertilized and damaged embryos were removed, and plates were individually wrapped in aluminum foil to prevent light exposure and incubated at 28˚C. At 6.25 h post-fertilization (hpf), embryos were transferred into new agarose-coated 6-well dishes containing either 50 mg/ml cycloheximide (Bennett et al., 2007;Poulain and Lepage, 2002) or an equivalent volume of DMSO (Roth, A994.2) diluted in embryo medium that had been incubated at 28˚C prior to transfer. Red color filters (Rosco, E106 Primary Red) were used to cover the dissecting microscope light source during the transfer to prevent transmission of VfLOV-dimerizing wavelengths and minimize BMP activation, and plates were wrapped in aluminum foil after transfer. At 6.75 hpf (~shield stage, 30 min after cycloheximide exposure), plates were transferred to a small 28˚C incubator containing the LED array (Figure 3-figure supplement 1B) and exposed to blue light for 30 min (6.75-7.25 hpf). 20 min after light exposure, when most BMP target genes are maximally induced (Figure 4), embryos were fixed and colorimetric in situ hybridization was carried out as described in the Fluorescence and colorimetric in situ hybridization section below. pSmad1/5/9, pSmad2/3, and pErk immunofluorescence staining For pSmad1/5/9, pSmad2/3, and pErk immunofluorescence staining, embryos were fixed in 4% formaldehyde in PBS at 4˚C overnight, then transferred to MeOH and stored at À20˚C for at least 2 h. See below and Figure 1-figure supplement 1E for imaging and quantification details.

pSmad1/5/9
Embryos were washed at least three times with PBST (phosphate buffered saline + 0.1% Tween-20), then blocked for at least 1 h at room temperature in blocking buffer (10% FBS (Biochrom, S0415), 1% DMSO, 0.1% Tween-20 in PBS). Embryos were incubated in 1:100 rabbit anti-phosphoSmad1/5/ 9 antibody (Cell Signaling Technology, 13820) in blocking buffer at 4˚C overnight. One wash with blocking buffer followed by 3-5 washes with PBST were carried out at room temperature, then embryos were blocked again with blocking buffer for at least 1 h. Embryos were incubated in 1:5000 goat anti-rabbit Alexa Fluor 488-conjugated secondary antibody (Life Technologies, A11008) in blocking buffer at 4˚C overnight. Embryos were then incubated in 1:5000 DAPI (Life Technologies, D1306; stock concentration: 5 mg/ml) in blocking buffer at room temperature for at least 1 h, then washed at least five times with PBST. Stained embryos were wrapped in aluminum foil and stored at 4˚C overnight before SPIM imaging.

pSmad2/3
Embryos were incubated in ice-cold acetone (Roth, 5025.5) for 7 min, then washed at least three times with PBST, blocked for at least 1 h in 10% FBS in PBST and incubated in 1:5000 rabbit anti-pSmad2/3 (Cell Signaling Technology, 8828) in 10% FBS in PBST at 4˚C overnight. Embryos were then washed at least five times in PBST, blocked again for at least 1 h in 10% FBS in PBST, and incubated in 1:500 goat anti-rabbit HRP secondary antibody (Jackson ImmunoResearch, 111-035-003) in 10% FBS in PBST at 4˚C overnight. Next, embryos were washed at least five times in PBST, then once in TSA 1x amplification buffer (TSA Plus Cyanine 3 System, Perkin Elmer, NEL744001KT). For staining, embryos were incubated in 75 ml 1:75 Cy3-TSA in 1x amplification buffer in the dark at room temperature for 45 min. After washing at least six times with PBST, embryos were incubated in 1:5000 DAPI (Life Technologies, D1306; stock concentration: 5 mg/ml) in PBST at room temperature for at least 1 h, then washed at least four times with PBST. Finally, embryos were wrapped in aluminum foil and stored at 4˚C overnight before SPIM imaging.

pErk
Embryos were washed at least three times with PBST, then transferred to ice-cold acetone for 20 min and washed at least three times with PBST. After blocking in 10% FBS in PBST for at least 1 h, embryos were incubated in 1:5000 mouse anti-pErk antibody (Sigma, M8159) in 10% FBS in PBST at 4˚C overnight. Embryos were then washed at least five times in PBST, blocked again for at least 1 h in 10% FBS in PBST, and incubated in 1:500 donkey anti-mouse HRP secondary antibody (Jackson ImmunoResearch, 715-035-150) in 10% FBS in PBST at 4˚C overnight. Embryos were washed at least five times with PBST, then once in TSA 1x amplification buffer. Next, embryos were incubated in 75 ml 1:75 Cy3-TSA in 1x amplification buffer in the dark at room temperature for 45 min. After washing at least six times with PBST, embryos were incubated in 1:5000 DAPI (Life Technologies, D1306; stock concentration: 5 mg/ml) in PBST at room temperature for at least 1 h, then washed at least four times with PBST. Stained embryos were wrapped in aluminum foil and stored at 4˚C overnight before SPIM imaging.
Fluorescence and colorimetric in situ hybridization BMP target gene probes were generated by amplifying full or partial coding sequences (CDS) from wild type TE zebrafish cDNA and cloning into pCS2+ or pCR-bluntII TOPO (ThermoFisher, 450245) vectors. Plasmids were linearized with the indicated restriction enzymes, column purified (Promega, A9282), and DIG-labeled probes were generated using the indicated polymerase (Roche, 11175025910).
High-confidence BMP target gene in situ hybridization probes: Note that the znfl2b in situ probe contained 47 SNPs compared to the reference genome (Danio rerio GRCz11).
The Whole-mount colorimetric in situ hybridization was carried out as described previously (Thisse and Thisse, 2008). Embryos were fixed in 4% formaldehyde in PBS, incubated at 4˚C overnight, then transferred to MeOH and stored at À20˚C for at least 2 h. Stained embryos were imaged in 2:1 benzyl benzoate:benzyl alcohol with an Axio Zoom.V16 microscope (ZEISS).
For fluorescence in situ hybridization (FISH), the same protocol was used until the blocking step, at which point embryos were blocked in FISH blocking buffer (2% blocking reagent (Roche, 11096176001) in 1x maleic acid buffer (100 mM maleic acid, 150 mM NaCl, 180 mM NaOH, 0.1% Tween)) for at least 2 h at room temperature with gentle rocking, then incubated at 4˚C overnight in 1:150 anti-DIG-POD (Roche, 11207733910). The following day embryos were washed at least five times with PBST. To develop signal, embryos were incubated in 75 ml 1:75 Cy3-TSA in 1x amplification buffer (TSA Plus Cyanine 3 System, Perkin Elmer, NEL744001KT) for 30 min at room temperature in the dark. Embryos were then washed at least five times with PBST, incubated in 1:5000 DAPI (Life Technologies, D1306; stock concentration: 5 mg/ml) with agitation at room temperature for at least 1 h (or overnight at 4˚C), then washed at least five times with PBST. One day after Cy3 incubation, embryos were imaged on a ZEISS Lightsheet Z.1 (see below and Figure 1-figure supplement 1E for imaging and quantification details). All FISH embryos shown in Figure 1F-Y were fertilized and fixed on the same day.

SPIM imaging of immunofluorescence staining and fluorescence in situ hybridization
Fixed embryos were mounted in 1% low melting temperature agarose (Lonza, 50080) using a glass capillary and imaged with a ZEISS Lightsheet Z.1 selective plane illumination microscope (SPIM). The imaging chamber was filled with water, and filters and light sheets were auto-aligned prior to imaging. For fluorescence in situ hybridization (FISH) and pSmad1/5/9 immunofluorescence (IF) experiments, embryos were positioned using the DAPI signal with the animal pole pointing toward the imaging objective to produce animal views; for ventral views, embryos in the correct orientation were rotated 90˚. For animal views, 50-90 z-slices with 7 mm between each slice were acquired per embryo, covering the entire blastoderm over a distance of 350-630 mm depending on embryo size. For ventral views,~70 z-slices with 7 mm between each slice were acquired per embryo, spanning roughly half of the embryo.
For pSmad2/3 and pErk IF, embryos were mounted in the orthogonal orientation compared to pSmad1/5/9 and FISH experiments, and three lateral images were acquired per embryo: one at the brightest region, a second rotated 120˚, and a third rotated 240˚.
All images were acquired with dual light sheet illumination using a W Plan-Apochromat 20x objective at 0.5x zoom and the imaging conditions described below.
SPIM imaging conditions: Maximum intensity projections were generated using the software ZEN (2014 SP1, black edition) and used for the analyses described below.

Mathematical modeling of target gene induction and decay kinetics
To estimate induction and decay of transcripts from the NanoString data (Figure 4), time-dependent pSmad1/5/9 and transcript changes were modeled mathematically. The change in the amount of endogenous (P e ) and optogenetically induced (P o ) pSmad1/5/9 levels can be described by the following general differential equations: The observed pSmad1/5/9 levels in uninjected embryos correspond to G(t), whereas the observed pSmad1/5/9 levels (P s ) in light-exposed Opto-BMP embryos correspond to the sum of G(t) and H(t). Therefore, the change in the amount of P s over time can be described by: Thus, the levels of optogenetically induced pSmad1/5/9 can be calculated by subtracting the pSmad1/5/9 levels in uninjected embryos from the pSmad1/5/9 levels in light-exposed Opto-BMP embryos: Similarly, changes in the endogenous transcript levels (T e ) and optogenetically induced transcript levels (T o ) over time can be described by the following general differential equations: The observed transcript levels in uninjected embryos correspond to K(t), whereas the observed transcript levels (T s ) in light-exposed Opto-BMP embryos correspond to the sum of K(t) and L(t). The change in the amount of T s over time can therefore be described by: Thus, the levels of optogenetically induced transcripts can be calculated by subtracting the transcript levels in uninjected embryos from the transcript levels in light-exposed Opto-BMP embryos:

Modeling method 1
The NanoString transcription data was first analyzed using the simplest model of induction and decay ( Figure 5): where P o represents the optogenetically induced pSmad1/5/9 input, T o the pSmad1/5/9-dependent target gene, s the induction rate constant, and l the decay rate constant of the induced gene. P o was obtained by fitting the measured pSmad1/5/9 immunofluorescence data H(t) ( Figure 3C, Figure 5L, and Figure 3-figure supplement 1J) with a polynomial of degree five using the function polyfit in MATLAB 7.10.0 (R2010a). The induction-decay model was simulated in COMSOL Multiphysics 3.5a in a 10 mm domain (representing approximately one cell) with no-flux boundary conditions and an initial concentration T o (0). For each experiment, the combination of parameters T o (0), s, and l was found that minimizes the sum of squared differences (SSD) between the simulations of the induction-decay model T o (t n ) and the data L(t n ) for all measured time points n. The minimization was performed numerically using a constrained optimization algorithm (Nelder-Mead, MATLAB 7.10.0) with zero for the initial guesses of T o (0), s, and l, and a maximum of 500 iterations. s and l were constrained between biologically plausible values of 0.00001/s and 0.1/s, and T o (0) was bounded between À100 a.u. and 100 a.u. R 2 values were calculated from the minimizing SSD (SSD min ) to assess the goodness of the fits by Fitted values for high-confidence BMP target genes, experimental repeat 1:

Modeling method 2
In a second approach (Figure 5-figure supplement 1G-W), the NanoString transcription data was fitted with the analytical solutions to the differential equation system which describes the changes in endogenous as well as optogenetically induced pSmad1/5/9 and transcript levels based on the simplest model of induction and decay after an optogenetic pulse of length t L (i.e., 30 min = 1800 s for all experiments). k 1 represents the activation rate of endogenous pSmad1/5/9, k 2 the decay rate constant of pSmad1/5/9, and k 3 the activation rate of optogenetically induced pSmad1/5/9. Optogenetic switch-like activation was modeled with the Heaviside step function . k 4 and k 5 represent the activation rate and decay rate constants of endogenously induced BMP-dependent transcripts, and s and l are the induction rate and decay rate constants of the induced gene.
The analytical solutions to this equation system are: The pSmad1/5/9 data was fitted with the computer algebra system Maple (Waterloo Maple Inc) using the function LSSolve to minimize the difference between the pSmad1/5/9 data in uninjected embryos and P e (t), as well as the difference between the pSmad1/5/9 data in light-exposed Opto-BMP embryos and P e (t) + P o (t) with the initial guesses d Pe ¼ 0 a:u:, d Po ¼ 0 a:u:, k 1 ¼ 0=s, k 2 ¼ 0:00167=s, k 3 ¼ 0=s, k 4 ¼ 00167=s and a maximum of 20000 iterations and an optimality tolerance of 0.3981071706 Â 10 À14 . The best fitting parameters d Pe ¼ À76:19 a:u, d Po ¼ 264:1 a:u:; k 1 ¼ 0:1429 a:u:=s; k 2 ¼ 0:000900=s; and k 3 ¼ 0:954 a:u:=s were then used for the simulation of the gene induction dynamics in the NanoString data.
The NanoString data was fitted in Maple using the function LSSolve to simultaneously minimize the difference between the NanoString data in uninjected embryos and T e (t), as well as the difference between the NanoString data in light-exposed Opto-BMP embryos and T e (t) + T o (t) with the initial guesses d Te ¼ 0 a:u:, d To ¼ 0 a:u:, k 4 ¼ 0=s, k 5 ¼ 0:00167=s, s ¼ 0=s, l ¼ 0:00167=s and a maximum of 10000 iterations and an optimality tolerance of 0.3981071706 Â 10 À14 .
Fitted values for high-confidence BMP target genes:
Quantification of pSmad1/5/9 immunofluorescence staining and fluorescence in situ hybridization  (Schindelin et al., 2012) with ventral to the left (brightest signal) and dorsal to the right (dimmest signal; for the very early pSmad1/5/ 9 images prior to clear onset of BMP signaling, embryos were oriented with the brightest side on the left and the dimmest on the right where obvious, but correspondence with ventral-dorsal is not clear in those early cases). A polygonal region of interest (ROI) was then manually drawn around the embryo and used to create a mask in order to remove image background (for FISH experiments, the Cy3 signal was used to draw the mask; for IF experiments the DAPI signal was used). The average pixel intensity in each column of pixels from ventral to dorsal was then acquired (pixel area: 0.46 mm Â 0.46 mm). For genes that are restricted to the margin (cdx4 and eve1), a second manually positioned circular ROI was used to exclude the non-margin region of the embryo (Figure 1-figure  supplement 1E).
For FISH experiments, non-probe-exposed control embryos for background subtraction were imaged and intensity profiles acquired as described above. The orientation of these background subtraction embryos was random. Images for background subtraction controls were acquired in the same imaging session as experimental FISH images.
After intensity profiles were acquired, absolute distance was converted into percent embryo length to account for embryo-to-embryo variability in size, and intensity measurements were averaged into bins of 0.5% embryo length using an automated routine (0 < bin 1 < 0.5%, 0.5 < bin 2 < 1%, etc.).
For FISH experiments, the average intensity at each position in all 10 non-probe-exposed background embryos was calculated. This spatial background average was subtracted from each experimental FISH raw intensity profile, and data from the first and last 5% embryo length was excluded because the averages at the most ventral and dorsal regions are composed of relatively few pixels and are therefore less reliable.
The profiles of individual embryos were normalized following the procedure in Gregor et al., 2007 using the model which relates the mean intensity profile c À x ð Þ of all data points for a given target gene to the intensity profile I n x ð Þ for an embryo n through the embryo-specific proportionality constant A n and the nonspecific background b n . A n and b n were determined by minimizing the objective function i X ðI n ðx i Þ À ðA n c À ðx i Þ þ b n ÞÞ 2 for the data points at all positions x i with the Nelder-Mead algorithm using the function fminsearch in MATLAB 7.10.0, the initial guesses 1 and 0 for A n and b n , a maximum of 10000 function evaluations, and a maximum of 5000 iterations. For display, each average profile was then divided by its maximum intensity ( Figure 1P-Y, Figure 7F-J, and Figure 7-figure supplement 3).
The Gaussian function Ae À xÀ ð Þ 2 & was fitted to each profile using a constrained Nelder-Mead algorithm in MATLAB 7.10.0 with a maximum of 10000 function evaluations, a maximum of 5000 iterations, the initial guesses 300, 20, and 10000, the lower bounds 300, -50, and 100, and the upper bounds 100000, 50, and 100000 for A, m, and &, respectively. Gene expression range was defined as . The resulting ranges from 9-10 embryos were averaged to define each gene's mean range.
For pSmad1/5/9 IF spatial quantification experiments, the average image background intensity was determined for each image using a small ROI in the corner outside of the embryo, and subtracted from each IF raw intensity profile. Since the averages at the most ventral and dorsal regions are composed of relatively few pixels and are therefore less reliable, data from the first and last 5% embryo length was not considered. The mean of the dorsal-most 5% at 2.75 hpf was then subtracted from all profiles. These profiles were then normalized as described above for the FISH data, assuming embryo-specific constant nonspecific background and proportionality constants that relate immunofluorescent staining intensity to protein concentration.
Number of embryos assessed in spatial quantification experiments: To quantify total pSmad1/5/9 IF intensity (Figures 3C-E and 6B, and Figure 3-figure supplement 1I,J), an ROI was manually drawn around the embryo in Fiji based on DAPI signal and used to create a mask in order to remove image background as described above. The average intensity within the ROI was then calculated.
For experiments shown in Figure 3C-D and Figure 3-figure supplement 1I, image background intensity was measured using a small ROI in the corner of each image outside of the embryo. The average image background was then subtracted from the embryo intensity measurements to generate background-subtracted intensities.
For shield-stage experiments shown in Figure 3C,E, Figure 5L, Figure 6B, and Figure 3-figure supplement 1J, the average intensity within a small ROI on the dorsal side was measured in uninjected embryos; for each time point, these values were averaged and subtracted from the embryo intensity measurements to generate background-subtracted intensities.
Number of embryos assessed in total pSmad1/5/9 IF quantification time course experiments:

NanoString RNA quantification
For the NanoString time course experiment in untreated ( Figure 2) and FGF/Nodal-inhibitor-treated embryos (Figure 7-figure supplement 2L-ZA), embryos from wild type TE incrosses were collected~15 min after mating commenced. Embryos were incubated at 28˚C, dechorionated using Pronase (Roche, 11459643001) at~1.5 hpf, and sorted into 10 agarose-coated 6-well plates, one plate per time point. Each plate had one well containing embryo medium and one well containing FGF/Nodal inhibitor. To keep temperature and therefore development steady, plates were only removed from the 28˚C incubator immediately prior to embryo collection. Every 30 min from 2.75 to 7.25 hpf, treated and untreated embryos were snap-frozen in liquid nitrogen.
For NanoString experiments quantifying responses to BMP signaling pulses using Opto-BMP (Figures 4, 5 and 6, Figure 5-figure supplement 1J-W, and Figure 6-figure supplement 1), embryos from TE incrosses were dechorionated using Pronase and injected at the one-cell stage with 5.2 pg opto-Alk3 + 5.2 pg opto-Alk8 + 8.9 pg opto-BMPR2a mRNA (Figure 3-figure supplement 1A). Control siblings were left uninjected, and embryos were sorted into agarose-coated 6well plates and incubated at 28˚C. 70-90 min post-fertilization at the 4-16 cell stage, unfertilized and damaged embryos were removed, and plates were individually wrapped in aluminum foil and incubated at 28˚C. At the appropriate time, individual plates were transferred to a small 28˚C incubator containing the LED array, exposed to light for the appropriate duration, and embryos were either snap-frozen in liquid nitrogen immediately (e.g., for the 10 min during exposure time point), or rewrapped in aluminum foil and returned to 28˚C incubation in the dark (e.g., for the 80 min postexposure time point). RNA was prepared as described for the RNA-sequencing experiment. 30 ml aliquots at 20 ng/ml were provided to Proteros GmbH (Planegg-Martinsried, Germany) for analysis using a customdesigned NanoString nCounter Elements TagSet with probes targeting high-confidence BMP target genes identified by RNA-sequencing, and housekeeping genes for normalization. Samples were measured using an nCounter SPRINT according to the standard protocol with a 24-30 h hybridization length. nSolver 4.0 software (https://www.nanostring.com/products/analysis-software/nsolver) was used to subtract background and normalize the RNA count data using the geometric means of the positive spike-in controls and the housekeeping genes eef1a1l1 and act2b, respectively. Lanes that failed quality control were repeated.
Number For experiments in which transcriptional responses to BMP signaling pulses are assessed (Figures 4, 5 and 6C-F, and Figure 6-figure supplement 1), it is necessary to determine changes in transcript levels compared to uninduced embryos. Because each of the three sets of Opto-BMP embryos had uninjected control siblings collected at the same time, average induction was calculated by first subtracting the uninjected transcript count from its corresponding injected sibling count, then by averaging the three subtracted counts (also see the section Mathematical modeling of target gene induction and decay kinetics above for a formal description of this procedure).

Calculation of spatial coefficients of variation
The spatial coefficient of variation ( Figure 7N) for each condition (untreated, bmp-overexpressing, +SB-505124, + SU-5402, and +SB-505124 and SU-5402) was calculated as follows: First, at each position x, the average normalized intensity and standard deviation s for all n genes quantified by FISH were determined ( Figure 7F-J). Next, the standard deviation was divided by the average normalized intensity at that position This was repeated for every position along the dorsal-ventral axis for all five conditions to calculate the spatial coefficients of variation for the 10 measured genes.

Statistical analyses
In the following experiments, significance was defined as a p-value 0.05 using an unpaired twotailed Student's t-test assuming equal variance in Excel.
To determine how light exposure at different developmental stages affects BMP signaling in Opto-BMP embryos, total pSmad1/5/9 immunofluorescence intensity was quantified in uninjected and Opto-BMP-injected embryos exposed to light at high (3.5-4 hpf) or shield (6.75-7.25 hpf) stage ( Figure   To determine how different light intensities affect BMP signaling in Opto-BMP embryos, total pSmad1/5/9 immunofluorescence intensity was quantified in uninjected and Opto-BMP-injected embryos exposed to low (70 lux) or high (3900 lux) intensity light for 10 or 20 min ( Figure 6B).
Low-and high-intensity light, Opto-BMP versus uninjected p-values ( Figure 6B): To determine whether BMP target gene expression domain boundaries differ significantly in untreated embryos, range was defined in individual embryos as described in the section Quantification of pSmad1/5/9 immunofluorescence staining and fluorescence in situ hybridization. Ranges were then averaged.
BMP target gene range comparison p-values ( Figure 1P- The shape of the temporal BMP target gene expression profiles assessed by NanoString in untreated and SU-5402/SB-505124-treated embryos can be well approximated by the modified cumulative distribution function of the normal distribution which was used for regression analysis using a constrained Nelder-Mead algorithm in MATLAB 7.10.0 with a maximum of 10000 function evaluations, a maximum of 5000 iterations, the initial guesses 1000, 5 h, 1 h, and 100, the lower bounds 100, 3 h, 0.05 h, and 0, and the upper bounds 10000, 7 h, 3 h, and 1000 for A, n, t , and b, respectively. The activation time of each BMP target gene was defined as the average time point at which the curves reached about two mean average deviations (i.e., 1:5 Á t ) from the inflection point n (Figure 2 and Figure 7-figure supplement 2L-Y).
id2a (Chong et al., 2005) and smad6a (White et al., 2017) were excluded from this analysis because they are maternally contributed.
To determine whether FGF/Nodal loss affects the timing of gene activation, activation times in untreated versus SU-5402/SB-505124-treated samples were compared (  To identify differences in BMP target gene expression in the absence of FGF/Nodal signaling, transcript counts from SU-5402/SB-505124-treated embryos were compared to counts from untreated embryos (Figure 7-figure supplement 2L-Y).
SU-5402/SB-505124-treated versus untreated p-values (Figure 7- For experiments in which transcriptional responses to BMP signaling pulses at high or shield stage were measured using NanoString (Figures 4 and 5M-Z), mRNA counts in Opto-BMP-injected embryos were compared to uninjected embryos.
High-stage BMP signaling pulse, Opto-BMP versus uninjected p-values ( Figure 4): For experiments in which transcriptional responses to low-and high-amplitude BMP signaling pulses of different durations were measured using NanoString ( Figure 6C-F and Figure 6-figure supplement 1), mRNA counts from uninjected embryos were first subtracted from Opto-BMPinjected siblings. Then the subtracted counts from light-exposed embryos were compared to subtracted counts from unexposed control embryos.
Low-and high-amplitude BMP pulses, exposed versus unexposed p-values ( Figure 6C

Data and code availability
The raw images, data, and source code for custom scripts used in this work are available from the corresponding author upon request. Image quantification data and differential gene expression analyses are available in the accompanying source data files and Supplementary file 1, respectively. The RNA-sequencing data has been deposited at the GEO repository (accession number: GSE135100).
. Transparent reporting form

Data availability
The RNA-sequencing data has been deposited at the GEO repository (accession number: GSE135100) and can be accessed at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE135100. Image quantification data is available in the accompanying source data files.
The following dataset was generated: Author (