Cell-to-cell heterogeneity in Sox2 and Brachyury expression guides progenitor destiny by controlling their movements

Although cell-to-cell heterogeneity in gene and protein expression has been widely documented within a given cell population, little is known about its potential biological functions. We addressed this issue by studying posterior progenitors, an embryonic cell population that is central to vertebrate posterior axis formation. These progenitors are able to maintain themselves within the posterior region of the embryo or to exit this region to participate in the formation of neural tube or paraxial mesoderm tissues. Posterior progenitors are known to co-express transcription factors related to neural and mesodermal lineages, e.g. Sox2 and Brachyury (Bra), respectively. In this study, we find that expression levels of Sox2 and Bra proteins display a high degree of variability among posterior progenitors of the quail embryo, therefore highlighting spatial heterogeneity of this cell population. By over-expression/down-regulation experiments and time-lapse imaging, we show that Sox2 and Bra are both involved in controlling progenitor motility, acting however in an opposite way: while Bra is necessary to progenitor motion, Sox2 tends to inhibit cell movement. Combining mathematical modeling and experimental approaches, we provide evidence that the spatial heterogeneity of posterior progenitors, with regards to their expression levels of Sox2 and Bra and thus to their motile properties, is fundamental to maintain a pool of resident progenitors while others segregate to contribute to tissue formation. As a whole, our work reveals that heterogeneity among a population of progenitor cells is critical to ensure robust multi-tissue morphogenesis.


Introduction
Cells are the functional units of living organisms. During embryogenesis, they divide and specify in multiple cell types that organize spatially into tissues and organs. Specification events take place under the influences of the cell's own history and of environmental clues. Over the last years, access to new technologies has revealed that embryonic cells often display an unappreciated level of heterogeneity. For instance, gene expression analysis suggest that, within the same embryonic tissue, cells which were thought to be either equivalent or different, are actually organized into a continuum of various specification states (1,2). The impact of this new level of complexity on morphogenesis has not been extensively explored due to the difficulty of experimentally manipulating expression levels within targeted populations of cells in vivo. The progenitor cells located at the posterior tip of the vertebrate embryo (called here posterior progenitors) constitute a great model to study how a population of stem cell-like cells being specified into distinct cell types. The use of fluorescent tracers in bird and mouse embryos revealed that these progenitors contribute to the formation of the presomitic mesoderm (PSM), the mesodermal tissue that gives rise to the muscle and vertebrae, and of the neural tube the neuro-ectodermal tissue that gives rise to the central nervous system (3)(4)(5)(6). These experiments have also shown that while some cells are leaving the progenitor zone, other progenitors remain resident in this area. Grafting experiments have later confirmed these properties by showing that these resident progenitors are indeed capable of self-renewal while giving progeny in different tissues (7,8). The posterior region therefore contains a population of different types of progenitors able to give progeny in the mesodermal or the neural lineages along the antero-posterior axis of the vertebrate embryo. Retrospective clonal analysis studies performed in the mouse embryo (9) confirmed heterogeneity of this progenitor population. Indeed, these studies revealed the existence of single progenitors able to give rise to only neural or mesodermal cells but also pointed out to the existence of a third type of progenitor able to generate both neural and mesodermal cells. These bi-potential progenitors, named neuro-mesodermal progenitors (NMPs), have later been shown to exist in zebrafish (10) and in bird embryos (48,50,51). In the early bird embryo (stage HH4-7), the future posterior progenitors are initially located in anterior epithelial structures: the epiblast and the primitive streak. At later stages (stage HH8 and onward), these progenitors are re-located caudally in the embryo in a dense mesenchymal structure that prefigure the embryonic tailbud where they will give rise to their progeny in the neural tube and PSM (48,(11)). To sustain the formation of tissues that compose the body axis, the posterior progenitor population must adopt an equilibrium between maintenance and specification/exit of the progenitor zone.
Two transcription factors, Sox2 (SRY sex determining region Y-box 2) and Bra (Brachyury), have been described for their roles in respectively neural and mesodermal specification during embryonic development (12,13). Sox2 is known to be expressed in the neural progenitors that form the neural tube where it contributes to the maintenance of their undifferentiated state. Its involvement in neural specification has also been revealed by a study showing that ectopic expression of Sox2 in cells of the presomitic mesoderm is sufficient to reprogram these cells, which then adopt a neural identity (14). The Bra protein was initially identified for its essential function in formation of the paraxial mesoderm during the posterior extension phase (12,15). Its crucial role in mesodermal specification has been demonstrated, in particular, by phenotypic study of chimeric mouse embryos composed of both Bra mutant and wild-type cells, and in which only wild-type cells are capable of generating mesodermal cells (16). More recent studies have shown that Sox2 and Bra proteins are expressed in posterior progenitors of developing embryos, indicating that activation of their expression takes place in progenitor cells before these cells colonize the neural tube or presomitic mesoderm (17,18). In addition, these studies have shown that both proteins are co-expressed in progenitor cells, an observation consistent with the presence of bi-potential progenitors in this tissue. Work done in mouse embryo and in vitro systems derived from embryonic stem cells indicates that Bra and Sox 2 influence the choice between neural and mesodermal lineages by their antagonistic activities on the regulation of neural and mesodermal gene expression (19).
The purpose of this study is to understand further the relations between specification and tissue morphogenesis processes within the progenitor region of the developing vertebrate embryo. In particular, we want to investigate which cellular mechanisms can underlie maintaining enough progenitor cells while others are contributing to the elongating paraxial mesoderm and neural tissues. Analyzing the relative level of protein expression, we show that the specification factors Sox2 and Bra are expressed with a high degree of spatial heterogeneity in the progenitor region of the quail embryo. Over-expression and down-regulation experiments show that the ratio of Sox2/Bra expression is key in maintaining progenitors posteriorly and in controlling the exit of progenitor cells in their destination tissues. Using time-lapse experiments, we show that progenitor cells display diffusive migration properties comparable to PSM tissue. We developed a mathematical model to explore how heterogeneous expression of Bra and Sox2 combined with their opposite action on non-directional motility can drive both progenitor maintenance and multi-tissue morphogenesis. By overexpression and downregulation in vivo experiments, we further validate that Sox2 inhibits progenitor motility whereas Bra promotes it. Finally, by using our modeling approach we propose that spatial heterogeneity plays an important role in maintaining progenitor tissue shape during axial elongation.

Results
Levels of Sox2 and Bra proteins display high spatial cell-to-cell variability in the posterior zone.
The transcription factors Sox2 and Bra are known to be co-expressed in cells of the progenitor zone (PZ) (17,18). As they specify from these progenitors, neural cells maintain Sox2 expression and downregulate Bra while mesodermal cells downregulate Sox2 and maintain Bra. Although Sox2 and Bra are recognized to be key players in driving neural and mesodermal cell fates, the spatial and temporal dynamics of these events remains to be elucidated. As a first step to address this question, we carefully examined expression levels of the two proteins in the PZ of the quail embryo stage HH10-11. As expected, analyses of immunodetection experiments revealed co-expression of Sox2 and Bra in nuclei of most, if not all, PZ cells ( Fig. 1 A-C) (n=8 embryos). Noticeably, we observed a high heterogeneity in the relative levels of Sox2 and Bra protein between neighboring cells. Indeed, in the progenitor zone, we found intermingled cells displaying high Sox2 (Sox2 high ) and low Bra (Bra low ) levels or, conversely, Bra high and Sox2 low levels as well as cells in which both proteins appear to be at equivalent levels. This cellular heterogeneity was very apparent when compared to the adjacent nascent tissues, i.e. the neural tube and the PSM, where Sox2 and Bra protein levels were found to be very homogenous between neighboring cells ( Fig. 1 D-F). We detected such cell-to-cell heterogeneity in progenitor as early as stage HH5-6, a stage corresponding to initial activation of Sox2 and Bra coexpression in the quail embryo (Supplemental Fig. 1). We also observed heterogeneous levels of Sox2 and Bra proteins in PZ cells of chicken embryo, indicating that it is not a specific feature of quail (Supplemental Fig. 2). To infer how Sox2 and Bra protein levels goes from being co-expressed in heterogeneous manner in the PZ to being expressed homogeneously in the nascent tissues we analyzed variations of their respective levels in a series of seven volumes (containing around 100 cells in each volume) located in a posterior to anterior path (from the PZ to the maturating tissues) corresponding to PSM or neural tube putative trajectories ( Figure 1G-G'). Data showed that expression level of Sox2 increases (+ 2.22 folds, n = 7 embryos) while that of Bra decreases (-3.81 folds) following the neural path ( Figure 1G). On the contrary, on the paraxial mesoderm path Sox2 level decreases (-2.12 folds, n =7 embryos) while Bra level first increases in the posterior PSM (1.14 fold, position 1 to 2) and decrease later on (-5.06 folds, position 2 to7) (Figure 1 G'). Next, to define whether the cellular heterogeneity found in the PZ depends more on variability of one of the two transcription factors, we quantified protein levels per nuclei of cells populating the PZ. By plotting Sox2 and Bra levels in individual cells we noticed a broader distribution for Sox2 levels (coefficient of variation of 41.8%) compared to Bra levels (coefficient of variation of 30.75%) (Fig 1 H), indicating that cell-to-cell heterogeneity in the PZ is preferentially driven by differences in Sox2 levels. To better quantify and visualize Sox2 and Bra heterogeneity we calculated the Sox2-to-Bra ratio (Sox2/Bra) for each cell of the PZ and compare it to neural tube and PSM. Comparison of Sox2/Bra values showed high divergences between the three tissues and confirmed the high heterogeneity previously observed in the PZ cells ( Fig 1I). It must however be noticed that these quantitative data revealed a broad range of cell distribution, highlighting, in particular, the presence of cells in the PZ displaying similar Sox2/Bra values as mesodermal or neural cells. We next asked whether the cellular heterogeneity caused by differences in Sox2 and Bra levels is present in the whole volume of the PZ or displays regionalization in this tissue. To do so, we analyzed the spatial distribution of the same Sox2/Bra values on optical transverse sections performed at anterior, mid and posterior positions of the PZ ( Figure 1J-J'''). This analysis confirmed that heterogeneous Sox2/Bra values are equally represented in the mid area of the PZ ( Figure 1J''). Cells with a high ratio (Sox2 High Bra Low ) were found to be more represented in the most dorso-anterior part of the PZ ( Figure 1J') and cells with a low ratio (Bra High Sox2 Low ) were found to more represented in the most posterior part of the PZ ( Figure 1J'''). This particular anteroposterior distribution was further confirmed by tissue expression analysis (Supplemental Fig. 3). However, variations of Sox2/Bra values were also noticed in these areas, indicating that cell-to-cell heterogeneity is present in the whole PZ.
Altogether, our data, highlighting significant variability in Sox2 and Bra protein levels within progenitors of the PZ, evidence an unexpected cell-to-cell heterogeneity of this cell population. Noticeably, despite an overall enrichment of Sox2 high cells in the dorsal-anterior part of the PZ and Bra high cells in the most posterior part, no clear spatial regionalization of these cells was detected, indicating that the PZ is composed of a complex mixture of cells displaying variable Sox2/Bra levels. This variability is further lost as cells enter the neural tube or the PSM.

Relative levels of Sox2 and Bra in PZ cells influence their future tissue distribution.
The fact that cell-to-cell heterogeneity caused by differences in the Sox2 and Bra levels is observed in PZ cells but not in PSM and the neural tube cells is suggestive of a role of these relative protein levels in the decision to leave or not the PZ and to locate in a specific tissue. To test this possibility, we developed functional experiments aiming at increasing or decreasing Sox2 and Bra levels in PZ cells.
To do this, we performed targeted electroporation of progenitors in the anterior primitive streak/epiblast of stage HH5 embryos to transfect expression vectors or morpholinos and further analyzed subsequent distribution of targeted cells, focusing on the PZ, the PSM and the neural tube ( Figure 2). As early as 7hrs after electroporation, we could detect the expected modifications of Sox2 or Bra expression in PZ cells for both overexpression and downregulation experiments (Supplemental Figure 4,5). We observed a significant decrease in the Sox2/Bra levels either by overexpressing Bra or downregulating Sox2 while this value increased when Sox2 is overexpressed or when Bra is downregulated (Figure 2 A,F). Consistent with previous studies (19), we found that cross-repressive activities of Sox2 and Bra contributed to amplify such ratio's modifications (Supplemental Figure 5). After transfection of expression vector or morpholino, we let the embryos develop until stage HH10-11 and examined fluorescent cell distribution in the different tissues. For this, we measured the fluorescence intensity of the reporter protein (GFP) in the PZ, the PSM and the neural tube and calculated the percentage of fluorescence in each tissue. We obtained reproducible data using control expression vector with less than 20% of the fluorescent signal found in the PZ (16.78% ±2.83), a little more than 20% in the PSM (22.64 % ± 3.30) and about 60% in the neural tube (60.57% ± 4.39) (Figure  2 B, E). We next found that overexpression of Bra leads to a marked reduction of the fluorescent signal in the PZ (1.17% ± 0.57) and to an increased signal in the PSM (33.04% ± 4.06) but has no effect on the neural tube signal (Figure 2 C, E). Elevating Bra levels is thus sufficient to trigger cell exit from the PZ and to drive cells to join the PSM. However, this is not sufficient to impede PZ cell contribution to form the neural tube. Similarly, we found that overexpression of Sox2 drives exit of the cells from the PZ (1.16% ± 0.67) favoring their localization in the neural tube (75.40% ± 4.57) without significantly affecting proportions of cells in the PSM (Figure 2 D,E). Spatial distribution of the fluorescent signals obtained using control morpholinos appeared very similar to that observed using the control expression vector (18.93% ± 3.06, 22.68 ± 4.09 and 58.38% ± 3.63 for the PZ, the PSM and the neural tube, respectively). We found that downregulation of Bra leads to exit of cells from the PZ (4.16% ± 1.57), favors cell localization in the neural tube (88.23%, ± 2.04) at the expense of the PSM (7.59% ± 1.81) ( These data thus point out that the relative levels of Sox2 and Bra proteins are key determinant of PZ cell choice to stay in the PZ or leave it to enter more mature tissues. Major changes of the Sox2 to Bra ratios, tending either towards higher or lower values, are sufficient to trigger cell exit from the PZ. These experiments also pointed to the critical influence of the relative levels of Sox2 and Bra in controlling the final destination of cells exiting the PZ, with Sox2 high (Bra low ) cells and Bra high (Sox2 low ) cells preferentially integrating the neural tube and the PSM, respectively. Intriguingly, we also observed that the preferential cell colonization of one given tissue is not always clearly associated with depletion of cells in the other (see discussion).

PZ cells are highly motile without strong directionality
To have a better understanding of how progenitors either stay or leave the PZ we needed to study their behavior dynamically using live imaging. We electroporated stage HH5 quail embryos with an empty vector encoding for nuclear GFP and preformed time-lapse imaging was further performed from stage HH8 to stage HH12. We focused our interest on the PZ but also on the PSM and the posterior neural tube in order to be able to compare migration properties between tissues ( Figure 3A, Supplemental Movie 1). Because these three tissues have a global movement directed posteriorly due to the embryonic elongation, we generated two types of cellular tracking: the raw movement, in which the last formed somite is set as a reference point and the "corrected" movement, in which the cellular movements are analyzed in reference to the region of interest ( Figure 3B). Tracking cell movements allowed for quantification of motility distribution, directionality of migration and time-averaged mean squared displacement ( Figure 3C-E) (n=7 embryos). First, we noticed that PZ cell raw motilities are higher than that of PSM or neural tube cell motilities ( Figure 3C, top panel). PZ raw directionality is also more pronounced in the posterior direction ( Figure 3D, upper panel). These results confirm that the PZ is moving faster in a posterior direction than surrounding tissues as previously measured using transgenic quails embryos (20). Analysis of local (corrected) motility reveals that PZ cells move globally as fast as PSM cells and significantly faster than neural tissue cells ( Figure 3C, bottom panel). The distribution of corrected PZ cell motilities is however different to the one of PSM cells as it showed slower moving cells (PZ corrected motility violin plot in Figure 3C is larger for slow values than the PSM counterpart and Supplemental Figure 6). After tissue correction, the directionality of PSM cell motion was found mostly non-directional as previously described (20) with a slight anterior tendency which is expected for a posteriorly moving reference ( Figure 3D, red plot in the lower panel). The distribution of corrected angles of PZ cell motilities is also globally non-directed with the exception of a slight tendency toward the anterior (to some extent more than the PSM cells), suggesting that our method is able to detect trajectories of cells exiting the region of interest ( Figure 3D, yellow plot lower panel). As PZ cell movement was found being mostly non-directional, we wanted to better characterized their diffusive motion by plotting their mean squared displacements (MSD), measured in each tissue over time, as it has been previously done for PSM cells (20). This analysis showed that progenitor MSD is linear after tissue subtraction, as intense as the PSM cell MSD and significantly higher than neural tube cell MSD, thus demonstrating the diffusive nature of PZ cell movements (Fig. 3E).
Together, these data evidenced that, in the referential of the progenitor region, PZ cell migration is diffusive/without displaying strong directionality (expect a slight anterior tendency), with an average motility that is comparable to PSM cells and significantly higher than that of neural tube cells.

Modeling spatial cellular heterogeneity and tissue morphogenesis
To understand better how heterogeneity in Sox2/Bra can control the choice of progenitors to stay or leave the PZ we designed an agent based mathematical model ( Figure 4). The main hypothesis of this model is that motility in a given progenitor is directly driven by the Sox2/Bra value; Bra is promoting non-directional motility whereas Sox2 is inhibiting it. We chose some initial conditions in which the neural tube, the PSM and the axial progenitor are already formed. This choice was motivated by the fact that our experiments focus on a developmental window, which is passed the formation of these tissues (HH8-HH12). At this period of development the main tissues deformations are in the anteroposterior and medio-lateral axis (21), thus we decided to design a 2D model (X,Y) evolving through time. We implemented cell numbers, proliferation, and tissue shape to be as close as possible to biological measurements (21) (Supplementary data). We set up the PZ cells to express random and dynamic Sox2/Bra levels with a defined probability to switch into a Sox2 high (Bra low ) state (neural tube state) or a Bra high (Sox2 low ) state (PSM state) ( Figure 4A). To take into account our biological observation that PZ cells are as motile as PSM cells we set up the threshold at which the motility is switching from high to low at a Sox2/Bra level of 1.6 (the level above which cells become committed to neural fate). Finally, to consider the physical boundaries existing between tissues, we integrated a non-mixing property between cell types. We first verified that our mathematical model (Supplemental Movie 2) recapitulates the basic properties of the biological system: in particular, it allows for recreating spatial heterogeneity of relative Sox2/Bra levels in cells of the PZ ( Figure 4B) and reproduces general trends in tissue motility and non-directionality ( Figure 4C,D). Simulations of our model also showed that relative cell numbers (taking into account proliferation) evolve as expected with a stable number of PZ cells and an increase in neural tube and PSM cells ( Figure 4E). To check if spatial heterogeneity of Sox2/Bra levels can self-organized in our model we made a simulation in which every progenitor start with equivalent levels of Sox2 (50%) and Bra (50%). We observed that spatial heterogeneity was emerging after few time-points and persisting all along the simulation suggesting that this feature is indeed able to self-organize independently of the initial levels of Sox2/Bra (Supplemental Movie 3). We next explored the model's ability to reproduce maintenance of progenitors and elongation of posterior tissues during the elongation process. By looking at different time points in the simulation, we observed that the progenitor region stays posterior to the neural tube. Neural tube and paraxial mesoderm both extend posteriorly ( Figure 4F, Supplemental Movie 2). These results show that our mathematical model is able to reproduce the main properties of our biological system. Interestingly, simulations also showed dynamic deformation of the progenitor region, which adopts asymmetric shapes and then gets back to symmetric, thus highlighting self-corrective properties of the system. To challenge further this model and test if it can recapitulate the experimental results we obtained by over-expression and down-regulation of function of Sox2 and Bra we explored the consequences of numerically deregulating the Sox2/Bra values on tissue and cell behavior. As a result, Bra High values increase PZ cell motility in the model ( Figure 4G), lead to more PSM cells ( Figure 4H The results of our numerical simulations are therefore coherent with our experimental data showing that Sox2/Bra levels control the balance between maintenance of progenitors in the PZ and continuous distribution of cells into the neural tube and the PSM. Together, the results obtained by our model strongly suggest that progenitor behavior can be guided by creating Sox2/Bra-dependent heterogeneous cell motility.

Sox2 and Bra control motility of PZ cells
Our mathematical model suggests that the control of cellular motility by Sox2/Bra levels is critical to influence progenitor behavior (staying or leaving the PZ to go to the neural tube or the PSM). To test this prediction we needed to define whether Sox2 and Bra are indeed involved in controlling progenitor motility in vivo. Therefore, we proceeded to over-expression and down-regulation experiments followed by time-lapse imaging in quail embryos ( Figure 5 A-F). We focused on monitoring the PZ to define how cells behave in this region, either staying or leaving this tissue. As for control embryos (Figure 3), we first monitored raw cell motilities (Supplemental Figure 7) and conducted subtraction of the tissue motion to gain insight into local motility and directionality ( Figure 5). We found that Braoverexpressing PZ cells display higher motility without significant differences in directionality when compared to control embryos. By contrast, when PZ cell overexpress Sox2, we detected a significant reduction of their motility compared to control accompanied by an anterior bias in angle distribution (Figure 5 B,C,D, Supplemental Movie 6). Bra down-regulation leads to similar significant reduction of cell motility, as well as a change in directionality toward the anterior direction 6. Conversely, Sox2 down-regulation did not result in significant effect on average cell motility or directionality, even though a tendency towards a slight increase in motility was noticed. (Figure 5 B,E,F, Supplemental Movie 7).
These data, showing that changing the respective levels of Sox2 and Bra is sufficient to modulate PZ cell motility/migration properties, highlight a key role for these transcription factors in controlling PZ cells movements with Sox2 and Bra inhibiting and promoting cell motility, respectively. In that sense, these results validate the importance of the regulation of motility by Sox2/Bra levels hypothesized by our mathematical model. When cells have high Sox2/Bra levels they migrate less and are left behind the PZ to be integrated in the neural tube. When cells have a low Sox2/Bra ratio they tend to migrate more, mostly in a diffusive manner, explaining how they leave the PZ to be integrated in the surrounding PSM tissues.
Modeling the importance of spatial heterogeneity in morphogenesis.
Our experimental and theoretical data point out that heterogeneity in Sox2 and Bra expression are key to control progenitor behaviors during axis elongation. However, we still do not know the importance of the spatial randomization of this heterogeneity (the fact that we observe very different Sox2/Bra levels in neighboring cells without clear pattern) in comparison to a heterogeneity that is spatially organized. To assess this particular point, we created a second mathematical model in which Sox2 and Bra levels, instead of being randomly distributed, are simply patterned in two opposite gradients as we have observed in the dorsal part of the PZ ( Figure 6A, Figure 1 J, and Supplemental Figure 3). In this area, Sox2 was found displaying an antero-posterior decreasing gradient while Bra displays an opposite antero-posterior increasing gradient. This new version of the model has been set up to reproduce similar dynamics in cell specification ( Figure 6B) and relationships between Sox2/Bra and motility as in the previous model (Supplemental methods). We observed that the motility of the PSM, and PZ of the gradient model are globally comparable to the spatially heterogeneous model ( Figure 6C, 4C) and mostly non-directional (Supplemental Figure 8). We observed maintenance of PZ cells caudally and elongation of the different tissues ( Figure 6D, Supplemental Movie 8) suggesting a gradient in Sox2 and Bra expression can also explain progenitor maintenance and distribution. Despite the observed similarities, we noticed several crucial differences between our heterogeneous and our biologically inspired gradient model. Indeed, the speed of elongation is less important in the graded simulation (0.8 a.u.) versus the spatially heterogeneous simulation (1.2 a.u.) ( Figure 6E). To check if, at the cellular level, resident progenitors are displaced differentially in the posterior direction between the two models we tracked cells at the center of the PZ and calculated the distances they travelled in the Y direction. Analysis showed resident progenitors move more posteriorly in the heterogeneous (mean distance of 3.37 a.u.) versus the gradient model (mean distance of 3.18 a.u.) suggesting that the heterogeneous model is more efficient at maintaining resident progenitor posteriorly ( Figure 6F). We also noticed less transient deformations and self-corrective behavior of the PZ in the gradient simulation when compared to the heterogeneous simulation (Supplemental Movie 9). By analyzing the PZ shape on longer time scales (beginning to end of the simulation) in each model, we found that, the shape of the PZ changes considerably in the gradient simulation. Indeed, it became larger (mediolateral) and shorter (antero-posterior), showing less conservation of proportions 41% (gradient) vs 86% (heterogeneous) ( Figure 6G). Finally, to test if the changes we observed could be due changes in the diffusivity of cellular migration we plotted the MSD through time for each of the models. We found that the MSD of the spatially heterogeneous model is higher that the gradient model ( Figure 6H) suggesting that the spatial heterogeneity in expression of Sox2 and Bra is enhancing the diffusive behavior of the PZ cells. This higher diffusivity can therefore bring more fluidity in the tissue in order to remodel more efficiently and to maintain its shape on longer time scales and allow for more posterior displacements of resident progenitors. Together, data obtained from our mathematical models argue in favor of spatial cell-to-cell heterogeneity being able to create an appropriate balance between two possible decisions for a PZ cell, i.e. staying in place or going away to contribute to adjacent tissues. Moreover, it argues that spatial cell-to-cell heterogeneity is important to allow tissue fluidity and remodeling in order to maintain PZ shape during axial elongation.

Discussion
Our results reveal the existence of spatial heterogeneity in expression levels of Sox2 and Bra proteins in the posterior progenitors of bird embryos. Experimental and theoretical approaches converge toward the point that spatial heterogeneity of Sox2/Bra ratio play essential roles in regulating progenitor maintenance and their allocation to the neural tube or PSM through the control of cellular motility. Altogether, our data lead us to propose the following working model. Progenitor cells expressing intermediate levels of Sox2 and Bra stay motile and remain resident within the posteriorly moving progenitor tissue. Cells expressing high levels of Sox2 reduce their motility and therefore are forced to leave the progenitor zone to integrate in the neural tube, while progenitors expressing high level of Bra "diffuse" more actively than other progenitors and exit the PZ to integrate the PSM. This model shed new light on how specification and morphogenesis is coupled during vertebrate embryo axis elongation and highlights the fact that heterogeneity can be a beneficial feature to ensure robustness in morphogenesis.
The fact that we can detect different levels of Sox2 and Bra expression between progenitors is indicative that axial progenitor cells are in different specification states, some progenitors are more engaged toward the mesodermal fate, some toward the neural fate while others are still in between these states. These different specification states could explain why in our gain and loss-of-function experiments preferential distribution of electroporated cells in the neural tube is not always paralleled by a decrease in participation in the PSM (or inversely) ( Figure 2E, J). Indeed a progenitor, which is for instance, already engaged toward a neural fate, might not be competent to become mesoderm anymore and might follow its path toward neural tube even if experiencing a sudden rise in Bra level. Single cell sequencing of the progenitor region of mouse embryos have revealed that different types of specification states co-exist within posterior progenitors (19). It is likely that the different states that we observe by visualizing different Sox2/Bra proteins ratio is also defined by differential gene expression of mesodermal and neural genes. To test this hypothesis further it would be interesting to analyze other neural and mesodermal genes and test if they have heterogenic pattern of expression in the progenitor region. Due to its technical limitation, single cell studies does not reveal the exact locations of the different progenitor states found within the posterior region. The first heterogeneity that we have observed is patterned spatially in a gradient along the antero-posteriorly axis (Sox2 high anteriorly, Bra high posteriorly) in the dorsal part of the region. This graded expression has been described in chicken embryo (22) (49) and is coherent with fate maps studies at earlier stages showing that the antero-posterior axis of the epiblast/streak gives rise to progeny along the medio-lateral axis (4,6). For instance, anterior cells expressing high levels of Sox2 can give rise to neural cells and more posterior cells expressing high levels of Bra to PSM (and eventually to lateral mesoderm to cells located even more caudally). However, our analysis also reveals a spatial heterogeneity between neighboring cells in the progenitor region (dorsally and ventrally). This finding is suggestive of a more complex picture where position in the progenitor region do not systematically prefigure final tissue destination. In this case, Sox2/ Bra ratio would be determinant in assigning progenitor fate independently of the initial position. This scenario involve much more cell mixing and could participate in explaining why prospective maps of this region have shown multi-tissue contribution (3)(4)(5)(6). In our functional experiments aiming at biasing PZ cells toward neural state (Pcig-Sox2, Bra-Mo) we have stronger effects in motility and tissue distribution than in experiments in which we aim at biasing PZ cells toward mesodermal fate (Pcig-Bra, Sox2 Mo). These differences can be explained and reinforced by our finding that the motilities of PZ cells are much more similar to PSM motilities than to neural tube. Therefore, a change of progenitor behavior toward a neural tube behavior is much likely to affect motility and progenitor distribution than a change toward a mesodermal state.
The action of graded signaling pathways such as Wnt, FGF and RA have been described to positively regulate the expression of Bra and Sox2 genes and to affect progenitor destiny (23)(24)(25)(26)(27)(28)(29)(30). It will be necessary and useful to test whether these pathways regulate Sox2 and Bra in posterior progenitors of quail embryos. However, no data allow us to say that the activity of these signaling pathways is sufficient to explain the spatial heterogeneity that we observed between neighboring cells in this region. Our data, as well as data from the literature, suggest the existence of a mutual repression mechanism of Sox2 and Bra in posterior progenitors (Supplemental Figure 5 and (19)). The expression of these transcription factors could therefore be controlled positively through signaling pathways and negatively by cross-inhibition activities. This synergy between signaling and cross-repressive activity could therefore allow for a temporal dynamic responsible for the spatial heterogeneity we observe. In our mathematical model, we propose that the expression of Sox2 and Bra is spatially and temporally dynamic within progenitors. This hypothesis needs to be investigated with specific reporter and live imaging. However if true, one could postulate that a progenitor that oscillate between a high Bra low Sox2 state and a low Bra high Sox2 state give rise to progeny in the neural or mesodermal lineages depending on these ratios at the moment of cell division. Studies in mouse and recent studies in birds suggest that a single progenitor clone (true NMP) can give rise to progenies in the neural tube and the paraxial mesoderm (9). Interestingly, the authors observe that the final positions of these progenies are different along the antero-posterior axis suggesting a sequential production: progenitor cells give rise to mesoderm and then switch to neuroectoderm (or the opposite). An oscillation between cellular specification states of the NMP could therefore be a possible explanation to these observations.
Recent studies have shown that during the course of axis elongation axial progenitor cells undergo an EMT before reaching their full potential and give rise to progeny in the neural tube and the paraxial mesoderm (23,31) (48). We observe that, if analyzed between stage HH5 and HH8, the progenitor region displays a low range of local movements in comparison to stages HH8 to HH11 (data not shown). Therefore, it is likely that local motility is very low when progenitors are still epithelial and become high and non-directional when the progenitors become more mesenchymal and are relocated caudally. However, it is interesting that even if the tissue is more mesenchymal at later stages it keep its global posteriorly directed tissue motion. Works in bird embryo have shown that the posterior movement of the progenitor zone is the result of physical constrains exerted by the neighboring tissues PSM and neural tube (21,32). PSM expansion can indeed exert pressure on axial tissues generate a general posterior movement of the progenitor zone. Our data indicate that the range of cellular motility of progenitors belonging to this moving region is key in determining their fate, high Bra expressing cells move actively and leave the region laterally whereas high Sox2 cells move less, and are left behind in the neural tissue. In our experiments, over-expression of Sox2 and down-regulation of Bra leads to an anterior bias in the direction of progenitor movements. However If we consider the posteriorly moving progenitor region, this would correspond to an absence of movement for neural progenitors. Indeed, it has been proposed that the progressive depletion of progenitors committed toward the neural fate occurs without excessive cell mixing (33) and therefore could not require active migration. Although the role of Sox2 on progenitor cell migration has, to our knowledge, never been described, recent work has demonstrated that Sox2 is implicated in the transition from neural progenitor to neural tube during chick embryo secondary neurulation (22). Concerning Brachyury, it has been previously shown that this transcription factor has a role in cell migration. Mouse cells that have a mutation in the Brachyury gene have lower migration speed than wild type cells when isolated and cultured, explaining part of the mouse embryonic phenotype (34). The molecular mechanisms that act downstream of Sox2 and Bra to promote or inhibit cell migration are still to be discovered.
From the result of our mathematical models, we can propose that both a spatially random and a patterned heterogeneity in Sox2 and Bra expression are able to maintain progenitors caudally and to guide their progeny in the neural tube and mesoderm. In our biological system we observe a superposition of spatially randomized and a patterned distribution of Sox2/Bra levels. It is therefore tempting to suggest that both systems (spatially random and gradient) are at work in the embryo. Interestingly, little is known about the role of spatial heterogeneity in morphogenesis. We propose that this random pattern allow for more cell rearrangements/tissue fluidity in the progenitor zone. This fluid like state and the opposite solid like state of the anterior PSM tissue have been shown to be key for zebrafish embryo axis elongation (35,36). Furthermore, the fact that we observe more self-correction in our model also suggests that spatial heterogeneity can provide plasticity to the system. Interestingly, several studies have shown that this particular region of the embryo is able to regenerate after deletion of some of its parts (37,38). In that regards, spatial heterogeneity could be easily re-organized in remaining cells in comparison to gradients. Indeed, if gradients of Sox2 and Bra are controlled by secreted signals, one could imagine that the absence of tissue could be more detrimental to the diffusion of this signal and therefore to re-patterning. Spatial heterogeneity in gene and protein expression is a common trait of living system and has been found in many contexts including cancer cells (39). The link between cellular spatial heterogeneity and the robustness of morphogenetic processes that we describe can therefore be relevant beyond the scope of developmental biology.

Quail Embryos and cultures:
Fertilized eggs of quail (Coturnix japonica) obtained from commercial sources, are incubated at 38°C at constant humidity and embryos are harvested at the desired stage of development. The early development of quail being comparable to chicken, embryonic stages are defined using quail (40) or chicken embryo development tables (41). Embryos are grown ex ovo using the EC (early chick) technique (42) 6 to 20 hours at 39°C in a humid atmosphere.

Electroporation :
We collected stage 4-6 quail embryos. The solution containing the morpholinos (1mM) and pCIG empty (1-2μg/μL) as a carrier or the DNA solution containing expression vectors Pcig , pCIG-Bra or pCIG-Sox2 (2-5μg/μL) were microinjected between the vitelline membrane and the epiblast at the anterior region of the primitive streak (45). The electrodes were positioned on either side of the embryo and five pulses of 5.2V, with a duration of 50ms, were carried out at a time interval of 200ms. The embryos were screened for fluorescence and morphology and kept in culture for up to 24 hours.

Image acquisition and processing:
Image acquisition for immunodetection was performed using a Zeiss 710 laser confocal microscope (20x and 40x objectives). Quantification of Sox2 and Bra levels in 3D have been made with Fiji or with the spot function (DAPI Staining) of Imaris. Expression have been normalized to DAPI to take into account loss in intensity due to depth, expression and ratios have been calculated and plotted using Matlab. Protein expression after gain and loss of function experiments has been done 7 hours after electroporation by analysing immunodetection signal levels within GFP positive progenitors and normalizing to endogenous expression. Fluorescence distribution in tissues has been acquired on a wide field microscope Axio-imager type 2 (Colibiri 8 multi-diode light source, 10X objective). The images of electroporated embryos were processed with the Zen software that allows the assembly of the different parts of the mosaic ("Stitch" function). The images were then processed with the "Stack focuser" plugin of the Image J software. The different tissues were delineated on ImageJ with the hands-free selection tool and then the images were binarized using the threshold tool. The total fluorescence intensity emitted for the cells transfected with the different constructs were measured and the sum of the positive pixels for the different tissues was calculated. The percentage distribution of fluorescence in the different tissues was then calculated.

Live imaging and cell tracking
Live Imaging has been done using Zeiss Axio-imager type 2 (10X objective) and was described here (31). Briefly, stage 7-8 electroporated embryos were cultured under the microscope at 38 degrees under humid atmosphere. Two channels (GFP and brightfield), 3 fields of views, 10 Z levels have been imaged every 6 minutes for each embryo (6 embryo per experiments). Images were stitched and pixels in focus were selected using Stack Focuser (ImageJ). X, Y drift was corrected using MultiStackReg adapted from TurboReg (ImageJ) (46). Image segmentation was done after background correction using Background Subtractor plugin (from the MOSAIC suite in ImageJ) and cell tracking was done using Particle Tracker 2D/3D plugin (ImageJ) (47). A reference point was defined for each frame at the last formed somite using manual tracking. Region of interest were defined manually and their posterior movement was defined by manual tracking of the tailbud movement. Subtraction of the tissue movement was done by defining the average motions of cells in the region. Violin plots were generated on Prism 8 (Graphpad). MSD and distribution of angles were calculated and plotted with a Matlab routine. Angle distribution was calculated from trajectories, weighted with velocities and plotted as rosewind plot using Matlab.

Mathematical modeling
We initially distribute 1100 progenitor cells, 1200 neural cells, and 3200 PSM cells, randomly in their respective areas. Each cell type is endowed with its proliferation rate: 11.49 hours for the progenitor cells, 10.83 hours for the neural cells and 8.75 hours for PSM cells. Each cell I is characterized by its ratio of Sox2/Bra R_I(t) between 0 and 1 (depicted as 0-2in Figure 4A to match with biological ratios) and its position in 2D (x_I(t),y_I(t)), each of these variables is time-dependent. In the heterogeneous case, we initially randomly attribute a ratio between 0.15 and 0.85 to the progenitor cells. At each time step, these cells update their ratio through a first order ODE using the function represented in Figure  5A (+ noise), and then update their position (x,y), depending on their ratio, by a biased/adapted random motion. Interaction properties between cells such as adhesion, maximum density, packing are implemented in the bias of the random motion, as detailed in the Supplementary Materials and Methods. In the simulations we represent a portion of the posterior body (1 unit=150 micrometers), block the movements of the cells in the most anterior region as we consider it a very dense area (somites, epithelium neural tube), and block their passage to either side of the PSM as we consider the lateral plate to be a solid structure. A-F : Sox2 (green) and Bra (red) expression has been analyzed at the cellular scale in the caudal part of stage HH11 quail embryo, either in the progenitor region (A-C) or in the nascent neural tube and PSM (D-F). Overlay images are presented in C and F. Note that cell-to-cell heterogeneity in Sox2 and Bra levels are apparent in the progenitor region, with neighboring cells having higher Bra (red arrow), higher Sox2 (green arrow) expression levels or both proteins at similar levels (yellow arrow) which is not apparent in the nascent neural tube and PSM tissues. G, G': Levels of Sox2 and Bra in different embryonic locations corresponding to neural tube maturation (G) and paraxial mesoderm maturation (G'). Images on the left display locations for the measurements (blue squares), charts on the right display the different levels measured. H: Distribution of normalized cell-to-cell expression for Sox2 and Bra in the progenitor region (n=8 embryos). Coefficient of variation are more important for Sox2 (41.48%) than for Bra (30.75%). I: Cell distribution of Sox2/Bra levels in the progenitor zone (n=9 embryos), the neural tube (n=7 embryos) and the PSM (n=8 embryos). J-J''': Representation of the Sox2 and Bra ratio (green to red) in digital transversal sections (40 µm)     A: Graphical representation of the mathematical function defining the ratio Sox2/Bra dynamics. In each progenitor cell, Sox2/Bra ratio oscillates randomly between 0.4 and 1.6, noise in the system ensures that some cells pass below 0.4 to specify into paraxial mesoderm (red) and some cells pass above 1.6 to become neural tube (green). Medium and low ratios (below 1.6) confer high-motility, high ratios (above 1.