Tracking early mammalian organogenesis – prediction and validation of differentiation trajectories at whole organism scale

Single-cell transcriptomics has significantly contributed to our understanding of cell type diversity across species, organs and developmental processes. Efforts to build cell atlases from different model organisms include extensive transcriptomic profiling of embryonic development, with profiling of the mouse being particularly relevant given its broad use as a model for mammalian development. Independent efforts now provide coverage from embryonic days (E)3.5-E6.5, E6.5-E7.5, E4.5-E7.5, E6.5-E8.5, E6.5-E8.25, E9.5-E13.5 and E10.5-E15.0 (Argelaguet et al., 2019; Cao et al., 2019; Chan et al., 2019; Grosswendt et al., 2020; Han et al., 2018; He et al., 2020; Ibarra-Soria et al., 2018; Lescroart et al., 2018; Mittnenzweig et al., 2021; Mohammed et al., 2017; Pijuan-Sala et al., 2019; Scialdone et al., 2016), complemented by detailed analysis of specific organs such as the brain (La Manno et al., 2021) and the heart (de Soysa et al., 2019) or with emphasis on specific germ layers (Nowotschin et al., 2019). Combining datasets to provide an integrated timecourse over a longer timespan has been achieved by overcoming the challenge of integrating different sequencing technologies (Qiu et al., 2022). Of note, single-cell atlases of normal development have rapidly been used as a reference to interpret mutant phenotypes, providing new insights into the cellular and molecular processes controlled by key developmental regulators (Barile et al., 2021; Clark et al., 2022; Grosswendt et al., 2020; Guibentif et al., 2021; Mittnenzweig et al., 2021; Pijuan-Sala et al., 2019).

Integrating and annotating transcriptomic profiles in different contexts (e.g. across species, technologies, experiments, time points, etc.) represents a foundational step towards the construction and leveraging of cell atlases, but poses substantial challenges due to the difficulty of distinguishing between highly similar or transitional cell populations arising along complex differentiation trajectories. Without lineage tracing, computational inference of differentiation trajectories may provide useful information for understanding the dynamics of cellular diversification, but this needs to be interpreted with caution. The increasingly large number of methods for trajectory reconstruction furthermore highlights the need for benchmarking strategies (Saelens et al., 2019). Many trajectory inference methods are based on the concept of pseudotime, a latent dimension representing the transition between progenitors and differentiating cells as a function of transcriptional similarity. Though related, pseudotime and experimental time are different concepts. Thus, incorporating real time should improve reconstruction of developmental processes (Tritschler et al., 2019). Inspired by the Waddington landscape (Conrad Hal Waddington, 1957) and Optimal Transport theory (Monge, 1781), the probabilistic framework Waddington-OT (W-OT) takes advantage of experimental time and stochastic modelling to estimate the coupling probabilities of cells between consecutive time points and to reconstruct differentiation trajectories (Schiebinger et al., 2019).

Here, we report a densely sampled single-cell RNA-sequencing (scRNA-seq) atlas covering mouse development from E6.5 to E9.5 in 6 h intervals. This new atlas includes our 116,000 previously published E6.5-E8.5 transcriptomes that cover gastrulation and the initial phase of early organogenesis, complemented by 314,000 new E8.5-E9.5 transcriptomes that bridge a crucial gap in development not captured in existing datasets, i.e. the period of major morphological and organogenesis changes that occur between E8.5 and E9.5. This includes embryo turning, emergence of definitive-type haematopoietic cells, and initiation of the heartbeat and circulation. This new combined E6.5-E9.5 atlas delivers the most comprehensive transcriptome dataset for mammalian gastrulation and early organogenesis to date. To address the challenge of reconstructing and annotating cell lineages, we combined expert curation with a variety of computational methodologies to generate a resource of broad utility for the developmental biology community. We investigated further the intricate process of haemato-endothelial development, highlighting the multiple origins of endothelial cells and the formation of blood progenitors (BP) in asynchronous waves in spatially distinct sites. Moving beyond atlas generation, we also used precise embryo dissections to profile spatially defined areas of the developing embryo. Finally, computational cell fate predictions were contrasted with cell fates observed in state-of-the-art cell grafting experiments, where individual E7.5 primitive streak segments were orthotopically grafted into recipient embryos and resulting cell fates analysed after 24 h of culture by spatial and scRNA-seq analysis. Our study provides a blueprint for cell state/cell fate analysis during key stages of mammalian development, paving the way for future interdisciplinary studies of cell fates and origins.

We previously reported an scRNA-seq atlas covering mouse gastrulation and the early initiation of organogenesis between E6.5 and E8.5 (Pijuan-Sala et al., 2019). To capture the crucial organogenesis period between E8.5 and E9.5 we undertook a new timecourse experiment and integrated the new sampling time points (Fig. 1A). Thus, 116,312 cells distributed across nine time points from the original atlas were complemented with 314,027 new cells distributed across four new time points (E8.75-E9.5) as well as one overlapping time point (E8.5) to facilitate data integration (Fig. 1B). Combined, the new extended atlas, ranging from E6.5 to E9.5, contains 430,339 cells across 13 time points spanning 3 days of mouse development (Fig. 1A-C).

All embryos for the new dataset were dissected before droplet capture, to (1) profile yolk sac (YS) cells separately and (2) dissect the embryo proper to provide single-cell data anchored by anatomical location (Fig. 1A). As illustrated below, sequencing the various tissue segments independently aids the disambiguation of transcriptionally similar cell states (Fig. 1D-F). A combination of computational approaches and manual curation (see Materials and Methods, Fig. S1-S6 and Tables S1,S2) allowed us to define 88 major cell states, more than double the number identified in the previous atlas and reflecting the rapid diversification of cell states during the 24 h time window from E8.5 to E9.5 (Figs 1G and 2). The new combined dataset has been made freely available through a user-friendly web portal to be explored and leveraged by the wider scientific community (see Data availability).

The computational reconstruction of developmental processes using single-cell transcriptomics remains a major challenge despite the large number of available methods. The timecourse experimental design of this atlas provides the advantage of incorporating developmental stage information, thus permitting the use of methods such as W-OT that anchor inferred developmental progression in real time rather than pseudotime (Schiebinger et al., 2019). In W-OT, developmental processes are modelled using a probabilistic framework that allows inferring ancestor and descendant probabilities between cells sampled at consecutive time points, thus allowing for time series analysis. We applied W-OT to the developing haematopoietic and endothelial systems, which are needed early during development to enable effective circulation once the heart starts beating at around E8.25-E8.5. Blood and endothelium arise in waves, at least some of which are thought to entail shared haemato-endothelial progenitors (reviewed by de Bruijn and Dzierzak, 2017; Elsaid et al., 2020). The first blood cells are so-called primitive erythrocytes, arising at E7.5 from mesodermal thickenings in the prospective blood islands of the YS. Next is a wave of definitive-type BPs that originate from the YS vasculature from E8.5. The haematopoietic stem cell (HSC) lineage is the last to emerge from the endothelium of major arteries of the embryo, the dorsal aorta and vitelline and umbilical arteries, starting from E9.5. Both YS definitive-type BPs and HSCs are derived from a specialised subset of endothelium, the haemogenic endothelium, through a so-called endothelial-to-haematopoietic transition. The extended mouse atlas covers the two YS-derived waves of blood emergence, as well as a variety of intra- and extra-embryonic populations of endothelial cells, providing a unique opportunity to explore the highly complex emergence and diversification of the initial haemato-endothelial landscape during mammalian development.

To encompass all haemato-endothelial lineages, a W-OT fate matrix was computed for all blood and endothelial cell populations in the landscape [erythroid, megakaryocyte, megakaryocyte-erythroid progenitors (MEP), erythroid-myeloid progenitors (EMP), BPs, haemato-endothelial progenitors (HEP), endothelial populations from YS (YS EC), embryo proper (EP EC) and allantois (allantois EC), as well as venous endothelium and endocardium] leaving all other cells in the extended atlas grouped as ‘other fate’. In brief, a W-OT fate matrix is a transition probability matrix from all cells to several target cell sets at a given time point, commonly the end point of the timecourse experiment (here E9.25 and E9.5). Following this strategy, three main W-OT studies were performed. As a first proof-of-concept, only YS blood and endothelial cells were considered with E9.25-E9.5 cells from these two populations defined as targets (Fig. 3), which recovered the known divergence between the primitive and definitive-type YS waves over time. Specifically, the W-OT inferred primitive wave mainly generates nucleated primitive erythrocytes, as well as a small number of macrophage and megakaryocyte progenitors (red trajectory in Fig. 3B) (Tober et al., 2007). The second, definitive-type YS wave starts with the emergence of EMPs followed by MEPs (McGrath et al., 2015) arising from yolk sac haemogenic endothelium (YS HE) as highlighted by a collection of known markers for all these populations (blue trajectory in Fig. 3B, and gene marker inspection in Fig. 3C and Fig. S4-S6). To further characterise gene expression changes along the inferred primitive and definitive hematopoietic trajectories in the W-OT YS landscape, we employed Slingshot (Street et al., 2018) to infer lineage pseudotimes and used Tradeseq (Van den Berge et al., 2020) to identify gene changes associated with each trajectory (Fig. S7A,B and Table S3). It is worth noting that cells expressing lymphoid and microglial-like progenitor markers were also detected during the YS definitive wave, albeit at a low frequency (Fig. S8).

For the second W-OT analysis, we wanted to explore the potential origins of endothelial cells together with blood cells, and thus considered the complete haemato-endothelial progenitor population, which resulted in a considerably more complex haemato-endothelial landscape as endothelial cells are found in multiple regions across the embryo (Fig. 4A,B, Fig. S9). For our third W-OT study, we merged blood cells with ‘other fates’ in the W-OT matrix, allowing us to focus explicitly on the complex endothelial landscape (third W-OT analysis). It is worth noting that, due to high levels of transcriptional convergence at endpoints during blood and endothelial differentiation, it is technically challenging to interrogate differentiation towards specific terminal cell types. Therefore, we opted to include multiple terminal blood and endothelial cell states in our W-OT analyses, allowing us to explore differentiation trajectories that produce molecularly convergent cell states from transcriptionally distinct starting populations. Importantly, the comprehensive metadata associated with our dataset makes it straightforward for the wider research community to explore alternative data analysis options that may be better aligned with their research questions.

Altogether, our W-OT landscapes revealed at least three putative endothelial differentiation trajectories that occur in different spatial locations and time windows during embryogenesis (Figs S10A and S11A-C). The first wave of endothelial production initiates early in the YS (Figs S7C and S11C). By contrast, the second wave(s) occurs in posterior regions of the embryo at intermediate time points (Fig. S11B,E) and the final wave initiates at later stages and occurs in anterior/medial sections (Fig. S11A,D). Dynamic gene expression changes were explored with Tradeseq (Van den Berge et al., 2020) along these inferred trajectories (Figs S7C and S11D-G) and detailed results of these analyses are presented in Tables S3-S5. In addition to revealing multiple spatiotemporally separated trajectories with different gene expression patterns, these analyses revealed a strong connection between somitic tissues and EP ECs, highlighted by differences in the UMAP coloured by log odds between the first landscape and the second/third landscapes (Fig. 3A versus Fig. 4A and Fig. S10A). This connection led to the identification of the endotome, a cell population not previously reported in mammalian embryos, which will be explored in the following section.

Finally, to identify distinctions between downstream endothelial populations we performed differential expression analysis (Fig. S10C). As expected, YS endothelium expressed Lyve1, venous endothelium showed high expression of genes associated with vasculogenesis and angiogenesis such as Clec1b and Cldn5 and the endocardium expressed markers indicative of Bmp (Id1 and Id3) and Notch signalling (Hey1). By contrast, embryo proper endothelium was relatively immature, as indicated by lower levels of endothelial genes such as Cdh5, Pecam1 and Dlk1. Furthermore, embryo proper endothelial cells that arise in different anatomical locations expressed distinct sets of Hox genes (Fig. S10D).

A small but clear subset of endothelium appeared to be part of a continuous differentiation trajectory originating from cells transcriptionally similar to endotome, a somitic subset defined in zebrafish embryos (Nguyen et al., 2014) (Fig. 4B and Fig. S11A,D). Added credence is given to this transcriptionally defined trajectory because both the contributing endotome and endothelial cells are derived from the same portions of the embryo (anterior/medial regions), demonstrating the utility of sub dissecting embryos before generating cell suspensions for sequencing. Furthermore, the mouse endotome cells are characterised by the expression of marker genes including Cxcl12, Pax3, Meox1, Foxc2, Pdgfra and Alcam, consistent with their recent discovery in zebrafish (Murayama et al., 2023; Nguyen et al., 2014; Tani et al., 2020) as well as by Hlf expression which has been associated with intra-aortic haematopoietic clusters (Yokomizo et al., 2019) (Fig. 4C). Further characterisation of the mouse endotome-like population (Fig. S12 and Tables S6 and S7) reveals it also expresses genes related to skeletal system development (Snai1, Sox9, Tbx1, Tcf15, Twist1, Irx5, Pkdcc, Foxp1), vasculature development (Apoe, Nr2f2, Col3a1, Col4a1, Ednra), muscle formation (Fzd2, Six1, Tcf15, Twist1) and VEGF signalling (Vegfb, Vegfc).

In the zebrafish embryo, endotome cells migrate towards the dorsal aorta and differentiate into endothelial cells that contribute to the niche for the emerging HSCs (Nguyen et al., 2014). Similarly, earlier studies on chick embryos showed a somite-derived cell population that contributes non-haemogenic endothelium, replenishing vascular cells that have undergone an endothelial-to-haematopoietic transition (EHT) (Pardanaud et al., 1996; Pouget et al., 2006; Sato et al., 2008). More recently, angioblast-like cells with similarities to somitic mesoderm-derived angioblast and endotome cells of chicken and zebrafish embryos were also reported in pluripotent stem cell-derived human axioloids (Yamanaka et al., 2023). In lower vertebrate models, the endotome was shown to originate at the ventral–posterior region of the sclerotome (Tani et al., 2020), which agrees with the transcriptional neighbours for the endotome population we discovered here in mouse embryos. It took our densely sampled and sub dissected scRNA-seq approach to discover a likely connection between the newly discovered mouse endotome and intra-embryonic endothelium (Fig. 4A-C, Figs S10A and S11A).

Intriguingly, in the Force Atlas representation of the haemato-endothelial landscape a subset of endotome-like cells present in later stage embryos (E9.25-E9.5) are placed next to EMP (Fig. 4D, red box). In addition, these cells cluster with EMPs (Fig. 4Di, cluster 20) when performing Louvain clustering using the top 50 principal components from the haemato-endothelial landscape. However, this association is only visually highlighted when a force-directed layout is generated on a subset of the atlas. By contrast, clustering, cell type annotation and UMAP embedding over the entire atlas, as well as the endothelial only landscape, demonstrates that these cells have a transcriptional identity aligned with earlier stage endotome cells that are transcriptionally distinctive compared with EMPs (Fig. 4Dii versus Fig. 2 and Fig. S11A).

Further exploration of these later stage endotome cells reveals that they lack a clear haemogenic signature (Fig. 4E, Fig. S13) and instead likely represent progenitors of endothelium and/or vascular mural and connective tissue cells, so-called vascular associated cells (VACs; Fig. 4B-E, Figs S11A and S13). Our densely sampled and regionally sub dissected atlas therefore suggests substantial transcriptional plasticity and developmental potential within mesodermal cells involved in the formation of intra-embryonic blood vessels. Future studies will be required to investigate how this plasticity extends to the specification of intra-embryonic HSCs within the vascular niche at E10.5.

The complexity of this atlas is both an opportunity and a challenge for computational inference of developmental processes at the whole organism scale. Importantly, the atlas also provides a framework for complementary experiments which can in turn expand the impact of the atlas. Of particular interest to developmental biologists is the question of whether trajectory reconstruction analysis can unveil the timing at which cells start to differentiate towards either one or a set of specific cell fates. Although two-dimensional representations of single-cell expression data are often taken as a starting point, such representations lend themselves to overinterpretation, thus highlighting the need for experimental validation. We therefore complemented our extended atlas with (1) scRNA-seq of carefully dissected subregions of the primitive streak at E7.5 and (2) orthotopic transplantation of the same streak regions into recipient embryos followed by extended in vitro embryo culture and both microscopic as well as molecular fate analysis.

Previous primitive streak cell transplant and labelling experiments produced a fate map of cells along the anterior-posterior axis of the gastrulating mouse embryo revealing reproducible region-specific contribution to all the major cell lineages of the embryo (Kinder et al., 1999). To link these data to the extended mouse atlas, we revisited these experiments and dissected the primitive streak of E7.5 (early bud; EB) mouse embryos into four sequential regions labelled A to D (with A being most posterior and D being most anterior), generated single-cell suspensions and profiled the regionally dissected cells using a modified Smart-seq2 scRNA-seq protocol providing deep coverage of each cell (Fig. 5A, Fig. S14). We mapped these primitive streak cells onto the extended atlas and employed label transfer to assign cell types (Materials and Methods). Approximately 27% of the isolated cells mapped to pluripotent cell types such as epiblast or primitive streak, 46% mapped to mesoderm, 24% to ectoderm and 3% to endoderm lineages (Fig. 5B).

There were clear differences across primitive streak regions (Fig. 5C,E). For example, only cells from the posterior-most region A mapped to lateral plate mesoderm (LPM) and primordial germ cells (PGCs) as well as expressing known posterior markers such as Msx2 (Sun et al., 2016) (expressed in the allantois and known to be involved in PGC migration) and Bmp4 (Fig. 5C-E, Fig. S15). In contrast, cells from the most anterior region D mapped to the node and expressed several brain and cilia-associated genes, including Riiad1 and Pifo (also known as Cimap3) (Fig. 5E, Fig. S15). Cells from regions B and C were similar to each other and mapped to paraxial, presomitic and somitic mesoderm as well as neural and gut cells. However, these cell types were not unique to regions B and C and were found across the entire posterior-anterior axis of the primitive streak (Fig. 5C,D, Fig. S15). Consistently, Robo1, involved in central nervous system and heart development, Fzd10, associated with neural induction, Prickle1, a limb development gene, and Gas1, a somitic gene also present in LPM, all showed no bias for a particular segment and were expressed throughout the entire length of the primitive streak (Fig. 5E). Mesoderm-associated populations such as paraxial, presomitic and somitic mesoderm can be seen as transcriptionally defined subsets during later development (Guibentif et al., 2021), but these were not yet observed in this EB stage-derived primitive streak dataset. In summary, transcriptional signatures of cells obtained from anterior to posterior primitive streak regions were notably heterogeneous and already showed a bias towards particular molecular states.

Having identified molecular differences along the anterior to posterior axis of the primitive streak, we next used the extended cell atlas to predict the fates of regions A to D computationally and examined how these compare with experimentally determined cell fates. To this end, using computational fate inference, the closest neighbouring cells of the E7.5 primitive streak cells were identified in the atlas, and downstream fates were inferred based on the W-OT framework outlined above. Based on an initial survey of all predicted fates, a subset of major predictions was selected and visualised in so-called fate plots (Fig. 6C, Fig. S16E), which revealed associations between predicted cell fates and specific portions of the primitive streak. Cells from the posterior-most region A, for example, were primarily associated with mesodermal fates such as allantois, somites and the non-neural ectoderm, whereas the notochord and neural tube fates clearly favoured the most anterior region D (Fig. 6C, Fig. S16E).

To obtain experimental cell fate data, E7.5 embryos were orthotopically grafted with primitive streak regions A to D, cultured for 24 h and the fate contributions of the donor regions analysed at E8.25 by microscopic assessment and scRNA-seq (Fig. 6A). To facilitate analysis of cell fates, transgenic embryos carrying a ubiquitous membrane-bound tdTomato (mTom; Muzumdar et al., 2007) were used as graft donors. Cultured whole embryos and sections were immunostained (Fig. 6B, Fig. S16B) and careful observation of the location of mTom⁺ cells showed clear differences in tissue distribution between the donor primitive streak regions (summarised in Fig. 6D and detailed in Table S8). Notably, region A of the primitive streak contributed predominantly to the most posterior vasculature of the embryo, including the paired dorsal aortae and so-called vessel-of-confluence (VoC; Rodriguez et al., 2017), and to the allantois and some YS cells. Regions B and C also contributed to the paired dorsal aortae, as well as to somites, and to the LPM. Region D contributed mostly to the neural tube and notochord (Fig. 6D). To analyse the cellular contribution of the grafted primitive streak regions at the transcriptional level, we flow-sorted the mTom⁺ single cells from additional batches of cultured embryos and performed scRNA-seq (Smart-seq2). After low-level pre-processing, cells were mapped onto the extended atlas as before. Compared with the cells directly isolated from the E7.5 streak, donor primitive streak-derived cells in the E8.25 cultured embryos mapped to almost twice as many cell types (38 versus 63, respectively; compare Figs 5B and 6E), highlighting their advanced differentiation and diversification.

Comparison between the microscopically and transcriptionally-observed fates of the grafted primitive streak regions showed good concordance, with only subtle differences in the contribution to distinct lineages (Fig. 6D,F; microscopically observed fates are based on embryo counts, whereas transcriptionally observed fates are based on individual cell counts). Importantly, comparison between observed and predicted fates showed largely concordant patterns across computationally inferred fates, molecularly (scRNA-seq) mapped fates and microscopically assigned fates (Fig. 6C,D,F). It is worth noting that scRNA-seq-based analysis of fates has substantially enhanced granularity over microscopic analysis, with fate assignment based on our extended atlas for example providing up to 88 different cell types/molecular states for fine-grained annotation. In summary, combining multi-disciplinary approaches, we established primitive streak cell ‘end fates’ at a single-cell resolution. Placing these in the context of the extended mouse gastrulation and organogenesis atlas allowed establishing the proof-of-concept for a fate-predictive algorithm, able to forecast fate trajectories of individual primitive streak cells.

To realise the full potential of single-cell atlas efforts for developmental biology research, molecular profiling datasets need to (1) provide sufficient sampling density to enable a time-series capable of capturing rapid developmental processes and (2) contain enough single cells sequenced at reasonable depth, well annotated and provided to the broader community through user-friendly web portals. Here, we report such a resource covering the crucial stages of mouse gastrulation and early organogenesis, from E6.5 to E9.5 sampled every 6 h in 13 individual time steps. This dataset transforms our previous effort by more than tripling the number of cells and more than doubling the number of defined cell states. Furthermore, we revisit classical embryo grafting experiments with transgenic and single-cell analysis tools, providing a foundation for future efforts aiming to fully reconstruct cell lineage trees.

Previous single-cell atlas efforts from us and others placed observed molecular states into the context of existing knowledge of mouse development (reviewed by Tam and Ho, 2020) and further integrated by Qiu et al. (2022, 2023 preprint), altogether contributing to a deeper understanding of the molecular heterogeneity accompanying lineage differentiation during mouse embryogenesis. Our present study advances beyond these efforts in several impactful ways. By incorporating W-OT, spatial information and encompassing the developing YS region, our analysis achieves a comprehensive characterisation of blood and endothelial formation spanning the crucial developmental window of E8.5-E9.5. Our dense sampling, deep sequencing and regional sub dissection of embryos allowed us to identify cell states not previously observed in early mouse development, such as the endotome and VACs. Moreover, our trajectory analyses suggest that the endotome-like progenitor population is related to downstream endothelial cells, suggesting a greater intricacy in endothelial differentiation than previously appreciated.

Similar cell types have previously been described in zebrafish and chicken embryos and recently in pluripotent stem cell-derived human axioloids (Nguyen et al., 2014; Pardanaud et al., 1996; Pouget et al., 2006; Sato et al., 2008; Yamanaka et al., 2023), highlighting the value and complementarity of using a variety of different model organisms. Moreover, the zebrafish studies suggested that endotome and VACs may contribute to a cellular niche that promotes intra-embryonic formation of blood stem/progenitor cells (Nguyen et al., 2014). In the mouse, the best-understood site of HSC emergence is the haemogenic endothelium of the dorsal aorta where pro-HSCs are generated from E9.5 (Rybtsov et al., 2014). These haemogenic endothelial cells are lateral plate- or splanchnopleuric mesoderm-derived and were shown to be replaced by somite-derived endothelial cells concomitant with the extinction of HSC generation (Pouget et al., 2006). Interestingly, another report suggests that 7 days later during mouse gestation, endothelial cells in fetal bone marrow undergo haemogenic transdifferentiation and produce BPs and differentiated cells (Yvernogeau et al., 2019). Pax3-Cre lineage tracing further suggested a somitic origin of those haemogenic cells.

Intriguingly, more recent unpublished Tbx6⁺ (Yvernogeau et al., 2020 preprint) and Pax3⁺ lineage tracing studies (Lupu et al., 2022 preprint) have been performed at embryonic time points (E7.5-E11.5) exploring developing anatomical regions (anterior/medial embryonic sections) that overlap with the emergence of endotome, VACs and endotome-derived endothelial cells in our extended atlas. These studies suggest that somitic mesodermal precursors give rise to endothelial cells of the limb and trunk region in developing mouse embryos, as well as stromal cells juxtaposed to the dorsal aorta, which may be like the VAC population we identify in this present study. In these studies, somite-derived endothelial cells make only a minimal contribution to dorsal aorta endothelium by E9.5-10.5. Future work will need to dissect the underlying connection we identify here between a putative endotome cell population, VACs and embryo proper endothelium. Such studies will: (1) enhance our understanding of early blood and endothelium development, (2) likely reveal principles relevant to cell plasticity and potential in other developmental contexts, and (3) provide new mechanistic insights that could be exploited to control differentiation, for example for directed differentiation of pluripotent cells for cell therapy applications.

By combining scRNA-seq with classical embryo grafting experiments, we show how single-cell atlases provide a powerful resource to revisit classical concepts of developmental biology. The results presented in this study are limited in scope but already illustrate how timecourse-guided trajectory reconstruction performed with computational methods such as W-OT represents a promising approach to dissect complex developmental processes. Our findings also support substantial fate bias present in the E7.5 primitive streak, based on deeply sequenced single-cell profiles. As we did not perform heterotopic transplants, our experiments could not assess fate plasticity nor full fate potential. Moreover, deeper and more comprehensive transplant experiments should be designed around new experimental tools not available when this work was performed, including extended in vitro embryo culture (Aguilera-Castrejon et al., 2021), as well as single-cell barcoding to permit reconstruction of single-cell phylogenies (Bowling et al., 2020).

A future thereby emerges where a confluence of complementary technologies will transform our understanding of early mouse development to a level previously only attained with non-mammalian organisms. Building on data resources such as the one reported here, mechanistic insights will continue to require perturbation experiments, with the important proviso that carefully chosen experimental perturbations have the added benefit to provide new insights into disease processes, in particular congenital defects associated with mutations in developmental regulator genes. Our extended mouse gastrulation atlas provides the developmental biology community with a significant new resource to probe novel hypotheses concerning cell fate acquisition and lineage commitment during embryogenesis.

All procedures were performed in strict accordance to the UK Home Office regulations for animal research. Mouse embryos were collected under the project licence number PPL 70/8406. Animals used in this study were 6-10 week-old females, maintained on a lighting regime of 14 h light and 10 h dark with food and water supplied ad libitum. Following wild-type C57BL/6 matings, females were killed by cervical dislocation at E8.5, E8.75, E9.0, E9.25 and E9.5. The uteri were collected into phosphate buffered saline (PBS) with 2% heat-inactivated fetal calf serum (FCS) on ice and the embryos were immediately dissected and processed for scRNA-seq. For each time point, four embryos were selected based on morphology and somite counts, in order to span the range expected for the given time point according to Davidson (1989) and processed individually. The exception is the E9.5 time point, where embryos were smaller than expected and only two embryos were collected with lower somite numbers. The YS was systematically separated from the rest of the embryo and processed as a separate sample. Of the four selected embryos at each time point, two were partitioned in defined anterior-posterior sections dissociated as separate samples and two were dissociated as bulk and the suspension then divided into two separate samples for 10x RNA-seq analysis. For the E8.5 partitioned embryos, they were divided into two halves, with the cut being made at the fourth somite level (i.e. right before the fourth somite level). The anterior portion (including headfolds, branchial arches and heart rudiment) and the posterior portion (including allantois, hindgut, primitive streak) were processed as individual samples. For E8.5 bulk-dissociated embryos, the single-cell suspension was divided into two separate samples for 10x RNA-seq analysis.

E8.75-E9.5 embryos were divided into three segments, with cuts made below the otic pit and below the heart (at the 10th-12th somite level). The anterior-most portion (including brain structures anterior to rhombomere 6 and branchial bars), mid-portion (including the heart and remains of vitelline vessels) and posterior portion (including allantoic structures, hindgut and posterior-most somites) were then singularised and further processed as separate samples. Single-cell suspensions were prepared by incubating the samples with TrypLE Express dissociation reagent (Life Technologies) at 37°C for 7 min under agitation and quenching in PBS with 10% heat-inactivated serum. The resulting single-cell suspension was washed and resuspended in PBS with 0.4% bovine serum albumin and filtered through a Flowmi Tip Strainer with 40 μm porosity (Thermo Fisher Scientific, 136800040). Cell counts were then assessed with a haemocytometer. Single-cell RNA-seq libraries were generated using the 10x Genomics Chromium system (version 3 chemistry), and samples were sequenced according to the manufacturer's instructions on the Illumina NovaSeq 6000 platform.

The mouse gastrulation atlas was processed exactly as described in Pijuan-Sala et al. (2019). This timecourse experiment contains 116,312 cells distributed across nine time points sampled across E6.5-E8.5 at 6 h intervals. Only two samples of this dataset (where a sample is a single lane of a 10x Chromium chip) contained pooled embryos staged across several time points. Cells from these samples are denoted as ‘mixed gastrulation’ in the metadata.

Raw reads were processed with Cell Ranger 3.1.0 using the mouse reference 1.2.0, mm10 (Ensembl 92) and default mapping arguments. The following steps of pre-processing were performed with R using the same functions, parameters and software versions broadly described in Pijuan-Sala et al. (2019): swapped molecule removal, cell calling, quality control, normalisation, selection of highly variable genes, doublet removal, batch correction and stripped nucleus removal. Therefore, singularity containers available in https://github.com/MarioniLab/EmbryoTimecourse2018 provide the necessary software for reproducing these steps.

Log transformed normalised counts obtained from Pijuan-Sala et al. (2019) and those generated here were integrated into an extended timecourse experiment across mouse developmental stages E6.5-E9.5 (13 time points). Highly variable genes (HVGs) were calculated using ‘trendVar’ and ‘decomposeVar’ from the scran R package, with loess span of 0.05. Genes that had significantly higher variance than the fitted trend (Benjamini–Hochberg-corrected P<0.05) were retained. Genes with mean log2 normalised count <10^–3, genes on the Y chromosome, the gene Xist and the reads mapping to the tdTomato construct (where applicable) were excluded. Hence, this expression matrix contained 23,972 genes and 430,339 cells. Batch correction with fastMNN function from scran (Lun et al., 2016) was performed as described above, resulting in 5665 genes and 75 batch-corrected principal components.

Metadata annotations such as embryonic developmental stages and cell types were assigned to the mixed time points (annotated as mixed gastrulation in the metadata) and newly generated E8.5 samples, respectively, using a strategy based on fastMNN. In this approach, UMI counts from both the reference and the query datasets are merged, normalised and log transformed together. Then, HVGs and top principal components are computed to subsequently use fastMNN for re-scaling the principal component analysis (PCA) space from both datasets. The annotations from the reference metadata are assigned to the query data as the mode among k nearest neighbours (KNN) between the query and the reference PCA subspaces using queryX function from Biocneighbours. The number of nearest neighbours is chosen depending on the resolution of transferred annotations. Mixed gastrulation time points were allocated to embryonic developmental stages using the 30 nearest neighbours queried from the top 50 batch-corrected principal components from the E6.5-E8.5 reference dataset. New E8.5 cell type annotations were assigned using ten nearest neighbours with respect to E8.5 cells from the reference atlas in the corresponding E8.5 subspace of the integrated batch corrected PCA described above.

The W-OT approach was conceived to model timecourse experiments of developmental processes as a generalisation of a stochastic process using unbalanced optimal transport theory. Thus, it allows estimating the coupling probabilities between cells of consecutive time points, while taking into account cell growth and death (Schiebinger et al., 2019). The transport maps of consecutive time points were constructed over the entire set of cells with W-OT (1.0.8.post1) using default settings, except for skipping the dimension-reduction step, and instead using the batch-corrected principal components as input as well as three iterations for learning the cell growth rate. Embryonic stages E9.25 and E9.5 were collapsed into a single time point for computing the transport maps because the somite count (a more accurate measure of developmental stage) overlapped between these embryos.

The transport maps explained above were used to estimate the full trajectories of every cell population present at E8.5. That is, for each cell population (i.e. Erythroid3) the coupling probabilities were used to reconstruct the sequences of ancestor and descendant distributions by pushing the cell set through the transport matrix backwards and forwards, respectively. Cells were allocated as descendants from a cell population at E8.5 by selecting those with maximum mass across all trajectories at time points E8.75-E9.5.

Mapping publicly available data from the atlas of brain development (La Manno et al., 2021) and both atlases of gut development (Nowotschin et al., 2019) against the extended gastrulation atlas was performed following the strategy based on fastMNN mentioned above. However, in this case the matrices of principal components from the query dataset were randomly split into subsets smaller than 10,000 cells for fastMNN and merged back to generate a single mapping output.

A combination of complementary strategies was used for defining final cell type annotations (Fig. S1). First, cell type annotations from Pijuan-Sala et al. (2019) were transferred within overlapping time points (E8.5) and cell descendants were estimated as described above (Fig. S1A,B). The landscape was then split into two subsets, a mesodermal and an ectodermal-endodermal landscape (the latter including neuromesodermal progenitors) (Fig. S1C), and subsequently clustered using the scanpy implementation of the Leiden algorithm (Traag et al., 2019) with a resolution of 5 (Fig. S1D). Highly variable genes and batch corrected principal components were recomputed on the subsets before clustering. Then, annotations from the original atlas expanded through mapping and estimation of cell descendants were manually refined by means of differential expression analysis between clusters (using findMarkers from the scran package version (Lun et al., 2016), as well as literature and visual inspection of gene markers resulting in major 88 cell type populations (Fig. 2). In addition, we used the FindAllMarkers (min_pct=0.25, logfc_threshold=0.25) function in Seurat (4.2.0) (Hao et al., 2021) to identify marker genes for all 88 cell types (Table S1).

The strategy for generating the haemato-endothelial landscapes was based on the so-called W-OT fate matrix. In brief, a W-OT fate matrix is a transition probability matrix (the rows of this matrix add up to 1) from all cells to a number of target cell sets at a given time point, commonly the ending point of the timecourse experiment. To fully cover the haemato-endothelium, three W-OT fates matrices were computed, one for each landscape presented in our results. One using E9.25-E9.5 YS blood and YS endothelial cells as targets, another one with YS blood as well as both YS and embryonic endothelial tissues, and a third one where only endothelial cell types were considered targets (Figs 3, 4 and Fig. S10, respectively). Then, to identify cells that form a differentiation trajectory towards any ofthe abovementioned cell fates, for every cell in the landscape the likelihood probabilities associated with these fates and the remaining cells at E9.25-E9.5 cell types (i.e. cells grouped as ‘other fate’) were compared using the log odds or ratio of probabilities. The log odds is defined as the logarithm of the ratio of probabilities of two different categorical and mutually exclusive outcomes log(p/1−p). These ratios provided a quantitative setting for more interpretable thresholds when selecting cells that potentially belong to a cell fate trajectory. These cells are then retained and used to generate new layouts including only potentially fated cells towards blood and endothelium: that is, selecting pluripotent and mesodermal cells above a reasonable log odds threshold, recomputing HVGs, correcting for batch effects and generating new force-directed graphs from inferred trajectories. Cells with log odds >0 were considered fate biased towards the set of selected population targets, such as in the proof of concept performed in Fig. 3. As the log odds is computed by summing over the cell probabilities for all population targets (here denoted as p), divided by the probability of all other fates grouped together (1-p); cells with negative values but close to 0, although with higher uncertainty, might well be contributing to the population targets of interest. Thus, for more complex landscapes the threshold was lowered to log odds >−1 (Fig. 4, Fig. S10) in order to exclude all cells not associated with haemato-endothelial fates with high confidence, while keeping a higher degree of uncertainty for those retained.

To generate the UMAP layout of the whole embryo, we used the top 50 batch corrected principal components to generate a BBKNN graph (Polanski et al., 2020) and then the scanpy implementation of UMAP, with parameter min_distance=0.99. To generate the force-directed layouts of the haemato-endothelial landscapes, we recomputed HVGs for each subset as well as batch correction of PCA manifolds (we again retained the top 50 principal components). We then built a KNN graph (K=50) and used ForceAtlas2 (Jacomy et al., 2014) implementation of force-directed layouts included in scanpy.

Slingshot (2.7.0) (Street et al., 2018) was used, in a semi-supervised fashion, to infer differentiation trajectories and produce average pseudotimes along various putative lineage trajectories that were highlighted during W-OT when generating the YS landscape (Fig. 3 and Fig. S7; primitive, YS definitive and YS endothelial trajectories) and the endothelial landscape (Figs S10 and S11; anterior/medial and posterior trajectories) using the getLineages and getCurves functions (default settings, starting clusters were provided). The first 45 batch corrected principal components were provided to getLineages as input data for the respective landscapes as well as cluster labels generated using the default Seurat (4.2.0; Hao et al., 2021) functions FindNeighbors (dims=1:50, reduction=‘PCA’) and FindClusters (resolution=3) for the respective landscapes. Cells that were members of relatively small clusters (<100 cells) were excluded by setting their cluster value equal to −1 when running getLineages. Before running getLineages, supervised filtering of cells based on metadata (stage, anatomy and predicted anatomy labels) was performed. Predicted anatomy labels for embryo proper cells in the haemato-endothelial landscape were determined by performing label transfer from neighbouring cells with specific anatomical labels (anterior, posterior, anterior section, medial section and posterior section) to cells from embryo proper sections using the embeddingKNN function from StabMAP (Ghazanfar et al., 2023) (type=’uniform_fixed’, k_values=5, cords=top 50 PCs). For the YS primitive blood trajectory inference, cells from stages E6.5-E8.25 were considered. For the YS definitive blood trajectory, cells from stages E7.75-E9.5 were considered. Finally, for the YS endothelial trajectory inference cells from the YS and pooled regions were considered and blood cells (EMP, MEP, BPs, erythroid, megakaryocyte progenitors) were excluded. For the anterior/medial endothelial trajectory inference, cells from stages E8.5-E9.5 were considered and cells from anatomical posterior section and posterior were excluded, as well as embryo proper cells with predicted anatomy labels of posterior section and posterior. For the posterior endothelial trajectory inference, cells from stages E7.5-E9.5 were considered and cells from anatomical anterior section, medial section and anterior were excluded, as well as embryo proper cells with the same predicted anatomy labels (anterior section, medial section and anterior). When multiple, similar lineages towards downstream blood or endothelial cells were predicted by getLineages, average pseudotimes were calculated using the slingAvgPseudotime function. These average pseudotimes were used in subsequent Tradseq analyses.

To identify genes with altered expression over the inferred trajectories the average pseudotimes were provided to Tradeseq (1.10.0; Van den Berge et al., 2020). Generalised additive models were fit (fitGAM, nknots=6, cellWeights=1) to genes that were highly variable [top 2000 highly variable genes were identified using VariableFeatures function in Seurat (4.2.0)] among the cells that were part of the inferred trajectories. The associationTest function was then used to identify which of these genes had altered expression over the inferred trajectories. ComplexHeatmap (2.15.3; Gu et al., 2016) was used to visualise the expression patterns of the genes with the highest 300 waldStat scores (P-value<0.01 and meanLogFC>2). Metascape (Zhou et al., 2019) was used to identify Gene Ontology (GO) terms that were enriched for clusters of genes that were associated with the anterior/medial and posterior endothelial differentiation trajectories (Tables S4 and S5).

To identify predicted ligand-receptor interactions among cells that were present in the YS landscape, the haemato-endothelial landscape and the anterior/medial sections of the cells in the haemato-endothelial landscape (Table S7), we used the computeCommunProb (default settings) and filterCommunication (min_cells=10) in CellChat (1.6.1; Jin et al., 2021).

Subsequent to the construction of a complete haemato-endothelial landscape, Louvain (Subelj and Bajec, 2011) clusters (as implemented in the R package igraph; https://igraph.org) were generated from the batch-corrected principal components obtained as explained above (Fig. 4C). Cluster 20 was of particular interest for investigating both shared and non-shared gene expression profiles between the anatomically distinct populations YS EMP and non-yolk sac (NYS) EMP, and a subset of endotome cells that is placed next to them (Fig. 4D,E). Differential gene expression analysis was performed using the function findMarkers from the scran package (Lun et al., 2016). The statistical test performed by findMarkers uses a general linear model and moderated t-statistics to perform differential expression, as implemented in the R package limma (Ritchie et al., 2015). The threshold established was FDR<0.05 and LogFoldChange>0.5. Furthermore, only genes among the top 20 rank genes are displayed in heatmaps (considering that some genes can have tied ranks, gene lists usually contain more than 20 genes). CCA as implemented in the R package CCA (Gittins, 1985; Kanti Mardia and Bibby, 1979) was applied to mean metacell expression values of YS EMPs and the endotome-derived VACs. Positively correlated genes were then identified by selecting only those with positive coefficients (canonical scores) over the two canonical variates and visualised using a scatter plot fitting a linear regression. Importantly, such analysis was performed using Metacells (Baran et al., 2019) to strengthen gene expression signals.

We identified metacells (i.e. groups of cells that represent singular cell states from single-cell data) with the goal of achieving a resolution that retains the continuous nature of differentiation trajectories while overcoming the sparsity issues of single-cell data. We adopted the SEACells implementation (Persad et al., 2023). Following the method guidelines, metacells were computed separately for each sample using approximately one metacell for every seventy-five cells. Following metacell identification, we regenerated the gene expression matrices summarised at the metacell level. Sample-specific count matrices were then concatenated and normalised together. Metacells were used exclusively to inspect the correlation between genes expressed in YS EMPs and endotome-derived VACs, as well as the lack of a clear haemogenic signature. However, metacells were computed for all the haemato-endothelial landscapes.

Sequencing reads were aligned against the mm10 genome following the GRCm38.95 genome annotation (Ensembl 95) using STAR version 2.5 (Dobin et al., 2013) with --intronMotif and 2-pass option mapping for each sample separately (i.e. --twopassMode basic) to improve detection of spliced reads mapping to novel junctions. Sequence alignment map (SAM) files were then processed with samtools-1.6 in order to generate the corresponding count matrix with htseq-count from HTSeq, version 0.12.4. Cells with total counts lower than 10,000 reads were filtered out from downstream analysis based on observed library size distributions. The mitochondrial fraction of reads per cell and library complexity were computed as part of the quality control (QC). Only two cells showed a mitochondrial fraction larger than 4% and the lowest number of genes detected was 3113. Thus, no cells were excluded during QC. Size factors were computed with scran and log-transformed normalised counts were obtained using the logNormCounts function from scater using size factors centred at unity before calculation of normalised expression values. Then, HVGs were extracted using ModelVar function from scran, according to a significant deviation above the mean-variance fitted trend (BH-corrected P<0.05). These genes were selected for computing a PCA. The top 50 PCs were retained for batch correction using reducedMNN, which is an updated version of fastMNN. The resulting batch-corrected PCs were used for clustering and visualisation in downstream analysis.

Following the strategy to transfer annotations from the atlas to other datasets previously described above, cell type annotations were assigned to primitive streak cells according to the 15 nearest neighbours between the two manifolds (query data and atlas reference). The annotation transfer was performed separately for each batch to avoid introducing an extra confounding factor when applying MNN batch correction. Also, before mapping against the atlas, genes with total counts lower than 100 were filtered out from the Smart-seq2 datasets due to the large differences in dropouts. The resulting annotations were projected onto the UMAP extended atlas embryo landscape by highlighting the closest cells on it. By retaining the closest cell in the atlas and selecting the cell fate with highest probability allocated in a W-OT fates matrix, we aimed to not only map the observed transcriptome at E7.5 but make predictions about likely cell fates at future time points (e.g. E8.5 and E9.5).

Mouse embryos were collected under the UK Home Office project licence number PP9552402. Animals used in this study were 6-16 week-old wild-type females of a mixed (CBAxC57BL/6) F2 background, maintained on a lighting regime of 14 h light and 10 h dark with food and water supplied ad libitum. Following timed mating crosses between transgenic mTmG [B6.129(Cg)-Gt(ROSA)26Sor^{tm4(ACTB−tdTomato,-EGFP)Luo}/J (Muzumdar et al., 2007)] males and wild-type females, females were killed 7 days after a vaginal plug was observed via a Schedule 1 method. The uterine tissue and the decidua were cut away to extract the E7 embryo. The Reichert's membrane was peeled off and embryos were staged according to Downs and Davies (1993). EB-stage embryos were taken for further analysis. For primitive streak dissections, the extra-embryonic part of the conceptus was cut off, the anterior and posterior embryonic portions were separated, and the posterior portion flattened out, taking care not to mix up the orientation. The posterior portion containing the primitive streak was cut at the halfway point and each half was halved again creating four similar-sized primitive streak segments. All embryo and primitive streak dissections were performed in Dulbecco's PBS (with CaCl₂ and MgCl₂, Gibco) supplemented with 10% FCS (Gibco), 50 U/ml penicillin and 50 U/ml streptomycin.

For scRNA-seq of primitive streak cells, individual primitive streak regions were collected in separate Eppendorf tubes in FACS buffer [PBS without CaCl₂ and MgCl₂ (Gibco) supplemented with 10% FCS (Gibco), 50 U/ml penicillin and 50 U/ml streptomycin]. After centrifugation at 200 g for 5 min, tissues were resuspended in 200-250 μl of TrypLE Express (Life Technologies) and incubated at 37°C for 7 min with regular agitation and dissociated by pipetting. Cells were washed with 1 ml of FACS buffer, centrifuged at 200 g for 5 min and resuspended in FACS buffer with Hoechst 33358 for sorting. In order to reduce cell loss from these small cell populations, wild-type adult mouse thymocytes were added to the samples and were gated out based on their typical FSC-SSC profiles and absence of tdTomato expression. Each primitive streak region yielded on average 200 cells. Cells were kept on ice throughout the procedure. Individual cells were sorted directly into 96-well plates (Starlab, E1403-1200) into 2.3 µl of lysis buffer containing SUPERase-In RNase Inhibitor 20 U/µl (Ambion, AM2694), 10% Triton X-100 (Sigma-Aldrich, 93443) and RNAse-free water. Plates were processed using a combination of Smart-seq2 (Picelli et al., 2014) and mcSCRB-seq (Bagnoli et al., 2018). For detailed protocol see Sturgess et al. (2021). Pooled libraries were run on the Illumina HiSeq4000 at Cancer Research UK Cambridge Institute Genomics Core.

mTmG transgenic embryos were dissected to harvest four equal primitive streak regions as described above. Each segment was deposited into a 50 μl dissection buffer drop on the inside lid of a 5 cm dish and labelled. One larger drop was deposited in the centre, into which up to five dissected wild-type embryos of the same stage (EB) were transferred. Primitive streak regions were grafted into orthotopic positions in the wild-type embryos using a pulled glass capillary needle attached to a mouth pipette. Cells from each donor region were divided over 2-3 recipient embryos with each recipient embryo receiving 50-100 primitive streak cells, mostly as a chain of cells rather than in suspension. Embryos were gently washed in a drop of culture medium and gently transferred into individual wells of a 96-well Ultra-Low adhesive plate (Corning) into 200 μl culture medium. Culture medium was 100% rat serum (Envigo, custom collected). Embryos were cultured at 37°C, 5% CO₂ for 23-26 h.

Embryos were taken out of the incubator and assessed for (the beginnings of) a heartbeat. Only embryos that lacked deformities and had a heartbeat were used for further analysis by microscopy or scRNA-seq to assess the fates of the grafted primitive streak cells. For microscopically observed fates, whole embryos/concepti were immunostained in U-bottom 96-well plates as follows. Whole embryos were fixed in 4% paraformaldehyde in PBS for 20-30 min at room temperature (RT) and washed three times in staining buffer (0.4% Triton X-100, 2% FCS in PBS). Embryos were incubated in individual Eppendorf tubes in staining buffer for 3-6 h at RT on a shaking rack. Reporter gene expression (mTom) was enhanced with anti-RFP antibody (Rockland, 600-401-279, 1:400) and embryonic vasculature visualised with anti-CD31 antibody (R&D Systems, AF3628, 1:200).  Primary antibodies were diluted in the staining buffer and 50-100 μl used for each well/embryo. Incubation lasted overnight at 4°C, gently shaking. The next day, embryos were washed in staining buffer for up to 5-8 h at RT on a shaking rack. Whole mounts were then incubated with the secondary antibody solution (Alexa Fluor 555, Invitrogen, A31572, 1:200; Alexa Fluor 647, Invitrogen, A21447, 1:200 ) with nuclear dye (Hoechst 33358) overnight at 4°C and subsequently washed again for several hours or up to 24 h with several changes of the buffer. Wholemount embryos were imaged to assess overall contribution of mTom in relation to embryo anatomy (Zeiss 780 inverted confocal microscope at 10×). Embryos were next embedded in OCT and snap-frozen on ethanol-dry ice mix followed by cryosectioning, retaining every section to capture all mTom-expressing cells. Sections were then assessed, and notes taken of all areas with mTom⁺ contribution and careful analysis and scoring for detailed contribution to a selection of most anatomically distinct cell lineages (Table S8). Immunofluorescence in sections was imaged using a Zeiss 780 inverted confocal microscope with 25× (LD LCI PA 25×/0.8 DIC), 40× (LD C-Apochromat 40×/1.1) and 63× (Plan Apochromat 63×/1.4) objectives with oil immersion. Images were stored and analysed in OMERO. Adobe Photoshop was used for final image processing and Adobe Illustrator to generate figures.

We thank all the people and the different funding contributing to this project. We thank Sarah Kinston for technical assistance. We thank Biomedical Services staff at Oxford for animal care. We acknowledge Kevin Clark in the flow cytometry facility at the MRC WIMM [supported by the MRC HIU; MRC MHU (MC_UU_12009); NIHR Oxford BRC; Kay Kendall Leukaemia Fund (KKL1057), John Fell Fund (131/030 and 101/517), the EPA fund (CF182 and CF170) and by the MRC WIMM Strategic Alliance awards G0902418 and MC_UU_12025] for providing cell sorting services. We thank Christoffer Lagerholm and Jana Koth of the Wolfson Imaging Centre Oxford [supported by the MRC via the WIMM Strategic Alliance (G0902418), the MHU (MC_UU_12009), the HIU (MC_UU_12010), the Wolfson Foundation (18272) and the Wellcome Trust (Micron 107457/Z/15Z)] for their assistance.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201867.reviewer-comments.pdf