The zebrafish Tuebingen strain was used as WT fish and for generating mutants. Meomesa, Mhwa, and Meomesa;Mhwa mutant embryos were genotyped and obtained as previously described (Xing et al., 2022). Embryos were maintained in Holtfreter’s water at 28.5°C. Developmental stages of the maternal mutant embryos were indirectly determined by observation of WT embryos (Kimmel et al., 1995), which were born at the same time and incubated under identical conditions. For the Nodal signaling inhibitor SB431542 (SB) treatment, one-cell stage embryos (10 min postfertilization) were incubated in Holfreter’s water with 1% DMSO (control) or 50 μM SB and harvested or observed at desired stages. All experiments were approved by Tsinghua University Animal Care and Use Committee.

Whole-mount high-resolution double fluorescence in situ hybridization (DFISH) was carried out using the protocol provided by J. Gage Crump lab, as described by Welten et al. (2006) and Zuniga et al. (2010), with two modifications: antibody concentrations were 1:500 anti-fluorescein-POD and 1:1000 anti-digoxigenin-POD (without preadsorbing the antibody with prehybridized zebrafish embryos), and incubation time with fluorescein tyramide was 3 h in the dark. The DIG-labeled ndr1 probe and fluorescein-labeled ndr2 were used for DFISH.

Embryos were embedded in 1% low-melting-point agarose and imaged via Zeiss light-sheet Z.1 microscopy using a W Plan-Apochromat 20× objective at 0.5× zoom. Embryos were positioned with lateral views, and two lateral images were acquired for each embryo: one at the brightest view, and the other, rotated 180°. Each lateral view image generated by online dual side fusion included ∼250 z-slices with 1 μm interval. Then, two lateral views were fused by Multiview Process using ZEN (2014 SP1, black edition). The final image file (a size of ∼3 GB) per embryo was manually rotated in Imaris X64 9.0 to acquire the maximum intensity projections of animal views with dorsal to the right (indicated by the enhanced signal in the dorsal margin of WT and Meomesa embryos), which were captured after reset of all channels for further measurements.

For quantification of FISH, the maximum intensity projections were opened in Fiji to manually draw a circular polygonal region of interest (ROI) covering the entire embryo margin to measure the raw intensity profile (in uncalibrated optical density values, OD) per embryo. The staining dots in the marginal region were dispersedly distributed with blank, which could be a favorable control for the background subtraction. So, the minimal intensity (the intensity of the blank) in the marginal region (in a smaller and rigorous ROI) of the embryo was removed from the total intensity. The number of embryos assessed in FISH with final total intensity of ndr1 or ndr2 is presented in Supplementary Table S7.

We assume the expression levels of ndr1 and ndr2 follow the linear regression model with maternal Eomesa, Hwa/β-catenin, Nodal/Smad2, and time: Y n d r 1 = α 1 X e o m e s a + α 2 X β − C a t e n i n + α 3 X S m a d 2 + α 4 X T i m e + α 0 , Y n d r 2 = β 1 X e o m e s a + β 2 X β − C a t e n i n + β 3 X S m a d 2 + β 4 X T i m e + β 0 .

To estimate the value of the coefficients α and β and the t-test of the coefficient, we performed a maximum likelihood estimation using the statsmodel library in Python (Seabold and Perktold, 2010).

For the analysis of the regression coefficients of Hwa/β-catenin and Nodal/Smad2 in ndr1, we assumed that the following regression relationship was satisfied: Y n d r 1 ( i ) = α 1 X e o m e s a ( i ) + α 2 X β − C a t e n i n ( i ) + α 3 X S m a d 2 ( i ) + α 0 , where i denotes the different time points.

Similarly, we assumed that the following regression relationship of maternal Eomesa and Nodal/Smad2 in ndr2 was satisfied: Y n d r 2 ( i ) = β 1 X e o m e s a ( i ) + β 2 X β − C a t e n i n ( i ) + β 3 X S m a d 2 ( i ) + β 0 , where i denotes the different time points.

It is assumed that we want to find the causal effect of the existence of maternal Eomesa, Hwa/β-catenin, Nodal/Smad2, and the value of time on the outcome of ndr1 and ndr2. For ndr1, to define the causal effect, we consider two situations: Situation 1 (real situation): Where the existence of Hwa/β-catenin, Nodal/Smad2, and the value of time are changed, and the expression of ndr1 is observed. Situation 2 (counterfactual situation): Where the existence of Hwa/β-catenin, Nodal/Smad2, and the value of time are unchanged (but the expression of ndr2 and the existence of Eomesa are consistent with real situation). The causal effect is the difference between the expression of ndr1 attained in the real situation and the counterfactual situation: E [ Y ( n d r 1 ) , r e a l , X β − C a t e n i n = 1 , X S m a d 2 = 1 ] − E [ Y ( n d r 1 ) , r e a l , X β − C a t e n i n = 0 , X S m a d 2 = 0 ] .

In other words, X causes Y if changing X leads to a change in Y, keeping everything else constant. Changing X while keeping everything else constant is called an intervention, and it is represented by a special notation, d o ( X ) .

Formally, causal effect is the magnitude by which Y is changed by a unit interventional change in X: E [ Y ( n d r 1 ) | d o ( X = 1 ) ] − E [ Y ( n d r 1 ) | d o ( X = 0 ) ] .

To estimate the target quantity of the given observed variable, a backdoor criterion was introduced. If all common causes of the existing X and the outcome Y are observed, then the backdoor criterion implies that the causal effect can be identified by conditioning on all the common causes: E [ Y ( n d r 1 ) | d o ( X = 1 ) ] = E w E [ [ Y ( n d r 1 ) | X = 1 ,   W = w ] ,

where W refers to the set of common causes (the expression of ndr2 and the existence of Eomesa). After the identification of the correct estimand for the target quantity based on the causal model and estimation of the target estimand, we performed a refutation test by adding a random common cause variable (the estimation method does not change its estimate after we add an independent random variable as a common cause to the dataset), replacing treatment with a random (placebo) variable (the estimated causal effect should go to zero when we replace the true treatment variable with an independent random variable), and removing a random subset of the data (the estimation method does not change its estimate after we add an independent random variable as a common cause to the dataset) using the DoWhy library in Python (Sharma and Kiciman, 2020).

ChIP assay was performed as previously described (Listerman et al., 2006) with minor modifications: 1) the embryos (1,000 embryos per sample) were cross-linked with 1.85% formaldehyde (Amresco, 0493) for 15 min; 2) Sepharose beads were replaced by magnetic beads (Sigma, 16–663X) for convenience. The antibodies were mouse anti-Myc (Santa Cruz, sc-40, 10 μg), mouse IgG (Biyuntian, A7028, 10 μg, as control), rabbit anti-β-catenin (Cell Signaling Technology, #8408, 20 μl), rabbit anti-Smad2/3 (Cell Signaling Technology, #3102, 30 μl) (Liu et al., 2011), and rabbit IgG (Biyuntian, A7016, 0.125 μg, as control). Following ChIP, the purified DNA was used for qRT-PCR using specific primers (Supplementary Table S4, S5) with some changes according to Hong et al, (2011). Percentage inputs (% input) was calculated by “% Input = 2^(-ΔCt [normalized ChIP]),” where ΔCt [normalized ChIP] was acquired from “(Ct [ChIP] - (Ct [Input] -Log2 (input dilution factor))” with the input dilution factor of 10 in the ChIP. Standard deviation (S.D.) and multiple t-tests were analyzed using GraphPad Prism 7 from two independent experiments.

For making transgene constructs, the putative regulatory regions of ndr1 and ndr2 were amplified using specific primers (Supplementary Table S6), and the resulting fragments were ligated to drive GFP expression (Figures 3, 4). The transgene constructs were purified and injected into embryos at the one-cell stage. GFP expressions in embryos at 4.3–5 hpf were observed and photographed under Olympus MVX10 fluorescence microscopy.

Based on our previous study, the expressions of ndr1 and ndr2 are hypothesized to be mainly regulated by three factors, that is, maternal eomesa, maternal hwa-activated β-catenin signaling, and Nodal autoregulation (Xing et al., 2022). To quantitatively explore the dynamic contribution of these three factors to the expression of ndr1 and ndr2 during several time points, we first performed intensity compute of ndr1 and ndr2 in the maternal eomesa (Meomesa), maternal hwa (Mhwa), and wildtype (WT) embryos, as well as in the control group with the inhibition of Nodal autoregulation (SB431542 treatment). The expression levels of ndr1 and ndr2 were detected by the DFISH, imaged using light-sheet fluorescence microscopy (LSFM), and measured by the software Fuji (Figures 1A and B and Supplementary Table S7), which show a similar expression pattern to the previous reports (Meno et al., 1999; Xing et al., 2022) and simultaneously provide the quantified data for further modeling.

We applied “linear regression” and “causal inference” to evaluate the causal between the independent and dependent variables (Pearl, 2010). To reveal the relationship between variables, we constructed a multiple linear regression model with ndr1 and ndr2 intensities as dependent variables, respectively, and Eomesa, Hwa/β-catenin, Smad2, and time as independent variables. For the linear regression model of ndr1 (Supplementary Table S1), the p-value of the regression coefficient of the independent variable Eomesa >0.2 has no significant effect on the dependent variable ndr1, indicating that Eomesa may be the common cause for the outcome variable ndr1. For ndr2 (Supplementary Table S2), the p-value of the regression coefficient of the independent variable Hwa/β-catenin > 0.2 has no significant effect on the dependent variable ndr2, indicating that Hwa/β-catenin may be the common cause for the outcome variable ndr2.

Then, a causal graphical model was constructed for ndr1 with Hwa/β-catenin, Smad2, and time as treatment variables and ndr2 and Eomesa as the common causes (Figure 1C). The causal inference was performed using the backdoor criterion, and estimators were constructed using the linear model. Refutation tests by adding a random common cause variable, replacing treatment with a random (placebo) variable, and removing a random subset of the data were used to check its robustness to assumptions (Supplementary Table S3). It is of note that the “placebo_treatment_refuter = 0” indicates the model passes the refutation test and works. Similarly, the ndr2 causal graphical model was constructed with Eomesa, Smad2, and time as treatment variables and ndr1 and Hwa/β-catenin as the common causes (Figure 1H). Both causal graphical models pass the refutation test (Supplementary Table S3), which not only proves the correlation between ndr1, Hwa/β-catenin, and Nodal/Smad2 and the correlation between ndr2, Eomesa, and Nodal/Smad2, but also suggests that Hwa/β-catenin and Nodal/Smad2 can be the cause of ndr1 expression levels, while Eomesa and Nodal/Smad2 can be the cause of ndr2 expression levels.

We further calculated the variation of regression coefficients of Hwa/β-catenin and Nodal/Smad2 on ndr1 at different time points, and we found that the regression coefficients and the percentage contribution of Hwa/β-catenin gradually decreased during development (Figures 1F and G); meanwhile, those of Nodal/Smad2 gradually increased, indicating that the control of ndr1 expression was alternately transitional from Hwa/β-catenin to Nodal/Smad2. For ndr2, the coefficient variations of Eomesa continuously and steadily increased from 3.7 hpf to 4.7 hpf (Figures 1I and J), but the regression coefficients and the percentage contribution of Eomesa gradually decreased at later stages and were accompanied by the gentle rise of Nodal/Smad2.

To better understand how ndr1 and ndr2 are mechanistically regulated by Eomesa, β-catenin/Lef, and Smad2, the binding activity and regions were screened out by chromatin immunoprecipitation–quantitative PCR (ChIP-qPCR) analysis in zebrafish embryos at 4.3–5 hpf. Using immunoprecipitated DNA with a size range of 0.3–1 kb as the template, 14 ndr1 regions (ndr1R1–14) and 15 ndr2 regions (ndr2R1–15) (Figure 2), ranging from 60 to 250 bp, were amplified with specific primer pairs (Supplementary Table S4, S5). Results revealed that the promoter region (ndr1R11) of ndr1 was significantly bound by Eomesa (tagged with Myc), β-catenin, and Smad2, while Eomesa also occupied the upstream region ndr1R10 and Smad2 occupied the upstream regions ndr1R4, ndr1R5, ndr1R10 and the gene body regions ndr1R12 and ndr1R13 at the ndr1 locus. Concerning the ndr2 locus, the promoter region ndr2R10 was significantly bound only by Eomesa, and several upstream regions and gene body regions appear occupied by Eomesa and Smad2. However, β-catenin lacked binding peaks in the tested regions of ndr2, which is consistent with the observation that the ndr2 expression level was not significantly affected in Mhwa mutants (Figures 1A and B). It is of note that the ndr1 and ndr2 genomic DNA reported to be bound by Eomesa (Xu et al., 2014) are included in ndr1R11 and ndr2R13 (Figures 2A and B), and the ndr1R12 region covers the first intron of ndr1 predicted to be bound by Smad2 (Fan et al., 2007).

To test the sufficiency of the identified major regions screened out by ChIP-qPCR (Figure 3) for initiating and maintaining ndr1 expression, several genomic fragments amplified from each gene with specific primers (Supplementary Table S6) were used to construct GFP reporter transgenes (Figures 3A and B). The transgene constructs were individually injected into one-cell stage WT embryos, and GFP expression was observed at 4.3–5 hpf. For the ndr1 locus, the distal enhancer “a element” reported before (Fan et al., 2007), the 1,022-bp proximal promoter with 5′UTR (p), which embodies ndr1R11 with two putative Eomes-binding sites (EBS) and two putative Lef-binding sites (LBS), plus the 604-bp ndr1 intron 1 fragment (In1), which contains ndr1R12 with a putative Smad-binding site (SBS), were selected and fused to the gfp coding sequence to make the Tg(ndr1-aIn1p:gfp), Tg(ndr1-ap:gfp), and Tg(ndr1-In1p:gfp) construct (Figure 3B). Although Tg(ndr1-In1p:gfp) did not include the previously identified distal enhancer “a element” of ndr1 (Fan et al., 2007), it allowed GFP to express in the blastodermal margin with an expanded domain in the dorsal margin in over 70% of WT embryos (Figures 3B and C), the pattern of which mimicked endogenous ndr1 expression. The ratio of embryos with strong GFP expression in correct domains was reduced in Meomesa or Mhwa mutant embryos, while no embryos showed strong GFP expression in Meomesa;Mhwa double mutants (Figures 3C and D). The ratio of WT and single mutant embryos with strong GFP was decreased in the presence of SB. These results suggest that the ndr1 regulatory regions used here harbor the necessary cis-elements required for responses to Eomesa, β-catenin, and Nodal autoregulation.

For the ndr2 locus, we selected four regions (a–d), which contain Eomes-binding peaks (ndr2R3 in the “a” region, ndr2R4 and ndr2R5 in the “b” region, ndr2R10 in the “c” region, and ndr211 in the “d” region) or Smad2-binding peaks (ndr2R11 and ndr2R12 in the “d” region), for making constructs (Figures 4A and B). The basal construct Tg(ndr2-c:gfp) consisting of the “c” region (proximal promoter with 5′UTR) and the gfp coding region and exhibiting universal GFP expression in the ectoderm (data not shown) was used to construct Tg(ndr2-dc:gfp), Tg(ndr2-adc:gfp), Tg(ndr2-bdc:gfp), and Tg(ndr2-abdc:gfp). Initial test results identified that only Tg(ndr2-abdc:gfp) gave rise to GFP expression recapitulating endogenous ndr2 expression in the blastodermal margin but with nonspecific GFP in the ectoderm (Figure 4C). Then, we compared Tg(ndr2-abdc:gfp) expression in different genetic backgrounds and in the presence of SB (Figures 4D and E). The ratio of embryos with strong or mild GFP expression in the blastodermal margin in Mhwa was comparable to that in the WT background, which was unaffected by SB treatment. However, only a small fraction (32%) of Tg(ndr2-abdc:gfp)-injected Meomesa embryos had weak GFP expression in the blastodermal margin, and this ratio was slightly reduced (to 17%) in the presence of SB. When injected into Meomesa;Mhwa double mutants, GFP was hardly detectable. These changes were very similar to those seen for endogenous ndr2 expression. Therefore, the a–d regions of ndr2 harbor necessary cis-regulatory elements driving ndr2 expression in the blastodermal margin. It is worth noting that ndr2 regulatory elements for silencing ndr2 expression in the ectoderm appear missing in the tested regions.