Noggin4 is a long-range inhibitor of Wnt8 signalling that regulates head development in Xenopus laevis

Noggin4 is a Noggin family secreted protein whose molecular and physiological functions remain unknown. In this study, we demonstrate that in contrast to other Noggins, Xenopus laevis Noggin4 cannot antagonise BMP signalling; instead, it specifically binds to Wnt8 and inhibits the Wnt/β -catenin pathway. Live imaging demonstrated that Noggin4 diffusivity in embryonic tissues significantly exceeded that of other Noggins. Using the Fluorescence Recovery After Photobleaching (FRAP) assay and mathematical modelling, we directly estimated the affinity of Noggin4 for Wnt8 in living embryos and determined that Noggin4 fine-tune the Wnt8 posterior-to-anterior gradient. Our results suggest a role for Noggin4 as a unique, freely diffusing, long-range inhibitor of canonical Wnt signalling, thus explaining its ability to promote head development.

Movie S1. Comparison of EGFP-Noggin4 and EGFP-Noggin2 FRAP dynamics in the intercellular space of the animal ectoderm of the Xenopus laevis early gastrula. Intercellular spaces (IS) with EGFP-Noggin4 and EGFP-Noggin2 are shown. Bleached regions are visualised by white squares (15×15 and 7.5×7.5, respectively). In the case of EGFP-Noggin4, broader area was bleached purposely to make its quick diffusion more evident. Fluorescence recovery is demonstrated both by microscopy imaging and by normalised FRAP curves. A violet vertical line on every plot is moving in correspondence with the upper photo changes. In the first column we present weights calculated according to image analysis of Western and Flag-BAP calibration (Fig. S6a). Next, total weights were calculated taking into account the dilution ratio (column 2). At last the concentration in the precipitation reaction volume was computed taking into account the IP volume (300 µl) and Wnt8 molar weight. The adsorption isotherm at Fig. S6b was plotted using final concentrations in the right column. ; red asterisk -main axis; yellow asterisk -secondary axis). (i and i') By contrast, injections by the same way of even 150 ng/blastomere of Noggin4 mRNA did not induce secondary axes at all (0%, n=120; asterisk -main axis).Abbreviations: bl -blastopore, dbl -dorsal blastopore lip, n -notochord, np -neural plate, s -somite. Bar everywhere is 100 µ.   Figure S2. Noggin4 cannot antagonise BMP and Activin/Nodal signalling but inhibits Wnt/β -Catenin pathway. (a) Comparison of Noggin1 and Noggin4 abilities to bind BMP4, using immunoprecipitation. (b) In contrast to Noggin1 and 2, Noggin4 does not cause reduction of endogenous phosphorylation of Smad1. Embryos were injected with Noggin1/2 or Noggin4 mRNA (100pg/embryo) and phosphoSmad1 was detected on Western-blot by anti-phosphoSmad1 antibody (Cell Signaling Technology, Inc.). Overall level of Smad1 was detected by anti-Smad antibodies (ab66737, Abcam). Alpha-tubulin serving as a loading control was revealed by anti-tubulin antibodies (DM1A, Sigma). In all these cases, Fab fagments of Anti-Rabbit antibodies conjugated to alkaline-phospotase (A3937, Sigma) were used for detection of primary antibodies. (c) Effects of Noggin1, Noggin2 and Noggin4 on expression of the BMP-specific luciferase reporter, BRE. Embryos at two-cell stage were injected with BRE plasmid, either alone (control) or mixed with Noggin1 or Noggin2 mRNAs taken in a fixed concentration, or Noggin4 mRNA taken in increasing concentrations. (d and e) Effects of Noggin4 on expression of the Smad2-specific reporter ARE activated by ActivinB or Xnr2. (f) The effect of Noggin4 on expression of the Wnt/β -Catenin-specific reporter TOPFlash activated by Wnt8 and β -Catenin . (g) qRT-PCR analysis of expression of the direct genetic targets of the canonical Wnt pathway: Axin2, HoxA1, HoxB1, HoxD1, Siamois and Xnr3, in the late gastrula embryos (stage 12) injected at 4-cell stage with Noggin4 mRNA or Noggin4 MO1. For each of the tested genes, the expression level in the wild-type embryos was taken as one unit.  Figure S3. Effects of Noggin4 mRNA and MO2 injections on tadpole head structures and expression of anterior markers, FoxG1 and Rx, at the midneurula stage. (a-c) 5 days tadpoles injected as indicated.

9/28
Nog4 mRNA     The dorsal ectoderm explants were excised and collected from injected and wild-type, non-injected, sibling embryos at stage 11, followed by RNA extraction. qRT-PCR with primers for Noggin4 and EGFP was performed for samples obtained from non-injected and injected embryos, respectively. The data were normalised relative to expression of two housekeeping genes, ODC and EF1α as described in. 5 The expression level of the endogenous Noggin4 in the wild-type samples was used for normalisation and regarded as one unit.

20/28
Supplementary Materials and Methods

Calculation of the free diffusion coefficient on the basis of FRAP curves
The advantage of the Xenopus embryo model is that it allows to measure FRAP kinetics in single intercellular spaces between large flat ectodermal cells, avoiding technical complications. This model allows us to consider diffusion as a one-dimensional process, and record FRAP kinetics in a comparatively short time: less than 1 min. in the case of Noggin4. Accordingly, these characteristics of the Xenopus embryo model allowed us to use the simple equation of one-dimensional FRAP model described by Ellenberg et al. 6 The exact solution of this equation has a complex form and contains special functions. However, for most practical purposes, the following simple empirical equation that differs from the exact solution by less than 5%, can be used instead of the exact equation: 6 where I 0 is the background value of the fluorescence intensity (FI) in the middle point of a bleached region; I m is the background value of FI outside of the intercellular spaces; I(t) is a FI value in the middle-point at time t; w is the length of the bleached region; I(∞) is the FI value in the middle-point after complete recovery. If we denote the normalised concentration of the bleached substance by C z , then the equation (S1.1) transforms to: This equation was used to describe FRAP data sets for all studied proteins: EGFP-Noggin1, EGFP-Noggin2, EGFP-Noggin4, the secreted version of EGFP-hep and EGFP-Wnt8. All data were renormalised to obtain array data of C z (t) values according to the equation (S1.2). Then a running average with the frame size of four points was calculated. Whole timeline for a group of observations was divided into 40 intervals. For every interval the average and the standard deviation were computed. Data sets were fitted to the equation: where b is a major parameter that determines the diffusion coefficient (D = bw 2 /2π); a and c are parameters, which should neglect possible mistakes in normalization and calculation of dispersion at the very beginning of the FRAP curve. The following constraints were used: where s is a number of starting frames in the observation set (values 2, 3 or 4 were used). Limit-memory version of Broyden-Fletcher-Goldfarb-Shanno algorithm with Bounds (L-BFGS-B) was used for fitting. In particular we used BFGS functions included into scipy.optimize library. 7 The quality of fitting was estimated by the determination coefficient that was computed as followed: As shown in Fig S10, FRAP of EGFP-Noggin4 protein could be fitted by the equation (S1.2) with the determination coefficient 0.79. However, this equation produces much worse fitting for the data obtained for other proteins, including EGFP-Noggin1, EGFP-Noggin2, the secreted version of EGFP-hep and EGFP-Wnt8. When individual samples of a specific set were fitted separately (data not shown), a fluctuation of the parameter b (that is proportional to D) was measured as a coefficient of variance (CV ): While CV was small enough for EGFP-Noggin4, this parameter strongly fluctuates in the cases of other proteins. Thus, if CV for EGFP-Noggin4 was about 20%, its values for the rest of proteins were much larger, reaching 115% for EGFP-Wnt8. This indicates that not only diffusion, but also other physical processes could contribute to the recovery of bleached proteins.

Model of FRAP process in the case of diffusion with adsorption on immobile binding centers.
To model diffusion of a protein in the condition of adsorption, the following kinetic scheme was studied: where A is a rapidly diffusing protein, B is an immobile binding interactant, C is the complex of A bound to B; subscript "B" denotes a "bleached" protein (A can be bleached both in free and bound states). It is also assumed that both reactions have the same forward (k 1 ) and reverse (k −1 ) reaction rate constants. Let D be the diffusion coefficient of A. B is immobile and let it be distributed uniformly with initial concentration σ . C is also immobile. Let concentrations of A, B and C be a, b and c, respectively. Since the forward reactions have the second order, the reverse reactions have the first order and b = σ − c − c b , it follows that: Let the protein A be distributed uniformly and be in equilibrium with the immobile molecules (B, C) before bleaching. As bleaching and diffusion do not influence on a global equilibrium, then we have: 1. sums a + a b , c + c b are constants, only the ratio changes with bleaching/diffusion;

([A]+[A b ])[B]
[C] 3. the concentration of free binding centers σ − c − c b is also constant.
By denoting σ − c − c b by the constant σ f ree , we can reduce the system of equations (S1.4) and get: For this linear system, one can obtain an exact solution in Fourier form. If we replace a b = Ae λt e iδ x , c b = Ce λt e iδ x in the system (S1.5) we obtain: The eigenvalues λ 1,2 of the characteristic matrix have the form: Both eigenvalues are real and negative, therefore the homogenic stationary state is stable (the bleached region disappears). Assuming the concentration of the bleached protein at the reactor boundaries is constantly zero, we take the following Dirichlet boundary conditions: Let an initial distribution of the bleached protein (both free and bound one) be of the form: Therefore, we get δ = π/l. Finally, after finding eigenvectors for the system (S1.6), for a common solution of the system (S1.5) we obtain: (S1.9)

22/28
where C 1 , C 2 are determined from the initial maximum concentrations a 0 and c 0 defined in the equation (S1.8). By definition, the fluorescent recovery in the middle-point of the reactor is a sum: Combining the equations (S1.9) and (S1.10), we get: (S1.11) I(t) has a form of a linear combination of two exponents with characteristic times −1/λ 1 , −1/λ 2 (slow and fast, respectively). For these times one can obtain asymptotic estimations. Let represent equation (S1.7) in the following form: Under conditions of our experiment Dδ 2 ≈ 1 and k −1 << 1, so w < v. Next, we expand equation (S1.12) to a Taylor series by w: Taking into account, that we get: Using the assumption that w << v, we finally get for λ 2 : In the initial notations defined by equations (S1.12) for two characteristic times we have: Applying it to our system, we expect k −1 to be around 0.01 s −1 , K -around 10 −7 M, l -around 10 µm, σ f ree -around 10 −8 M and D -around 10 −8 cm 2 /s. As the concentrations and constants are in these ranges, the characteristic times defined in the equations (S1.13) have exact physical meaning. The first characteristic time is defined mostly by complex lifetime (1/k −1 ) and can be designated as "adsorption characteristic time". The second characteristic time is defined by spatial and diffusion characteristics (l 2 /D) and can be designated as "diffusion characteristic time". At large t values the solution is reduced to a single exponent with the large characteristic time. Under the experimental conditions "adsorption characteristic time" is much larger than "diffusion characteristic time". As a result, long-term characteristic time of FRAP kinetics reflects a lifetime of the complex.
Note, that without loss of generality, asymptotic equations (S1.11) can be used when boundary conditions differ from a simple sinusoidal form (S1.8). We do not prove it rigorously, however we have tested the equations (S1.12),(S1.13) on the set of numerical solutions of the equations (S1.5). These tests show that the equations (S1.12),(S1.13) well describe numerical data under different parameters and under an initial Gaussian distribution of the bleached proteins.
The rigorous solution of the equations (S1.5) could be found in the papers of Tardy and McGrath 8,9 but working with such a form is uncomfortable. In the first study 8 the equations (S1.5) described the diffusion process of actin monomers accompanied with actin polymerisation and similar conclusions were drawn. Two regimes were isolated: "diffusion regime", when the fluorescence kinetics is determined by diffusion, and "turnover regime" when kinetics is determined by the equilibrium between unbound and cross-linked actin. In our case two regimes with different characteristic times also could be selected.
The common trends observed for completely different systems significantly reinforce the result. Also one may note that the model similar to (S1.5) was used for the determination of binding parameters (theoretically) in the Carrero's review. 10 The most self-consistent theory of reaction-with-adsorption systems was developed in McNally's group. 11 Analysing the equation (S1.5) authors theoretically isolated three possible simplified regimes: "pure-diffusion dominant", "effective diffusion" and "reaction dominant". Our model objects (adsorbed morphogens) resemble the last one because the FRAP characteristic time is determined mostly by the dissociation rate constant.

Effective diffusion coefficient
To introduce an effective diffusion coefficient correctly, one may focus on a single time dependency of C z for EGFP-Noggin4 and the theoretical curve that is shown at Fig. 4b. The theoretical curve was plotted according to the equation (S1.2) using parameters w = 14, D = 1.1µm 2 /s. Obviously, the term 4πD/w 2 in the equation (S1.2) equals to the reverse time of (1 − 1/e) recovery (Fig. 4b). The value (1 − 1/e) equals approximately 2/3, so we denote this reverse time as τ 2/3 . When the curve could not be fitted by the equations (S1.1),(S1.2) correctly, then the effective diffusion coefficient could be introduced as: (S1.14) Suggest (e 2 − 1)/4π ≈ 1/2; then the simple form for the effective coefficient is: The effective diffusion coefficient determined by the equations (S1.15) appears to be equal to the real diffusion coefficient when diffusion is Fickian. Thus, this equation could be applied to an arbitrary FRAP kinetics. Accordingly, we used the equation (S1.15) to interpret FRAP kinetics when a diffusion was accompanied by adsorption. Notably, a similar equation was used to assess "effective" or "apparent" diffusion coefficients of proteins that interact with a cell membrane. 10 An observation error for effective diffusion coefficient calculated by the equation (S1.15) can be calculated by this: For our experiments, we used the value of 2 µm as observation error of bleached region size (∆w) and as an observation error of recovery time (∆τ).
To determine τ 2/3 we should fit the recovery kinetics of adsorbed proteins with a correct equation. As we show in Supplementary Materials and Methods, 1.2, the recovery kinetics for proteins that are adsorbed at the cell surface should have a two-exponential form. The equation (S1.11) could be rewritten using parameters of a fit (a, c, τ 1 , τ 2 ): Thus, we fit our FRAP data with this double-exponential curves using same L-BFGS algorithm with constraints. We considered τ 1 to be the diffusion-driven recovery time and it could not be greater than 50 s, whereas τ 2 was considered to be the adsorption-driven recovery time and it could not be less than 50 s. All data fitted by this way are presented at Fig. 5. This kind of a fit is much better for all proteins than a fit presented at Fig. S10, because the values of the determination coefficients are around 0.9 for all experiments. Notably, the double-exponential equation fits Noggin4 experimental points slightly better than Ellenberg's equation. This indicates that Noggin4 nevertheless demonstrates a weak binding to some components of the extracellular matrix, without a significant effect on protein diffusivity. Theoretical curves plotted over FRAP data (Fig. 5) were used for determination of τ 2/3 , whereas the size of a bleached zone, w, was set by FRAP procedure. These τ 2/3 and w values were used to compute D E .

Mathematical analysis of FRAP data
3D structure of human Noggin homodimer in a complex with BMP2 was described earlier. 12 To determine the shapes of Noggin2 and Noggin4 homodimers, we used x-ray coordinates of human Noggin (PDB ID: 1M4U; doi:10.2210/pdb1M4U/pdb) as a reference structure and the method of homological modelling based on MODELLER software to construct the 3D models of these homodimers. 13 To test the stability of homological models of Noggin proteins, all structures from this stage were simulated in the molecular dynamics. The coordinates of constructed homodimers were put into the boxes of size 10××14 nm and then the systems were solved with the sufficient amount of TIP3P water molecules and Na + /Cl − ion pairs were added to concentration of 150 mM. All-atom OPLS force field was used in topology. 14 MD calculations were performed in GROMACS package 15 with "stochastic dynamic" integrator and Parrinello-Rahman barostate. The step of 2 fs was used in simulations as well as constraints for all valence bond lenghts. The structure equilibrated after 100 ns of simulation was used in further calculations. As a result, the shapes of Noggin2 and Noggin4 were determined by using built models (Fig. S11).
Free diffusion coefficients for all proteins were calculated by means of Hydropro software, 16 which considers only the shape of a molecule but not intramolecular mobility. This type of estimation is rather rough and, thus, we did not consider detailed structure of the linker between EGFP and the core protein. In our model EGFP protein was simply placed near the N-terminus of every Noggin monomers in the dimer (see Fig. S11 left).
As we were unable to find any information about viscosity within the intercellular space of Xenopus laevis embryos, we performed calculation of diffusion coefficients at different values of viscosity. The obtained dependency of the diffusion coefficient on viscosity is linear in double-logarithmic coordinates (see, as an example of these calculations, the diagram for EGFP-Noggin4 on Fig. S11b). By means of this diagram, we determined that the diffusion coefficient of EGFP-Noggin4, which was measured with FRAP technique in a living embryo (7.1 µm 2 /s, see Table 1), corresponded to the diffusion coefficient of the same protein diffusing freely in liquid with the viscosity value of 0.05 Ps. This viscosity is very close to the viscosity of blood 3 and thus could be considered to be quite close to the viscosity within the intercellular space. In turn, this indicates that, by contrast to other Noggins, whose theoretic diffusivity exceeds their actual diffusivity in a living embryo more than in 40 times (see Table 1), the diffusion of Noggin4 in living embryo occurs almost freely, without significant interactions with HSPGs or other molecules in the intercellular space.

Estimation of the Noggin4-Wnt8 complex formation parameters by the in vivo imaging
Here we took advantage of the model described in Supplementary Materials and Methods, 1.2 to estimate the dissociation constant of the Noggin4-Wnt8 protein complex in a living embryo on the basis of FRAP data.
Let EGFP-Noggin4 be the rapidly diffusing protein A and let TagRFP-Wnt8 bound to HSPG be the immobile binding interactant B (S1.3). The idea of the developed approach is based on the fact that, as it is demonstrated above, the dissociation constant is determined by the ratio of amplitudes of the fast (diffusion, A 1 ) and the slow (adsorption, A 2 ) components of FRAP kinetics (see Fig. 4c). In particular, it results in different shapes of the theoretical curves built by means of the equations (S1.9) and describes FRAP kinetics of EGFP-Noggin4 in the presence of equal concentration of "immobile" TagRFP-Wnt8 at different dissociation constants of their protein complex (Fig. S13). In turn, having a set of experimental points of FRAP kinetics and fitting them to a curve built by means of the equations (S1.9), one may resolve an opposite task: to estimate the dissociation constant, which corresponds to this fitting curve.
The fluorescence intensity kinetics after photobleaching of the area where Noggin4 and Wnt8 are co-localised is shown at Fig. 9. This experimental data were fitted with the equations (S1.9). One can address the main system of equations (S1.5) and its initial conditions (S1.8) to observe all the parameters, which should be found. The concentration of Wnt8 and Noggin4 was measured from images and substituted into σ f ree , a 0 in these equations. We take D of free Noggin4 (see Table 1) as the free diffusion coefficient D in these equations. Then, the rest parameters k 1 , k −1 and their ratio (the dissociation constant) could be assessed by fitting. In Fig. S13 a series of curves which were plotted using our known parameters and different K values is shown. These experimental curves are similar to the red one in Fig. S13. Whereas the diffusion-driven kinetics (A 1 ) is much smaller than the adsorption-driven one (A 2 ), it is possible to make only an upper estimate of A 1 . It means that the value of K could be of the order of 10 −7 M or less. On Fig. S15 we show the fit of the averaged experimental curve with K = 10 −7 M. The fitting procedure was the same as described in Supplementary Materials and Methods, 1.1.

Estimation of the endogenous Noggin4 concentration.
At first, we assessed the concentration of the endogenous Noggin4 protein in the intracellular space of living embryos. To this end, we estimated with qRT-PCR by how many times the concentration of EGFP-Noggin4 mRNA in embryos that were analysed in our FRAP experiments was higher than the concentration of the endogenous Noggin4 mRNA. To perform this analysis, we collected five replicated groups of the animal cap ectodermal explants (10 explants in each group) excised from different batches of the early gastrula wild-type embryos and from their siblings preliminary injected with EGFP-Noggin4 mRNA in the same way as it was done for FRAP experiments. The ectodermal explants were taken given the fact that, unlike the wild-type embryos, their microinjected siblings may contain Noggin4 templates not only in the ectodermal layer but throughout the embryo. The relative concentrations of EGFP-Noggin4 and the endogenous Noggin4 mRNA in these groups of explants were compared by qRT-PCR with the following pairs of primers to EGFP and Noggin4 (these pairs of primers were preliminary selected as providing very close to identical efficiency of PCR from EGFP and Noggin4 templates): EGFP Forward 5'-GGCGACGTAAACGGCCACA and EGFP Reverse 5'-CTGCACGCCGTAGGTCAGG;

25/28
Nog4 Forward 5'-TGGTTGGTGCAGAGGGCCAG and Nog4 Reverse 5'-CATGCCAGGTGGCCAGGAGC As a result, we established that the concentration of the endogenous Noggin4 mRNA was on an average 20 times lower than the concentration of the synthetic EGFP-Noggin4 mRNA in embryos used for the FRAP experiments (Fig. S16).
Thus, assuming that the effectiveness of the endogenous Noggin4 mRNA translation could not be lower than that of the longer EGFP-Noggin4 mRNA, we concluded that the concentration of the endogenous Noggin4 protein could not be less than 1/20 of EGFP-Noggin4. Given that the concentration of the latter protein was estimated as 1 µM (see above), one may estimate the concentration of the endogenous Noggin4 protein in the intercellular space as 5·10 −8 M. For simplicity, we supposed the total concentration of Frizzled8 receptors that are able to bind Wnt8 to be close to that of Noggin4 -5·10 −8 M -throughout the neurectoderm. In contrast, concentration of Wnt8 forms a posterior-anterior gradient; in its posterior maximum being also close to 5·10 −8 M but declining to the anterior end of the embryo according to the logistic sigmoid function, which fits the experimentally revealed the anterior-posterior gradient of Wnt/β -catenin signaling in the presumptive neural plate.
Finally, we took into account that Noggin4 and Frizzled8 bind to Wnt8 with the constants of 10 −7 M (see above) and 10 −8 M, respectively.
Given these initial conditions, we analysed if the change in Noggin4 concentration could exert influence on the established concentration gradient of Wnt8-Frizzled8 complex in the embryo, whose overall length was taken as 1000 microns. To this end, the following mathematical model was developed. Let two simultaneous competing reactions exists: Let the concentrations of free Wnt8, Noggin4 and Frizzled8 be w, n and f , respectively. Let the full initial amount of every protein be known. Let denote the full amount by w 0 , n 0 and f 0 . So we get an equilibrium: For the full amount of Wnt8, we obtain: w + (n 0 − n) + ( f 0 − f ) = w 0 Combining last two equations we get: Finally, we have the system of two equations with two unknown variables: n and f . Thus, we can determine the amount of the Wnt8-Frizzled8 complex ( f 0 − f ), which is proportional to the intensity of Wnt8 signalling.
As it is indicated above, we assumed that the concentration of Wnt8 declines from 5·10 −8 M at the posterior end of the embryo to zero at the anterior end according to the logistic sigmoid function: where ∆x is gradient steepness (Fig. 7). Then, given that the

Calculation of K d by the theoretical analysis of the Western Blotting data
The model used for analysis of WB data was based on an assumption that the following reactions take place: