Introduction

Morphogens pattern tissue by acting over space in a concentration dependent manner (Wolpert 1969; Kicheva and Briscoe 2023). Biochemically, morphogen molecules may either regulate transcription directly (Driever and Nusslein-Volhard, 1989), or through a specific transducer following a cascade of intracellular interactions - a signaling pathway. The latter involves additional, potentially non-linear steps of signaling. A prominent example of morphogen patterning is the vertebrate ventral neural tube where the morphogen Sonic hedgehog (Shh), produced from the notochord and floor plate, forms an extracellular gradient. Through its transcriptional effector Gli, this Shh gradient specifies different neural progenitor domains along the ventral-dorsal axis (Jessell 2000, Figure 1A). In this system, Gli activity profiles (cellular responses) regulate the expression of multiple homeobox transcription factors (fate genes) (Briscoe et al., 2000). These fate genes regulate their own and each other’s expression as a genetic regulatory network (GRN); this network finally settles into one of several discrete and stable states that correspond to different neuronal fates (Lek et al., 2010; Balaskas et al., 2012). On the basis of population-level observations of the Shh pathway and transcriptional network, different aspects of the Shh signaling/Gli activity profile have been proposed to control cell fate choices: e.g., threshold levels (presumably by different affinities between Gli and fate gene enhancers, Ericson et al., 1997; Oosterveen et al., 2012), duration of response (sustained Gli levels under the adaptive negative feedback of the Patched [Ptc] receptor, Dessaud et al., 2007) and the time-integrated response level (Dessaud et al., 2010). As these response dynamics are integrated through the GRN towards fate decisions, inherent noise in gene expression (Rosenfeld et al., 2005) may generate heterogeneity in the response-fate correlation, despite the ability of the GRN to increase robustness (Balaskas et al., 2012). In addition, our previous work reveals extensive cell movement (positional noise) and a role of cell sorting in refining the boundaries of fate domains in the zebrafish neural tube after cells are specified (Xiong et al., 2013), which are mediated by an adhesion code downstream of cell fate choices (Tsai et al., 2020). These studies suggest a possible limit of precision in morphogen signal interpretation in fate specification. To test this possibility, single cell observations using live imaging of transgenic reporters to simultaneously monitor both morphogen response and cell fate choices in vivo are required (Xiong and Megason, 2015). The zebrafish neural tube, with established reporter lines of Shh signaling and neural progenitor fates (Huang et al., 2012; Xiong et al., 2013; Tsai et al., 2020), provides a feasible model system to perform this analysis.

Figure 1.

Results

Live tracking of individual neural progenitor fate choices in zebrafish embryos

To enable a direct, single-cell quantification of Shh response dynamics and fate choice, we developed an imaging protocol using double transgenic reporters along with a nuclear or membrane localized EBFP2 for cell tracking (Figure 1B). We combined the tgBAC(ptch2:kaede) reporter whose Kaede expression reports Shh response (Huang et al., 2012) and a second transgenic reporter (e.g., tg(nkx2.2a:mgfp) [Ng et al, 2005], tg(olig2:gfp) [Shin et al, 2003], Figure 1C) whose expression levels indicate different fate choices (Xiong et al., 2013). Movies with good spatial-temporal coverage of the neural plate/tube (e.g., Movies S1) were analyzed by single cell tracking in GoFigure2, an open-source software we previously developed for image analysis (Xiong et al., 2013). By manually placing a seed at the center of the same cell at each time-point in the dataset, we generated a series of spherical segmentations inside the cell, which capture the fluorescence intensity of each reporter and the position of the tracked cell over time. Together, this time-dependent positional and reporter intensity information represents what we call a ‘neural progenitor track’ (Figures 1D-E, Movie S2). We analyzed over 200 randomly chosen tracked cells to cover the ventral progenitor types in the zebrafish neural tube: the medial floor plate cells (MFPs), lateral floor plate cells (LFPs), motor neuron progenitors (pMNs) and “More Dorsal” cells (e.g., p2,p1,p0 progenitors) (Figure 1A; Xiong et al., 2013).

These cell fates form distinct domains highly conserved among vertebrates and prescribe the pattern of neurogenesis (Jessell 2000; Lewis and Eisen, 2003). In zebrafish, MFPs along with the notochord are the source cells of Shh; LFPs, pMNs and More Dorsal cells depend on Shh signaling for correct specification (Park et al., 2004; Huang et al., 2012). Taking advantage of the reporter expression and the stereotypic position at the end time point of each track, we identified and assigned different fates to the tracked neural progenitors (Xiong et al., 2013). We imaged cells from both the anterior spinal cord (approximately somite levels 4-9) as well as the posterior spinal cord (approximately somite levels 9-14). The shape, size and morphogenetic movements of both the Shh sending and receiving tissues vary along the AP axis (Xiong et al., 2013), thus allowing us to sample a diversity of input-output (i.e., Shh response-fate choice) relationships.

Heterogeneity is pervasive in Shh response dynamics and fate choice in single cells

Our set of ∼200 single cell tracks provide a first glimpse into the dynamics of Shh interpretation in vivo. The main new insights come from being able to observe the dynamics of two nodes of the GRN in the same live cell, namely a Shh response reporter (ptch2:kaede) and a fate reporter. To compare our data with previous studies, we first analyzed the average Shh response behavior of cells within each fate group (calculated by simple averaging of single cell tracks within each of the different fates). In both the anterior and posterior spinal cord in zebrafish (Movie S3 and S4, respectively), the average LFP shows the fastest increase of Kaede intensity, while the average pMN shows intermediate increase and the average More Dorsal cell is the slowest (Figures 2A,C). Accordingly, the olig2:gfp reporter is highly expressed in the average pMN but comes to a plateau after initial expression in the average LFP (Figures 2B,D). These data are consistent with observations that these cell types require different levels of Shh for specification (Ericson et al., 1997; Park et al., 2004) and that in amniotes the p3 progenitors (equivalent to LFPs in fish) initially express the pMN marker olig2 but then shift to express nkx2.2 with continued and/or increasing Shh exposure (Jeong and McMahon, 2005; Dessaud et al., 2007).

Figure 2.

We next ask whether single cells behave stereotypically like the averaged cell. Strikingly, at the single cell level, we observed high variability in Shh-signaling reporter dynamics across cells of the same fate and some overlap of cells that finally take on different fates (Figure 2A’-D’). There are examples of cells that have similar ptch2:kaede traces but distinct fate reporter traces (e.g., in Figure 2E, track 1 vs. 2, 3 vs. 4), as well as cells with very different ptch2:kaede traces but similar fate reporter traces (e.g., in Figure 2E, track 1 vs. 3). These data show that single cells have heterogeneity in the correlation between Shh response and fate choice. Cells that choose the same fate could have distinct Shh responses and cells that have very similar Shh responses may make different fate choices.

Quantifying reporter activity dynamics from fluorescent intensity tracks

Given the apparent heterogeneity and overlap in the single cell track data, are there certain features in the reporter dynamics that predict fate outcome more accurately? In other words, are the apparently very similar ptch2:kaede tracks (as seen in Figure 2E) different in a certain way that could explain their distinct fate choices? Our tracks offer us an opportunity to address this question by quantifying and systematically classifying Shh signaling activity, namely the transcription rate of the ptch2 reporter, from the fluorescence intensities. However, the intensity measurements must be handled carefully and interpreted with consideration of experimental limitations to provide useful insights.

First, while the fluorescence intensity levels measured in single cells by image segmentation are in principle proportional to the reporter protein concentration (Figure 3A), the acquisition and analysis processes may introduce several errors: bleaching, bleed-through, saturation, and segmentation measurement error. Our sampling frequency is even lower here than our previous work, in which other imaging settings were similar and quantification of fluorescence signals indicated that bleaching was not significant (Xiong et al., 2013), suggesting that the bleaching factor can be ignored here. To account for bleed-through, which refers to the phenomenon of signal from one fluorophore contributing to another during simultaneous acquisition, we applied a percentage subtraction from each channel after determining the bleed-through ratio using all data points (methods). To account for saturation, which refers to signal intensity exceeding the dynamic range of detection causing an artificial signal plateau, we removed data points where the fraction of saturated pixels exceeds 5% within segmented cells (methods). To handle noise that arises from random fluctuations of segmentation measurement between time points (technical errors), we applied a moving-average smoothing to the tracks over windows of 6 time points (36min, see Figure 6C for effects of varying this time window).

Figure 3.

This operation causes some short-time scale features of the tracks (which are otherwise not possible to analyze due to technical errors in the measurements) to be lost, but it preserves trends on longer time intervals (Figure 3B). We further confirmed the mitigation of the impact of technical errors by re-tracking the same cells, where the average variation of intensity measurement over 6 time points is consistently <5% (data not shown). Together this pipeline transforms raw tracks to an intensity track that reasonably approximates the reporter protein concentration dynamics in single cells.

Second, in order to estimate the cellular Shh response/Gli activity dynamics (reflected by the transcription dynamics at the ptch2:kaede reporter), we investigated the relationship between fluorescent intensity (I) and the transcription rate of kaede (x) (Elowitz and Leibler, 2000). Since Kaede and GFP proteins can be considered as stable on the time-scale of our experiments (Ando et al., 2002; data not shown), assuming a constant coefficient (c1) for the rate of Kaede translation and maturation as a function of kaede mRNA level (m), the intensity change rate of Kaede (Figure 3C) is given by:

m is in turn a function of transcription and mRNA stability, given in the form of half-life (t1/2):

Combining (1) and (2) we have a simplified relationship between measured fluorescence intensity and Shh response (represented by x):

c is a constant. Note that the fluorescent signal is time delayed due to lags caused by transcription, translation, and maturation of the fluorescent protein, therefore our raw tracks do not show reporter responses (transcription) in real time. Importantly, use of destabilized RNA or proteins would not reduce this lag, but only diminish reporter intensity thus making their measurement noisier. Furthermore, because of the time and detection sensitivity differences between Kaede and GFP, we are prevented from cross-correlating the two reporter channels on the same time axis. Despite these limitations, equation (3) still allows us to estimate the Shh response dynamics in each single cell and compare that between different cell fates within each dataset.

To validate equation (3) and estimate reporter mRNA half-life (t1/2) in vivo, we first applied a heat-shock pulse on transgenic embryos with a hsp70l:EGFP transgene (Figures 4A; Halloran et al., 2000, constructed with an SV40 polyadenylation sequence like the ptch2:kaede transgene). This causes a defined pulse of transcriptional output at the reporter (Rieger et al., 2005). According to equation (3), the GFP intensity change that follows should reflect the exponential decay of the gfp mRNA. Indeed, single cell tracks we generated provide an excellent fit with the model’s prediction (Figure 4B), allowing us to estimate the mRNA half-life to be around 1.7 hours. Next, we performed in situ hybridization of kaede after sequential chemical inhibition of Shh response using cyclopamine (thereby stopping new kaede mRNA synthesis) and compared the overall mRNA signal with controls for signal level over a series of time points (Figure 4C; Huang et al., 2012). This method suggests a ∼3-hour mRNA half-life. The two estimations are thus congruent on an in vivo half-life of fluorescent reporter mRNAs somewhere between 1-4 hours (effects of varying this parameter are assessed in Figure 7B). These data support the validity of equation (3) and enable us to use it and the estimated half-life to infer the dynamics of the ptch2 transcription rate (i.e. Shh signaling response) in our single neural progenitor tracks (Figure 3D). Note that as equation (3) contains second derivative with respect to time of the intensity, this fitting is prone to measurement noise and requires good time resolution. However, due to our smoothing operation, sharp and short-time response changes have been filtered out. Our subsequent comparisons only consider those features that span over at least 0.5-hour windows (e.g., the peaks and valleys seen in Figure 3D) which are comprised of at least 5 time-point measurements and describe trends that can readily be seen in the raw intensity track (Figure 3B). This compromise minimizes the impact of technical errors in the analysis at the expense of model resolution.

Figure 4.

Systematic analysis of Shh response heterogeneity and association with fate outcome

The processed tracks allow us to gain deeper insights into the characteristics and sources of heterogeneity. First we performed K-means clustering of ptch2:kaede dynamics in all cells (Figure 5A). The resultant clusters are then compared in terms of their contribution to the LFP, pMN and more Dorsal fates (Figure 5B). We found that large clusters all show fate heterogeneity between “adjacent” fates (LFP and pMN, pMN and more Dorsal), such as Clusters 1,2,3,5 and 7. In contrast, a given cluster does not contribute both to LFP and more Dorsal, with the only exception in Cluster 11. These data show that the heterogeneity between Shh responses and fate choices is pervasive across distinct types (i.e. clusters) of response profiles, and primarily presents itself around progenitor domain boundaries. Over a larger spatial scale (such as between LFP and more Dorsal), the response-fate relationship can be considered predictive and accurate at single cell level.

Figure 5.

To explore the sources of signal heterogeneity, we compared the average and variation of Shh response by fates in different datasets (Figure 5C). The pMN tracks showed the most variability in all datasets while LFP and more Dorsal tracks are more similar on average and less variable. A striking feature that stands out is the similarity between pMN and LFP tracks in the posterior dataset, where the average pMN response appears as strong as the LFPs. To account for the heterogeneity in timing of response initiation, we used model-fit response levels (Figure 3D) to identify the response initiation time and re-aligned the tracks by this time point (Figure 5D). The re-aligned pMN average tracks are more similar to each other in the anterior datasets and show reduced variability, suggesting that some of the variability can be attributed to differences in timing of reporter activation. The posterior pMN tracks remain distinctly heterogeneous, suggesting a major contribution of variability is associated with the anterior-posterior level of the neural tube, which carries differences in tissue size and shape, cell movement, and morphogen source, among others.

Correlation of key aspects of Shh response with fate outcome

To further compare the correlation of some key dynamic features of morphogen response to cell fate outcome between the anterior and posterior neural tube, we used the model-fit responses to calculate the maximum transient response level (the peak Shh response within the track), the average response level (average Shh response over 10-15hpf), and the response time (the amount of time with an above-basal Shh response) in the tracked time windows. To give a visual indication of how well each metric correlates with the graded fate outcome, we ranked single cell tracks by the value of each metric high-to-low to make a vertical “French flag”: the sharper the flag, the stronger the correlation. The sharpness is quantified by the percentage of correct fate prediction by the metric thresholds. A value near 50% (coin toss) would suggest the metric does not correlate with fate choice at all, whereas 100% would suggest the metric makes perfect prediction. For example, when we use two late-stage cell position metrics, the lateral-medial/dorsal-ventral (LMDV, this distance refers to the axis along which cells are being patterned which is initially lateral-medial because patterning begins at the neural plate which then gradually transitions to dorsal-ventral as the neural tube is formed) should do well in the flag correlation while the anterior-posterior (AP) position should do poorly since we are trying to predict fates across the DV axis (Figures 6A). Indeed, we find that the LMDV position is highly correlated with cell fates. For the anterior tracks, the LMDV distance is consistently >90% correlated while for the posterior tracks the LMDV distance is more poorly correlated at earlier stages but becomes better at later times consistent with a necessary role for cell sorting in sharpening fate domain boundaries over time (Xiong et al., 2013; Tsai et al., 2020). In contrast and as expected, the AP distances never correlate with fate choices in the anterior or posterior spinal cord (all around 50%). The diversity of tracks and tissue geometry both within and across the anterior and posterior regions of the spinal cord (Figure 6B) enables us to search for a common metric that best correlates the Shh response dynamics with fate choice.

Figure 6.

We compared the same sets of cell tracks by the maximum transient level, total time, and average levels of Shh response (Figure 7A), when ventral neural progenitors are known to become specified (Lewis and Eisen, 2003; Park et al., 2004; Xiong et al., 2013). For the anterior tracks, all metrics are >85% correlative with fate choices, consistent with previous population level studies in amniotes (Ericson et al., 1997; Dessaud et al., 2007; Dessaud et al., 2010).

Figure 7.

Interestingly, for the posterior tracks, only maximum transient level correlates with fate choice for >80% of the tracks and the other two metrics fail to make 70%. This difference of correlation percentage is robust to wide variations in values of parameters used in our quantification (e.g. estimated Kaede mRNA half-life, smoothing parameters used in intensity calculation and transcription activity fitting) (Figures 7B-C). Other metrics (except LMDV cell position at late times, which is one of our fate identifiers) do not correlate better than the maximum transient level when both anterior and posterior tracks are considered (Figure 7D, data not shown). No single metric makes perfect predictions.

Discussion

The imaging, cell tracking, and data analysis reported in this study open a fresh perspective for looking at the important classic problem of morphogen interpretation. While it is reassuring that our data (obtained with a distinct technique) on average corroborate numerous previous studies on Shh mediated ventral neural tube patterning, the increased resolution into single cells and wide temporal coverage from when cells start responding to Shh to cell fate-specification reveal a high degree of previously unrecognized heterogeneity. This single cell variability could be due to other factors that differ between cells and their cross-talk with the Shh gene regulatory network (GRN), such as Notch signaling (Huang et al., 2012), or just simply noise that arises in the GRN (Rosenfeld et al., 2005). These possibilities are difficult to distinguish with the existing data, but they pose a challenge to morphogen patterning in general: how to ensure pattern precision despite this high heterogeneity in signal responses?

One possibility is to code morphogen response in a robust way under the given heterogeneity. Through quantifying the transcriptional dynamics at the Shh reporter, we are able to directly assess the correlation between certain proposed information-coding schemes and fate outcome. In the anterior spinal cord, maximum transient level, response time, and average level are all predictive of fates. In the more posterior spinal cord however, our results show that the temporal dimension of Shh response (which is part of the response time and the average level) no longer explains cell fate choices as well as the maximum transient level. The cause of this difference might be that there are intrinsic or environmental differences between neural progenitors in the anterior vs. posterior spinal cord, which affect their morphogen interpretation. Supporting this idea, the tissue geometry, and the range and movement of cells in zebrafish neural tube vary greatly along the body axis (Xiong et al., 2013; this study). We also found that the posterior LFP and pMN tracks still show very similar Shh response dynamics even after response initiation time re-alignment, distinct from the tracks from anterior datasets, suggesting an overall noisier environment in the posterior neural tube in zebrafish. This may be because posterior neural progenitors undergo more cell mixing, or because the spatially small Shh gradient as a result of the posterior neural tube’s small size introduces more positional noise in signaling. Overall, simple response-fate models are insufficient to account for all the heterogeneities observed or to make precise fate predictions.

While it is tempting to construct more complex models using our single cell data, certain limitations remain that warrant caution. First, our live imaging is limited by the brightness of reporters and the sensitivity of detection. Dynamics that fall below our imaging capacity, especially the Shh response activities before and in the early hours of our tracks, are not as well resolved. Second, due to the constraint of a live embryo undergoing morphogenesis, certain cells cannot be fully tracked throughout the imaging time window. Third, even though most ventral cells are specified and no longer require Shh by the end of our tracks (17-18hpf, Huang et al., 2012; Xiong et al., 2013), some cells may remain unspecified and their fate reporter could further change beyond our detection. These limitations could give misleading results in overly fine-grained comparisons (e.g., response level within each 1-hour window, data not shown), where a much smaller number of tracks and time points could dominate and thus bias the analysis.

Nevertheless, the results presented here along with our previous study (Xiong et al., 2013) raise an important consideration for the understanding of morphogen interpretation: given the pervasive heterogeneity of single cells in position, timing, and signaling, how is patterning precision achieved? We speculate that precision at the single cell level may not be necessary as the heterogeneity can average out on the population level, and thus important patterning features such as the size of each progenitor domain can remain robust even under single cell noise. An additional morphogen gradient such as BMP from the dorsal neural tube could provide information in parallel with Shh to improve fate choice precision as a function of position (Zagorski et al., 2017), which is inherently challenging for a single morphogen gradient to convey (Zagorski et al., 2024). After cell fate choices were made, remaining positional errors can be corrected by cellular mechanisms such as cell sorting to refine the final pattern (Tsai et al., 2020). Finally, additional mechanisms such as apoptosis may provide means of correction. Further testing of these ideas would require improved imaging capacities and access to multiple nodes of the GRN with more accurate live reporters (e.g., fluorescent protein knock-ins).

Acknowledgements

We thank D. D’India for fish care, K. Mosaliganti, T. Tsai and members of Matlab community for sharing scripts, A. Green, B. Appel and J. Kuwada for transgenic lines, W. Ma, T. Mitchison, A. Schier, and Megason lab members for discussions. This work is supported by NIH R01 GM107733 to S.G.M and Natural Sciences and Engineering Research Council of Canada (NSERC) to P.H. F.X. also acknowledges the support of a UKRI-EPSRC Frontier Research Grant (EP/X023761/1, originally selected as an ERC Starting Grant).

Author Contributions

F.X., A.R.T. and S.G.M. conceived this study. F.X. and A.R.T. performed the experiments and analyzed the data. S.N. performed track clustering and variability analysis. T.W.H. performed hsp70l:EGFP experiments and contributed to data analysis. P.H. contributed data of ptch2:kaede line. F.X. and S.G.M. prepared the manuscript with input from all authors. The authors declare no competing financial interest.

Materials and methods

Zebrafish strains and maintenance

Tg(shh:gfp) (Shkumatava et al., 2004), tgBAC(ptch2:kaede) (Huang et al., 2012), tg(nkx2.2a:mgfp) (Ng et al., 2005), tg(olig2:gfp) (Shin et al., 2003), tg(mnx1:gfp) (Flanagan-Steet et al., 2005) and tg(hsp70l:EGFP) (Halloran et al., 2000) have been described.

Homozygote parents from two different lines were paired to produce double transgenic embryos for imaging. Natural spawning was used and time of fertilization was recorded at the single cell stage of each clutch. Embryos were kept in 28°C incubators/chambers during imaging and other times except room temperature during injections, dechorionating and mounting. Staging was recorded using morphological criteria and aligned to the normal table (Kimmel et al, 1995). All fish are housed in fully equipped and regularly maintained and inspected aquarium facilities.

Fish-related protocols have been approved by the Institutional Animal Care and Use Committee (IACUC) at Harvard Medical School.

Microinjections of mRNAs

A pMTB construct (Xiong et al., 2013) containing mem-EBFP2 or h2b-EBFP2 was used for mRNA synthesis with the mMESSAGE mMACHINE system (Ambion). For mosaic injections, one blastomere of 8 to 16-cell stage embryos was injected (Nanoject) with approximately 1nl of 20ng/μl mem-EBFP2 or h2b-EBFP2. One round of screening was applied immediately after mosaic injections to eliminate damaged embryos and/or ones that missed the injection (no retention of co-injected Phenol Red label). Before imaging, injected embryos were screened for health.

Timelapse confocal imaging

Live imaging was performed using a Zeiss 710 confocal microscope (objective: C-Apochromat 40X 1.2 NA) with a home-made heating chamber maintaining 28°C. Embryos were mounted using the dorsal mount (Megason 2009) with a stereoscope (Leica MZ12.5). The mounting steps for imaging the neural plate/tube have been described in detail in Xiong et al., 2013. 405nm laser was used for photo conversion of Kaede and excitation of EBFP2 at the same time. 488nm and 561nm were then used in a separate track for the GFP and Kaede(Red) signals (See also Figure 1B). A typical stack covers an imaging space of 303μm(x), 303μm(y) and 120μm(z), divided by 700X700X80 voxels and is taken within 6 minutes. A typical movie lasts ∼10 hours containing ∼100 stacks at 6 minute intervals. Data analysis

Track creation and fate assignment

For the present study 31 timelapse datasets were collected (29 contain ptch2:kaede with olig2:gfp:15; nkx2.2:mgfp:6; shh:gfp:4; mnx1:gfp:4. 2 contain olig2:dsRed and nkx2.2:mgfp). 21/31 allow cell tracking and 13/21 have high coverage. Tracks were generated for 7/13 and 4/7 were further processed. The datasets acquired in Zeiss .lsm format were first converted to Megacapture format for import into GoFigure2 software (www.gofigure2.org, see also Xiong et al., 2013; Xiong et al., 2014). Inside GoFigure2, segmentation and tracking were performed manually on each dataset on randomly selected ventral neural tube cells in the labeled embryo (e.g., Figures 1D-E). To calculate the cell’s positions in relation to the embryo, notochord segmentations were generated every 10 time points (1 hour in real time) for each presented dataset. Segmentation and track tables exported from GoFigure2 were further processed using custom scripts (provided in accompanying supplemental data files) to generate the single slice labeled movies (Movie S2). The movies were watched individually to determine the fate choice of the cell using a combinatorial criteria including fate reporter expression and final position in the ventral neural tube, as described in Xiong et al., 2013. First, at 16hpf, MFPs were identified as the centered, apically constricted, triangle-shaped column of cells immediately dorsal to the notochord (as illustrated in Figure 1A). This fate identification is checked with shh:gfp movies in which MFPs are GFP+. Next, LFPs are identified as the two columns of cells flanking both sides of the MFP (Huang et al., 2012). This fate identification is confirmed with nkx2.2:mgfp movies (e.g., Movie S1). Next, pMNs are identified as the cells dorsal to the floor plate cells that are strongly olig2:gfp positive. Finally, More dorsal cells are those that are further dorsal to the pMN domain and GFP-in olig2:gfp movies or that are located >40% ventral dorsal distance at 16hpf in movies without olig2:gfp signal. The segmentations, tracks and cell fate information were then assembled through a preprocessing pipeline in Matlab (Mathworks).

Track processing

The output of the pipeline is a set of Matlab matrixes containing organized raw data for downstream analysis (e.g., NewFinalWrkSpace_olig2_9.mat, see supplemental data files). This organized raw data was then processed through a series of filters (provided in accompanying supplemental data files) to account for multiple experimental factors that introduce small errors in intensity measurements including bleed-through, saturation, and noise (e.g. SaturationMask.m: a first filter that removes aberrations in intensity calculation arising from saturated pixels in the image. BleedthroughCorr.m: correcting bleed-through signal between two channels. Bleed-through correction runs automatically after the algorithm assesses the degree of bleed-through and then applies a percentage subtraction on the values to further reduce errors in intensity calculation. SmoothData.m: a smoothing filter that handles short-time fluctuations that come from segmentation variations and acquisition noise. Different utility scripts were coded to allow user interactions with both the raw and filtered data, listed in AnalysisScript_List.xlsx.

Track alignment

ptch2:kaede reporter tracks were aligned at the point at which the corresponding promoter activity traces crossed an activation threshold value and maintained it for three timepoints (18 min). The activation threshold was determined by comparing promoter activity values from more dorsally-fated cells, in which the reporter is generally inactive, and those from LFP- and pMN-fated cells, in which the reporter is activated. A threshold value of 0.025 could distinguish the two sets of traces, with only ∼60% of more dorsally-fated cells crossing the threshold at any point within the observation period compared to 100% of LFP- and pMN-fated cells.

Track clustering analysis

ptch2:kaede reporter traces were clustered using the k-means clustering function in MATLAB. Each track was truncated to the set of timepoints for which >80% of reporter traces had observations. After truncation, tracks were only included in the analysis if they had observations for >80% of remaining timepoints. After filtering, 182 out of 216 traces were included in the clustering analysis. The elbow method was used to determine an appropriate number of clusters. The number of clusters was varied from 5 to 25, and for each value, the total Within Cluster Sum of Squares (WCSS) value was calculated. For each cluster the WCSS value corresponds to the sum of square distances for each point in the cluster from the centroid of the cluster. The total WCSS value decreases with an increasing number of clusters (because individual clusters become tighter) and begins to plateau at approximately 15 clusters. This value was therefore used for the analysis.

Full analysis protocols, scripts and sample datasets are available for download here: https://github.com/xionglab/NT_heterogeneity