The BMP signaling gradient is interpreted through concentration thresholds in dorsal-ventral axial patterning.

Bone Morphogenetic Protein (BMP) patterns the dorsal-ventral (DV) embryonic axis in all vertebrates, but it is unknown how cells along the DV axis interpret and translate the gradient of BMP signaling into differential gene activation that will give rise to distinct cell fates. To determine the mechanism of BMP morphogen interpretation in the zebrafish gastrula, we identified 57 genes that are directly activated by BMP signaling. By using Seurat analysis of single-cell RNA sequencing (scRNA-seq) data, we found that these genes are expressed in at least 3 distinct DV domains of the embryo. We distinguished between 3 models of BMP signal interpretation in which cells activate distinct gene expression through interpretation of thresholds of (1) the BMP signaling gradient slope; (2) the BMP signal duration; or (3) the level of BMP signal activation. We tested these 3 models using quantitative measurements of phosphorylated Smad5 (pSmad5) and by examining the spatial relationship between BMP signaling and activation of different target genes at single-cell resolution across the embryo. We found that BMP signaling gradient slope or BMP exposure duration did not account for the differential target gene expression domains. Instead, we show that cells respond to 3 distinct levels of BMP signaling activity to activate and position target gene expression. Together, we demonstrate that distinct pSmad5 threshold levels activate spatially distinct target genes to pattern the DV axis.


Introduction
During embryonic patterning, the molecular identity of unspecified cells is determined by the location of each cell within the embryo [1]. Thereby, a cell's position becomes translated into a specific cell fate. This positional information is provided by gradients of secreted signaling molecules called morphogens, which induce the specification of multiple cell fates [2][3][4][5]. The genetic programs underlying different cell fates are activated in distinct regions of tissue within the morphogen gradient. The proper position of gene expression boundaries is essential during development because these domains determine the differentiation and relative abundance mechanisms positioning target gene expression. We determined that cells respond to at least 3 distinct levels of pSmad5 to activate different target genes, and these threshold levels of pSmad5 can precisely position gene expression boundaries in the embryo.

Identification of target genes directly patterned by the BMP gradient
BMP signaling is essential for ventral tissue specification, but it remains unknown how graded BMP signaling is interpreted by cells of the embryo to generate the distinct gene expression domains that pattern the DV axis. Only a limited number of genes directly regulated by BMP signaling have been identified in the zebrafish gastrula. To identify the genes responding to the BMP gradient, we determined which genes are directly activated by pSmad5 from all BMPdependent gene expression during gastrulation. BMP-dependent gene expression was determined by performing RNA sequencing (RNA-seq) on wild-type and bmp7 mutant embryos (bmp7a sb1aub ) at 2 time points when the BMP signaling gradient patterns ventral tissues: early (shield) and mid-gastrula (70% epiboly) stages [22,36,37]. All ventral tissue specification is absent in bmp7 mutants, and the activation of Smad5 is abolished (Fig 1A, S1A and S1B Fig) [20,21,38]. We identified 1,559 genes that were differentially expressed (false discovery rate [FDR] < 0.05) in bmp7 mutant compared to wild-type embryos at an early gastrula stage (S1C Fig) and 852 genes that were differentially expressed at mid-gastrulation ( Fig 1B). Most differentially expressed genes were down-regulated in bmp7 mutants compared to wild-type embryos at both early and mid-gastrula stages, including many known markers of ventral tissues reflecting a loss of ventral cell fates ( Fig 1B, S1C Fig).
To identify the subset of BMP-dependent genes that are directly regulated by BMP signaling, we treated bmp7 mutant embryos at 4 hours post fertilization (hpf) with cycloheximide (CHX), a translation inhibitor, and then injected BMP2/7 recombinant protein into the intercellular space of the blastula (Fig 1C). First, we chose 4 hpf because the zygotic genome has been activated by this time point [39][40][41], but the DV axis has yet to be patterned, and all cells of the embryo remain competent to respond to BMP signaling [22]. Second, translation was inhibited with CHX to prevent the expression of secondary targets but not direct targets of BMP signaling. Finally, embryos were injected with BMP2/7 protein to activate BMP signaling and induce robust phosphorylation of Smad5 in bmp7 mutant embryos ( Fig 1D).
Total RNA was isolated 1.5 hours postinjection from bmp7 mutant embryos treated with CHX with or without BMP2/7 protein injection for RNA-seq analysis. We identified 363 genes that were differentially expressed (FDR < 0.05) in embryos injected with BMP2/7 protein ( Fig  1E). In BMP2/7 injected and CHX-treated embryos, a known direct target of BMP signaling, foxi1 [42], was confirmed to be expressed by in situ hybridization 1.5 hours after injection (S1D Fig). We compared the 274 genes that were up-regulated by BMP signaling after CHX treatment to the genes that were BMP-dependent during gastrulation. We found 57 genes that are both directly up-regulated by BMP signaling and endogenously expressed during gastrulation when ventral tissues are specified (S1 Table). Gene Ontology (GO) analysis of the 57 target genes found enrichment of terms for cell fate specification and tissue differentiation (Fig 1F). The BMP target genes were also enriched for GO terms for DNA-binding transcription factor activity (S1E Fig), which together are consistent with roles in specifying ventral cell fates. Fifteen genes were found to be both down-regulated by BMP signaling after CHX treatment and down-regulated in wild-type embryos compared to bmp7 mutants indicating that pSmad5 can directly inhibit gene expression (S2 Table). There is evidence that BMP-activated Smad transcription factors directly repress transcription via recruitment of repressors and chromatin modifiers [43][44][45]. Ten of these down-regulated target genes are known to be expressed within dorsal ectodermal tissue, which is consistent with a role in direct down-regulation by BMP signaling [23,[46][47][48][49][50][51][52][53][54]. Thus, we identified genes that directly read out the BMP signaling gradient to initiate the genetic cascade that specifies different ventral cell fates [37,42,[55][56][57]. We now have a comprehensive list of genes directly responding to the BMP gradient during DV patterning.

Ventral BMP target genes are expressed in at least 3 distinct domains
Next, we determined where the genes reading out the BMP signaling gradient are expressed along the DV axis of the embryo. Target genes responding to different aspects of the gradient would be predicted to be expressed in different domains of the embryo. To sort the target genes based on their expression pattern, we analyzed a previously published single-cell RNA sequencing (scRNA-seq) dataset of mid-gastrula zebrafish embryos using the Seurat method, which predicts the spatial position of individual cell transcriptomes [58,59]. Locations of single-cell transcriptomes were inferred by Seurat based on the co-expression of known landmark genes and mapped onto an embryonic grid that was divided into 64 bins: 8 bins across the DV axis and 8 bins across the animal-vegetal (AV) axis ( Fig 2D). The predicted expression patterns of the BMP direct target genes within the 64 bins across the embryo are shown as a heat map (Fig 2A-2D, S2-S5 Figs). To visualize the predicted expression profiles across the DV axis, we summed expression in all bins within each of the 8 positions along the DV axis ( Genes were then clustered based on their predicted expression profiles across the DV axis. Fifty of the 57 up-regulated target genes were found in the scRNA-seq dataset, and 27 genes showed predicted ventrally restricted expression (S2-S4 Figs). These ventrally restricted target genes were further sorted based on the bins where they were predicted to be expressed across the DV axis. Five genes had predicted expression enriched within the first 3 to 4 bins, 6 genes were predicted to be expressed within the first 5 bins, and 19 genes within the first 6 to 7 bins (Fig 2D, S2-S4 Figs). Target genes with uniform expression or that also included dorsal expression were assumed to have multiple signaling inputs and excluded from our analysis (S5 Fig, S1 Table). The down-regulated target genes that were sequenced in the scRNA-seq dataset displayed either dorsal enrichment or were not enriched along the DV axis (S6 Fig, S2 Table). To further investigate how the pSmad5 gradient is interpreted into spatially distinct gene expression domains, we focused on genes that are induced exclusively by BMP signaling. Some of the BMP target genes are known to have multiple signaling pathways contributing to the total expression pattern. For example, the expression of eve1 is known to be activated by both BMP as well as by Fibroblast Growth Factor (FGF) signaling [60][61][62]. Similarly, we also avoided examining genes that have significant expression remaining in bmp7 mutant embryos (S1 Table). Animal view, dorsal to right of mean pSmad5 intensities at mid-gastrula stage (8 hpf) WT (n = 4) and bmp7 mutant (n = 3) embryos. (B) Differential gene expression in WT and bmp7 mutants at 8 hpf using RNA-seq. Significantly up-regulated genes in WT compared to bmp7 mutants are shown in red, and significantly down-regulated genes are shown in blue. All other genes are shown in gray. A subset of known BMP-dependent genes is highlighted. See Table A in S1 Data for underlying data. (C) Schematic of method to isolate RNA for sequencing from CHX-treated bmp7 mutant embryos injected with BMP2/7 protein. (D) Representative image of pSmad5 intensities of all nuclei of an individual bmp7 mutant (n = 9) and bmp7 mutant injected with 10 pg of BMP2/7 protein (n = 9). Animal pole facing up. (E) Differential gene expression using RNA-seq of CHX-treated bmp7 mutants versus CHX-treated bmp7 mutants injected with BMP2/7 protein. Significantly up-regulated genes with BMP2/7 protein injection are shown in red, and significantly down-regulated genes are shown in blue. All other genes are shown in gray. A subset of BMP direct target genes is highlighted. See Table B in S1 Data for underlying data. (F) GO term analysis for biological processes of the 57 direct target genes. See Table C   For representative BMP targets, we chose 1 gene from each cluster to investigate how the gradient of BMP signaling directs target gene expression: sizzled in cluster 1, foxi1 in cluster 2, and bambia in cluster 3. To validate the predicted expression patterns of the 3 target genes, we performed fluorescent in situ hybridization (FISH) for each gene on early gastrula (7 hpf) wild-type embryos (S7A Fig), when the BMP gradient is patterning ventral fates [22,37]. Each embryo was subdivided into 128 bins equally spaced across the AV and DV axes (S7B Fig). Expression intensity for each lateral half of the embryo was averaged into 64 bins, normalized, and displayed similarly as the Seurat heat map (Fig 2A'-2C'). The measured and predicted expression profiles of the 3 target genes are similar across the DV axis (Fig 2A"-2C") validating this approach for examining the DV expression domain. While Seurat accurately predicted differences in expression patterns of the 3 genes across the DV axis, sizzled and bambia were incorrectly predicted to be uniformly expressed along the AV axis. Using the quantified FISH for these target genes and other landmark genes will improve the accuracy of Seurat's predicted expression pattern across the AV axis of the embryo.
To more precisely determine where the 3 target genes are expressed across the DV axis, we quantified the FISH expression intensity (Fig 2E-2G

Distinct pSmad5 levels and gradient slopes correspond to different target gene expression boundaries
Next, we determined where the boundaries of the 3 BMP target genes are located along the pSmad5 gradient. We took sibling embryos of those stained for FISH and quantified pSmad5 immunostaining at an early gastrula (7 hpf) stage ( Fig 3A). To identify the pSmad5 position that corresponds to each target gene expression boundary, we overlaid the DV position of the target gene boundaries onto the pSmad5 gradient profile of the sibling embryos. The pSmad5 gradient position at each expression boundary could indicate a distinct pSmad5 gradient slope or threshold that cells need to reach to induce differential expression of each gene. If gene expression is determined by a pSmad5 threshold level, then cells will not express a target gene until they reach the particular threshold level of pSmad5. Alternatively, if gene expression boundaries are positioned by distinct gradient slopes, then cells will activate target genes in response to a particular steep or shallow slope of the gradient independent of specific pSmad5

PLOS BIOLOGY
levels. It is also possible that 1 gene may respond to slope, while the others respond to distinct thresholds, for example. The same mechanism need not apply to all 3 target genes.
We found that expression of the 3 target genes correlates with significantly distinct pSmad5 gradient levels ( Fig 3B). The expression boundary of sizzled corresponds to 60% of maximum pSmad5 intensity (Fig 3D). While the expression boundary of foxi1 corresponds to 25% of maximum pSmad5 intensity (Fig 3E), that of bambia is very low at only 7% of pSmad5 maximum intensity (Fig 3F). The slope of the gradient changes across the DV axis in addition to the level of pSmad5. To determine if the expression boundaries correlate with distinct slopes of the pSmad5 gradient, we overlaid the boundaries of the target gene domains onto the slope of the pSmad5 gradient. The 3 target gene domains also correspond to 3 significantly distinct pSmad5 gradient slopes ( Fig 3C). The expression boundary of sizzled corresponds to a gradient slope of 1.

Test of gradient slope to position gene expression boundaries
The boundaries of target gene expression across the DV axis correspond to both distinct pSmad5 concentrations and pSmad5 gradient slopes. To directly test the ability of either pSmad5 concentration thresholds or gradient slope to position gene expression in the gastrula embryo, we utilized mutants that have a modified pSmad5 gradient shape and measured corresponding shifts in target gene expression (Fig 4A and 4B). If cells across the DV axis respond to distinct levels of pSmad5 to activate target gene expression, the boundaries of target gene expression will correlate with the same pSmad5 levels even if the gradient shape is altered ( Fig  4B). Alternatively, if cells respond to the shape of the gradient, the target gene boundaries will correlate with the same gradient slope regardless of pSmad5 level ( Fig 4B).
To investigate the spatial shifts in target gene expression when the pSmad5 gradient shape is altered, we used chordin mutant early gastrula embryos. Chordin is a BMP ligand antagonist that acts as a dorsal sink for BMP, thus shaping BMP signaling activity during gastrulation [23,35,63]. In chordin mutant embryos, the shape of the BMP signaling gradient is altered where the highest levels of BMP signaling activity expand laterally, and the slope of the gradient is shallower in the lateral regions ( Embryos from crosses between chordin-/-and +/-fish were immunostained at an early gastrula stage (7 hpf) for pSmad5, while siblings were assayed by FISH for the 3 target genes (S9C- S9E Fig). The pSmad5 gradient and FISH domains were quantitated blindly, followed by genotyping for the chordin mutation. The gradient of pSmad5 is expanded laterally in chordin mutants compared to the wild-type siblings (Fig 4C), as previously shown [23]. The expression domains of the 3 target genes were also significantly expanded laterally in chordin mutants ( Fig 4D-4F, S9B Fig).
To test whether a similar pSmad5 level delineated the boundary of the expression domains, possibly acting as a concentration threshold to provide positional information to cells, we determined the position of the target gene expression boundary on the pSmad5 gradients of wild-type and chordin mutant siblings. We found that a similar pSmad5 level corresponds to the boundary of sizzled, foxi1, and bambia expression in both wild-type and chordin mutants ( Fig 4G-4I, S9F Fig). The sizzled boundary corresponds to 56.8% and 56.5% of the maximum pSmad5 intensity in wild-type and chordin mutants, respectively. The foxi1 boundary corresponds to 16.5% and 20.8% in wild-type and chordin mutants. The bambia boundary corresponds to 7.2% and 8.0% in wild-type and chordin mutants. The similar levels of pSmad5 delineating the expression boundaries of these 3 target genes in wild-type and chordin mutants provides strong support for a concentration threshold model but does not eliminate the other models.
We also determined the pSmad5 gradient slope at the boundaries of the 3 target genes in both wild-type and chordin mutants. We did not find a consistent pSmad5 gradient slope at the expression boundaries for sizzled and bambia (Fig 4J and 4L

Test of signal duration to position gene expression boundaries
The gradient of BMP signaling forms from mid-blastula to early gastrula stages [22,23,36]. Embryonic cells are exposed to the BMP2/7 ligand for over 4 hours during this time period. It is not known if the activation of target genes requires differential duration to BMP signaling activity or if cells activate target gene expression once the pSmad5 threshold level is reached (Fig 5A). In a signal duration model that determines ventral target gene expression boundaries, the most ventrally restricted genes would require the longest signal duration (e.g., sizzled), whereas the most broadly expressed ventrolateral target genes would require the shortest signal duration (e.g., bambia) ( Fig 5A). To address the role of signal duration to pattern ventral cell fates, we tested the requirement of BMP ligand exposure to activate target gene expression. If cells respond to different durations of signal, then genes that require longer signal exposure will not be expressed after a pulse of BMP signaling. If cells respond to concentration thresholds, then a pulse of a high BMP2/7 concentration will activate all 3 target genes (Fig 5B).
To test the role of signal duration on the immediate response of gene expression, we injected BMP2/7 protein into bmp7 mutant embryos at 4 hpf and fixed embryos 30 minutes after injection for FISH and pSmad5 immunostaining (Fig 5C). We detected robust pSmad5  If target genes are activated by different pSmad5 levels, then all 3 target genes will be expressed following exposure to high levels of BMP signaling. If differences in signal duration activate BMP target gene expression, then a gene that requires a short signal duration will be expressed, but genes requiring longer signal durations will not be expressed. (C) Experimental schematic of a bmp7 mutant embryo injected with 5 pg of BMP2/7 protein that is fixed 30 minutes postinjection for pSmad5 immunostaining or FISH. (D) Representative immunostaining of pSmad5 intensities of an uninjected bmp7 mutant (n = 8) and a bmp7 mutant injected with 5 pg of BMP2/7 protein (n = 9). Animal pole is facing up. activation after 30 minutes of BMP2/7 ligand exposure ( Fig 5D). Strikingly, expression of sizzled (Fig 5E), foxi1 (Fig 5F), and bambia (Fig 5G) was also observed 30 minutes postinjection. This suggests that spatially distinct target genes can be rapidly activated following exposure to BMP ligand. Specifically, the activation of these target genes does not require the full 4 hours of BMP signaling that cells are exposed to endogenously.
While even shorter signal durations would be unlikely to be physiologically relevant, we nevertheless investigated gene expression responses at shorter durations of 10, 20, and 30 minutes in this assay. We found that sizzled is expressed within 10 minutes of activating BMP signaling in bmp7 mutant embryos (S10A and S10B Fig). The expression of foxi1 and bambia was first observed 20 and 30 minutes after BMP2/7 injection, respectively (S10C and S10D Fig).
While bambia with the broadest expression domain would be expected to require the shortest signal duration, we found that it was expressed latest among the 3 genes at 30 minutes in response to BMP signaling (S10 Fig).
To determine if genes in the same domain share the same transcriptional kinetics, we measured the expression of another broadly expressed target gene, ved (S11A Fig). The expression profiles of ved closely resemble bambia in both wild-type and chordin mutant embryos indicating that ved is activated by the same low pSmad5 threshold level (S11B Fig). However, unlike bambia, ved is rapidly activated in response to BMP signaling (S10D and S11C Figs). While differences in transcriptional kinetics have been suggested to underlie target genes activated by Nodal in zebrafish patterning [64], differences in the transcriptional kinetics of these BMP target genes do not correlate with domain size.
The mechanism of duration-dependent signaling can also include a genetic regulatory network that creates spatially distinct domains of expression [65]. To determine the role of secondary transcriptional regulation in defining the expression domains, we treated bmp7 mutant embryos with CHX at 4 hpf before injecting BMP2/7 protein (S12A Fig). Again, all 3 target genes were rapidly activated with the CHX treatment after 30 minutes of BMP ligand exposure (S12B-S12D Fig). Thus, BMP target genes can be expressed in distinct domains of the embryo independent of distinct durations of BMP signaling or feedback through genetic regulatory networks. Therefore, we conclude that different durations of BMP signaling activity are not directly positioning the expression of these target genes within the embryo.

Test of signal concentration to activate different target genes
We have shown that the 3 target genes are expressed rapidly upon exposure to high levels of BMP (Fig 5). If the target genes are responding to concentration thresholds alone, then exposing cells to different levels of BMP should activate target genes expressed in distinct domains regardless of signal duration or gradient shape. To expose cells to more stable and uniform levels of BMP, we manually disassociated cells from bmp7 mutant animal caps at 4 hpf and incubated the disassociated cells with 20 ng/ml or 5 ng/ml BMP2/7 protein for 2 hours. Incubation of 20 ng/ml and 5 ng/ml BMP2/7 recapitulated endogenous high and low levels of BMP signaling, respectively, found in wild-type embryos (Fig 6A, S13 Fig). Cells from bmp7 mutants expressed sizzled when incubated with the high level of BMP2/7 protein but not the lower level (Fig 6B, S14A Fig), as predicted in the concentration threshold model. Also consistent with the concentration threshold model, bambia is expressed in response to both high and low levels of BMP2/7 protein (Fig 6C, S14B Fig). Furthermore, these results show that a 2-hour duration of (E-G) Representative FISH in bmp7 mutants uninjected or injected with 5 pg of BMP2/7 protein for sizzled (E) (n = 11 uninjected, n = 11 injected), foxi1 (F) (n = 10, n = 11), and bambia (G) (n = 10, n = 11  a low BMP signaling level cannot induce the expression of a high-threshold gene that responds within 10 minutes in the embryo (S10B Fig). Thus, the activation of a high-threshold target gene does not display a duration-based response to BMP signaling.
We further quantitated in this assay system the number of cells above each predicted pSmad5 threshold for sizzled and bambia expression and the number of cells expressing each gene (Figs 3D-3F and 4G-4I). The predicted pSmad5 threshold for sizzled is 60 A.U., and the predicted pSmad5 threshold is 7 A.U. for bambia. There is a similar proportion of cells with pSmad5 levels above the predicted sizzled threshold (60 A.U.) and cells expressing sizzled in bmp7 mutants treated with 5 ng/ml and 20 ng/ml BMP2/7 (Fig 6D). Also, similar proportions of cells are above the predicted pSmad5 threshold level for bambia (7 A.U.) and are expressing bambia in the bmp7 mutants treated with 5 ng/ml and 20 ng/ml BMP2/7 (Fig 6E). Together, our data support a model where distinct concentration thresholds of BMP signaling activate spatially distinct target genes in DV axial patterning during gastrulation.

BMP morphogen interpretation: Multiple mechanisms for 1 signaling pathway
Multiple mechanisms have been reported for how a gradient of BMP signaling activity can provide positional information to pattern a tissue. In fact, there is conflicting evidence for how cells perceive gradients of BMP signaling depending on the context. In Drosophila, Dpp (the BMP homolog) is required for patterning and proliferation of the wing imaginal disc [66,67]. Cells in the wing disc have been suggested to sense the shape of the Dpp gradient [14] and the relative temporal change of Dpp signaling [68]. There is similar conflicting evidence for BMP patterning the dorsal neural tube. Both signal duration [10] and ligand identity [11] have been proposed to pattern neuronal identities. The mechanisms responsible for cellular interpretation of the BMP signaling gradient could vary in different contexts or by individual target genes. Understanding the mechanism of BMP morphogen interpretation in DV patterning will allow us to compare mechanisms across species and systems.
Here, we provide evidence that the BMP morphogen gradient acting to pattern DV axial tissues is interpreted in a concentration-dependent manner to activate 3 target genes underlying ventral cell fate specification. We precisely measured endogenous pSmad5 levels and gene expression activation within the embryo at a single-cell resolution. We found that the genes directly patterned by the BMP signaling gradient are expressed in at least 3 distinct domains of the embryo, which correspond to at least 3 different pSmad5 gradient levels that are required to activate representative target genes from each domain. Further, these pSmad5 levels act as thresholds that cells must reach to activate target gene expression and position gene expression boundaries in the embryo, regardless of gradient shape or exposure time. Together, our data support a model whereby the BMP signaling gradient is interpreted as a concentration-dependent morphogen providing positional information to pattern gene expression along the embryonic DV axis.

Multiple expression domains directly patterned by the BMP gradient
The genes directly patterned by the BMP signaling gradient remained unknown prior to this work. Identification of the genes reading out the BMP gradient during DV patterning is not only critical to study the mechanism of morphogen patterning we report here, but also valuable for investigations of the ventral cell types specified by this morphogen gradient. Known markers of ventrally derived tissues were found to be directly regulated by BMP signaling such as tp63 specifying epidermis [55,56], foxi1 specifying otic placode [37,42], and tfap2c specifying neural crest precursors [57], as well as genes with unknown roles in cell fate specification. We identified many genes encoding components of BMP and other signaling pathways. Interestingly, BMP signaling directly activates the canonical Wnt receptors fzd4/5, the Nodal signaling cofactor foxh1, and the Retinoic Acid binding protein crabp2b. This pathway specific feedback may be critical for robust gradient formation and integration of multiple signaling pathways during embryonic patterning. While our analysis identified the initial gene expression readout of the BMP gradient, further work is needed to resolve the genetic regulatory network committing progenitors to specific ventral cell fates.
To determine the number of positional values established by the BMP signaling gradient, we reanalyzed a scRNA-seq dataset [58] using the Seurat analysis algorithm [59] to infer the DV expression domains of the 57 direct targets identified. Based on these results, the target genes were sorted into 3 DV clusters and validated by analysis of a representative of each cluster as expressed in 3 distinct embryonic domains. Although we identified 3 distinct domains, the remaining target genes could further partition into additional domains across the DV axis. These broad overlapping domains will undergo further refinement through regulatory feedback and interaction with other signaling pathways later in development to specify cell fates. For example, the BMP target gene and preplacodal ectodermal marker dlx3b is initially broadly expressed in cluster 3 (S4 Fig). By late gastrulation (9 hpf), expression of dlx3b is restricted by factors also induced by BMP signaling, tfap2a/c, foxi1, and gata3 to form bilateral stripes that separate specification of epidermis and preplacodal ectoderm [37]. Future studies will have to address how other broadly expressed domains are refined into sharp, spatially discrete domains.

Concentration, not gradient shape or duration, positions expression boundaries
Our mutant analysis demonstrates that cells do not interpret the pSmad5 gradient slope to activate 3 spatially distinct target genes. The slope of the pSmad5 gradient undergoes a 2-fold decrease in chordin mutants compared to wild-type embryos in the ventral-lateral region (25 to 75 degrees) and a 2-fold increase in the dorsal-lateral region (125 to 155 degrees) (Fig 4J-4L), but the boundaries of target gene expression remain strongly correlated with specific pSmad5 levels (Fig 4G-4I). Previous analysis of the BMP signaling gradient over time has revealed that the gradient is highly dynamic from mid-blastula to mid-gastrula stages [22,23]. Nuclear pSmad5 levels rapidly increase at the most ventral positions (0 to 25 degrees) from 4.7 to 6.7 hpf. However, the pSmad5 gradient region determined to pattern target gene expression (70 to 110 degrees) is remarkably stable during this time [23]. Our lab previously found no significant difference in pSmad5 levels in this region during this time period, although the gradient slope undergoes a 2-fold increase [23]. Stable pSmad5 levels over time may allow cells to continually read out the gradient to reduce cell-to-cell variability, and steepening of the slope could allow for sharper gene expression boundaries to form.
We also showed that cells do not require integration of BMP signaling over a prolonged duration for initial patterning of the DV axis. Further genetic regulatory networks may act over time to refine the gene expression domains. We found that cells rapidly activate multiple, distinctly expressed target genes upon a 30-minute exposure to the BMP ligand. This is consistent with previous studies from our lab showing that BMP signaling prior to gastrulation is not required for DV patterning [22]. Specifically, reinitiating BMP signaling activity at 6 hpf in a BMP-deficient embryo rescues DV patterning. However, reinitiating BMP signaling activity at 6.5 hpf in a BMP-deficient embryo failed to rescue. Therefore, while cells do not require integration of BMP signaling before gastrulation, there is a critical window of time in early gastrulation during which cells are responding to BMP signaling. Furthermore, all cells across the DV axis are exposed to BMP signaling for the same length of time, while the gradient forms [22,23]. During mid-blastula stages (before 4 hpf), BMP signaling is activated at low levels everywhere, before signaling is cleared from the dorsal side and the gradient steepens [22]. Together, these data show that cells do not require differences in BMP signal duration to activate spatially distinct target genes.
While BMP signal duration or gradient slope has been shown to pattern tissues in other contexts, we find that distinct pSmad5 threshold levels pattern the embryonic DV axis. It will be important to investigate the molecular mechanism by which these thresholds activate target genes to compare BMP interpretation across contexts. The molecular mechanisms that establish different pSmad5 thresholds will be particularly interesting because BMP-responsive Smad transcription factors bind DNA with weak affinity [69,70]. Therefore, the classic model where differential affinity of transcription factors to regulatory elements produces spatially distinct target gene expression patterns alone cannot underlie BMP interpretation. To increase affinity and selectivity for DNA, Smad proteins bind to other high-affinity DNA-binding transcription factors [71]. Differential DNA-binding of Smad to target gene regulatory elements may be mediated by interactions with these cofactors. In Drosophila, the cofactor Zelda is required for BMP target genes to be expressed in the correct domain [72] and may represent such a Smad cofactor. However, a Zelda ortholog is not present in vertebrates, so other cofactors may play this role. Future studies on the association of Smad5 with cofactors or chromatin modifiers will be essential to further uncover the mechanism establishing concentrationdependent morphogen interpretation.

Zebrafish
Adult zebrafish (Danio rerio) were kept at 28˚C in a 13-hour light/11-hour dark cycle, and procedures were approved by the University of Pennsylvania Institutional Animal Care and Use Committee (IACUC). All zebrafish husbandry were performed in accordance with institutional, national ethical, and animal welfare guidelines. The embryos used for experiments were between 0 and 8 hpf, with some phenotypes tracked 1 to 2 days post fertilization. Embryos were collected and raised at 28˚C in E3 solution. In this study, sex/gender is not relevant since zebrafish sex determination takes place after 25 days post fertilization [73]. The zebrafish lines used were WT (TU) (RRID: ZIRC_ ZL57), chordin tt250 (RRID: ZDB-ALT-980413-523, ZIRC_ZL61), and bmp7a sb1aub (RRID: ZFIN_ZDB-ALT-050216-2, ZIRC_ZL1390). Adult chordin tt250 homozygous fish were generated by injecting chordin mRNA into 1-cell stage chordin-/-embryos to rescue the embryonic ventralization and then raised to adulthood. Adult bmp7a sb1aub homozygous fish were used to produce all bmp7a sb1aub homozygous embryos. Adult bmp7a sb1aub homozygotes were generated by injecting bmp7a mRNA into 1-cell stage bmp7a-/-embryos to rescue the embryonic dorsalization and then raised to adulthood.
Images were analyzed with a custom MATLAB algorithm to identify individual nuclei center points and extract pSmad5 intensities from within each nucleus [23,75], which were normalized based on the standard calibration bead intensity. Resulting embryos were aligned across the DV axis and conformed using Coherent Point Drift. Population means were generated after genotyping for wild-type and heterozygous sibling controls, since all imaging and analysis was performed blinded. Mean profiles were generated by averaging pSmad5 intensities of cells in a 30-μm band. Three-dimensional (3D) embryo-wide displays of mean pSmad5 were generated by projecting all nuclei on a sphere divided into 4,800 equilateral triangles and nuclei within each triangle averaged together. The pSmad5 gradient slopes were obtained by fitting a lowess fit to the average 3D data's spherical coordinates phi and theta using the "fit" function in MATLAB.

Fluorescent in situ hybridization, imaging, and quantification
Embryos were fixed in 4% paraformaldehyde at 4˚C and gradually dehydrated in methanol. Whole-mount in situ hybridizations were performed using DIG-labeled anti-sense RNA probes (made with labeling kit: Roche, 11175033910, Basel, Switzerland) to sizzled [76], foxi1 [37], and bambia. Probes were visualized with anti-DIG-Horseradish Peroxidase (Roche, 11207733910, Basel, Switzerland) developed with TSA Plus Cyanine 3 kits using a 1:50 dilution (Perkin Elmer, NEL744001KT, Waltham, Massachusetts, USA), and nuclei were stained with 1:1,000 dilution of Sytox Green. Embryos were cleared, mounted, and imaged as described for immunofluorescence. Images were analyzed using the same MATLAB algorithm, except fluorescent intensity was extracted from a 25-pixel sphere from the center point of each nucleus to include the cytoplasmic staining.

Cell disassociation cultures
At 4 hpf, bmp7 mutant embryos were dechorionated and placed into 1X Modified Barth's Saline (MBS) as previously described [77]. One hundred animal caps were removed with forceps by cutting the blastoderm at approximately 50% of its height, and the collected cells were diluted to a final concentration of 5 × 10 5 in 1X MBS containing Gentamicin (50 ug/ml; Gibco, Carlsbad, California, USA). The cells were disassociated by quickly vortexing, and 5 or 20 ng/ml of hBMP2-hBMP7 heterodimer (R&D Systems 2339-BM) was added to the tube for 2 hours before fixing with 4% paraformaldehyde. To perform immunofluorescence and FISH, cells were transferred to glass slides by cytospin at 750 rpm for 5 minutes. The pSmad5 immunostaining was performed and analyzed as described above with the slides mounted in Fluoromount-G (Southern Biotechnology Associates, Birmingham, Alabama, USA) for imaging. For FISH, cells were fixed in 4% paraformaldehyde for 10 minutes and allowed to air-dry for 1 hour before performing the FISH and analyzed as described above with the slides mounted in Vectashield mounting medium (Vector Labs, Burlingame, California, USA) for imaging. Fluorescent signal was normalized to the median of the top and bottom 5% of cells in wild-type early gastrula embryos imaged on slides.

RNA sequencing and analysis
Two replicates of 40 wild-type and bmp7 mutants were collected at early gastrula (shield, 6 hpf) and mid-gastrula (70% epiboly, 8 hpf). Three replicates of 50 bmp7 mutant embryos treated with CHX with or without injected BMP2/7 protein were collected 90 minutes after injection. Total RNA was extracted from dechorionated embryos with Trizol and purified with phenol-chloroform. Libraries were prepared with Illumina TruSeq stranded polyA-selection mRNA kit (Illumina, San Diego, California, USA) by the Next-Generation Sequencing Core at the University of Pennsylvania. Libraries were analyzed using the Agilent BioAnalyzer (Agilent, Santa Clara, California, USA) and Kapa Biosystems library quantitation kit (Roche, Basel, Switzerland) before sequencing on a HiSeq4000. Sequence reads were aligned to the zebrafish GRCz11 genome assembly with RNA-seq Unified Mapper (RUM) [78]. Differential expression was determined with EdgeR, and values were normalized to counts per million. GO analysis of BMP target genes was performed using the PANTHER classification system for the enrichment analysis (http://pantherdb.org/) and Fisher exact test with Bonferroni correction and P < 0.05 to determine terms that are statistically significant [79,80].

Seurat analysis of single-cell RNA sequencing dataset
Previously published scRNA sequencing from 6,100 individual cells dissociated from 75% epiboly embryos were used [58]. Cells with less than 2,000 genes sequenced were excluded from the analysis. Genes sequenced in fewer than 3 cells were also excluded from the analysis. Locations of the 47 landmark genes used are shown in [59]. Expression of target genes were mapped into bins using the Seurat package version 1.2 [59]. Data were normalized in Seurat.

Statistical analysis
Statistical tests were performed on GraphPad Prism software, and Student t tests (2 groups) or analysis of variance (3 groups) were performed. Error bars represent standard deviation. Figure legends indicate the number of n values for each analysis. Each experiment was performed at least 2 times.

Ethics statement
All animal procedures and protocols were performed in accordance with the approved IACUC protocols (#803105 and #804214) of the University of Pennsylvania and with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health.  Table D in S1 Data for underlying data. (C) Differential gene expression of WT and bmp7 mutants at early gastrulation (6 hpf) using RNAseq. Significantly up-regulated genes in WT compared to bmp7 mutants shown in red, and significantly down-regulated genes are shown in blue. All other genes are shown in gray. A subset of known BMP-dependent genes is highlighted. See Table E in S1 Data for underlying data.  Fig 1B, 1E and 1F and S1B, S1C, and S1E