Single-cell, whole-embryo phenotyping of mammalian developmental disorders

Abstract

Mouse models are a critical tool for studying human diseases, particularly developmental disorders¹. However, conventional approaches for phenotyping may fail to detect subtle defects throughout the developing mouse². Here we set out to establish single-cell RNA sequencing of the whole embryo as a scalable platform for the systematic phenotyping of mouse genetic models. We applied combinatorial indexing-based single-cell RNA sequencing³ to profile 101 embryos of 22 mutant and 4 wild-type genotypes at embryonic day 13.5, altogether profiling more than 1.6 million nuclei. The 22 mutants represent a range of anticipated phenotypic severities, from established multisystem disorders to deletions of individual regulatory regions^4,5. We developed and applied several analytical frameworks for detecting differences in composition and/or gene expression across 52 cell types or trajectories. Some mutants exhibit changes in dozens of trajectories whereas others exhibit changes in only a few cell types. We also identify differences between widely used wild-type strains, compare phenotyping of gain- versus loss-of-function mutants and characterize deletions of topological associating domain boundaries. Notably, some changes are shared among mutants, suggesting that developmental pleiotropy might be ‘decomposable’ through further scaling of this approach. Overall, our findings show how single-cell profiling of whole embryos can enable the systematic molecular and cellular phenotypic characterization of mouse mutants with unprecedented breadth and resolution.

Main

For more than 100 years, the laboratory mouse (Mus musculus) has served as the quintessential animal model for studying human diseases¹. For developmental disorders in particular, mice have been transformative, as a mammalian system that is nearly ideal for genetic analysis and in which the embryo is readily accessible⁶.

At its inception, mouse genetics relied on spontaneous or induced mutations resulting in visible physical defects that could then be mapped. However, gene-targeting techniques later paved the way for ‘reverse genetics’ (that is, analysing the phenotypic effects of intentionally engineered mutations)⁶. Through systematic efforts such as the International Knockout Mouse Consortium, knockout (KO) models are now available for thousands of genes⁷. Furthermore, with genome editing^8,9, it is increasingly practical to delete individual regulatory elements¹⁰.

Phenotyping has also grown more sophisticated. Conventional investigations of developmental syndromes typically focus on one organ at a specific stage (for example, combining expression analyses, histology and imaging to investigate a visible malformation)¹. The Mouse Clinic, involving a battery of standardized tests, reflects a more systematic approach¹¹, but phenotypes detected through such tests (for example, behavioural and electrophysiological) may require years of additional work to link them to molecular and cellular correlates. Furthermore, it is often the case that an intentionally engineered mutation results in no detectable abnormality¹². In such instances, it remains unknown whether there is truly no phenotype, or whether the methods used are simply insufficiently sensitive. In sum, phenotyping has become ‘rate limiting’ in mouse genetics.

Single-cell molecular profiling offers a potential path to overcome such barriers. As a first step, we and others have applied single-cell RNA sequencing (scRNA-seq) to profile wild-type mouse development at the scale of the whole embryo^{3,13,14,15,16,17,18}. Applying scRNA-seq to mouse mutants, several groups have successfully unravelled how specific mutations affect transcriptional networks and lead to altered cell fate decisions in individual organs¹⁹. However, there is still no clear framework for analysing such data at the whole-embryo scale.

scRNA-seq of 101 mouse embryos

We set out to establish whole-embryo scRNA-seq as a scalable framework for the systematic molecular and cellular phenotyping of mouse genetic models. We collected 103 mouse embryos, including 22 different mutants and 4 wild-type strains at embryonic day (E)13.5, and generally 4 replicates per strain (Fig. 1a). Mutants were chosen to represent a spectrum of phenotypic severity ranging from established pleiotropic disorders to KOs of individual regulatory elements.

Fig. 1: Single-cell transcriptional profiling of 103 whole mouse embryos staged at E13.5.

a, Categories of mutants (left) analysed by whole-embryo profiling with sci-RNA-seq3 (right). WT, wild type. KI, knock-in; DKO, double KO. −/Y, hemizygous b, Embeddings of pseudobulk RNA-seq profiles of MMCA embryos in PCA space with visualization of the top three PCs. Datapoints are coloured by background strain of the embryo. The black dashed ovals highlight two major groups corresponding to FVB versus other backgrounds. Embryo 104 is a clear outlier. c, Embeddings of pseudobulk RNA-seq profiles of MOCA³ and MMCA embryos in PCA space defined solely by MOCA, with MMCA embryos (grey) projected onto it. The top two PCs are visualized. Coloured points correspond to MOCA embryos of different stages (E9.5–E13.5), and grey points to MMCA embryos (E13.5). The dashed line highlights five MMCA embryos that are co-localized with E11.5 or E12.5 embryos from MOCA. Three are from Scn11a-GOF (embryos 33, 34 and 36), Carm1-KO (embryo 101) and wild-type (embryo 41) C57BL/6 strains. d, Number of cells profiled per embryo for each strain. Centre lines show medians; box limits indicate 25th and 75th percentiles; replicates (n = 3 for wild-type C57BL/6, n = 4 for all others) are represented by dots. Genotypes are listed by median cell number in ascending order. e, Three-dimensional UMAP visualization of wild-type subset of MMCA dataset (215,575 cells from 15 embryos). Cells are coloured by major-trajectory annotation. PNS, peripheral nervous system.

We grouped mutants into four categories (Fig. 1a and Supplementary Table 1). The first category, pleiotropic mutants, consisted of embryos with KOs of developmental genes expressed in several organs (Ttc21b KO, Carm1 KO and Gli2 KO), and two mutations of the Sox9 regulatory landscape suspected to have pleiotropic effects (Sox9 topological associating domain (TAD) boundary knock-in; Sox9 regulatory inversion (INV))^5,20,21,22. The second category, developmental disorder mutants, consisted of embryos intended to model specific human diseases (Scn11a gain of function (GOF), Ror2 knock-in, Gorab KO and Cdkl5 −/Y (hemizygous))^23,24,25. The third category consisted of embryos with mutations of loci associated with human disease (Scn10a/Scn11a double KO, Atp6v0a2 KO, Atp6v0a2^R755Q and Fat1 TAD KO)^26,27. The fourth category consisted of embryos with prospective deletions of cis-regulatory elements, including of TAD boundaries near developmental transcription factors (Smad3, Tbx5, Neurog2, Sim1, Smad7, Dmrt1, Tbx3 and Twist1)⁴. As a positive control, this fourth category includes a ZRS distal enhancer (zone of polarizing activity regulatory sequence) KO mutant, which specifically fails to develop distal limb structures²⁸. Except for Scn11a GOF, all mutants were homozygous.

To validate staging, we leveraged our previous mouse organogenesis cell atlas (MOCA), which spans E9.5 to E13.5 (ref. ³). After doublet filtering, we profiled 1,671,245 nuclei (16,226 ± 9,289 per embryo; 64,279 ± 18,530 per strain; median unique molecular identifier count: 843 per cell; median genes detected: 534; 75% duplication rate). Below we refer to this dataset as the mouse mutant cell atlas (MMCA).

Applying principal component analysis (PCA) to ‘pseudobulk’ profiles of the embryos resulted in two groups corresponding to genetic background (Fig. 1b), with FVB embryos clustering separately from C57BL/6J, G4 and BALB/c embryos. However, embryos corresponding to individual mutants did not cluster separately, suggesting that none was affected with severe, global aberrations. A single outlier (embryo 104) was aberrant with respect to cell recovery (n = 1,047) and appearance (Extended Data Fig. 1a).

To validate staging, we leveraged our previous mouse organogenesis cell atlas (MOCA), which spans E9.5 to E13.5 (ref. ³). PCA of pseudobulk profiles of 61 wild-type embryos from MOCA resulted in a first component (principal component 1 (PC1)) strongly correlated with developmental age (Fig. 1c). Projecting pseudobulk profiles of the 103 MMCA embryos to this embedding resulted in most MMCA embryos clustering with E13.5 MOCA embryos along PC1, consistent with accurate staging. However, five MMCA embryos seemed closer to E11.5 or E12.5 MOCA embryos. Four of these were retained as their delay might be explained by their mutant genotype, whereas one wild-type embryo (C57BL/6; embryo 41) was designated a second outlier. We removed cells from the two outlier embryos (embryos 104 and 41) as well as cells with high proportions of reads mapping to the mitochondrial genome (>10%) or ribosomal genes (>5%). This left 1,627,857 cells, derived from 101 embryos (Fig. 1d).

To facilitate data integration, we projected cells from all genotypes to a wild-type-derived ‘reference embedding’ (Methods and Extended Data Fig. 1b,c). Altogether, we identified 13 major trajectories, 8 of which were further stratified into 59 sub-trajectories (Fig. 1e, Extended Data Fig. 2a and Supplementary Table 2), generally covering the expected cell trajectories at this stage of development. These were also generally consistent with our annotations of MOCA, albeit with some corrections as described elsewhere^17,29. Greater granularity for some cell types is probably a consequence of the deeper sampling of E13.5 cells in these new data (Extended Data Fig. 2b).

Mutant-specific variation in cell-type composition

In analysing these data, we pursued three approaches: quantification of gross differences in cell-type composition (this section); investigation of more subtle differences in the distribution of cell states within annotated trajectories and sub-trajectories; and analysis of the extent to which phenotypic features are shared between mutants.

To systematically assess cell-type compositional differences, we first examined the proportions of cells assigned to each of 13 major trajectories. These proportions were mostly consistent across genotypes (Extended Data Fig. 3a), but some mutants exhibited substantial differences. For example, compared to wild-type C57BL/6, the proportion of cells in the neural tube trajectory decreased from 37.3% to 33.7% and 32.6% in the Gli2-KO and Ttc21b-KO strains, respectively, and the proportion of cells in the mesenchymal trajectory decreased from 44.1% to 37.1% in the Gorab-KO strain. These changes are broadly consistent with the gross phenotypes associated with these mutations^20,25,30, but are subject to the caveat of substantial interindividual heterogeneity (Extended Data Fig. 3b). We also observe differences in major-trajectory composition between the four wild-type strains. For example, wild-type FVB and G4 mice consistently had fewer mesenchymal and more neural tube cells than wild-type BALB/c and C57BL/6 embryos (Extended Data Fig. 3c). We further checked for technical effects (for example, experimental batch) that might confound cell-type recovery rates (Extended Data Fig. 4a–c).

We next sought to investigate compositional differences at the level of sub-trajectories. For each combination of background and sub-trajectory, we carried out regression to identify mutations that were nominally predictive of the proportion of cells falling in that sub-trajectory (uncorrected P value < 0.05; beta-binomial regression; Methods). Across 22 mutants, this analysis highlighted 300 nominally significant changes (Fig. 2a and Supplementary Table 3). Owing to the limited number of replicate embryos per strain, our power to definitively call such changes is limited, particularly in the smaller trajectories (Methods and Extended Data Fig. 4d). Nevertheless, two patterns are noteworthy, as follows.

Fig. 2: Cell composition changes for individual mutants across developmental trajectories.

a, The heatmap shows log₂-transformed ratios of cell proportions between each mutant type (y axis) versus the pooled reference (consisting of wild type and other mutants from the same strain; cell counts from replicates were merged) across individual sub-trajectories (x axis). Only nominally significant results are shown (Methods). The pie colour and direction correspond to whether the log₂-transformed ratio is >0 (blue, clockwise) or <0 (red, anticlockwise); the pie size and colour intensity correspond to the scale of the log₂-transformed ratio. log₂-transformed ratios that were >2 or <−2 were manually set to 2 or −2. The number of cells assigned to each sub-trajectory and selected thresholds of proportions (red vertical lines) are shown above. b, Three-dimensional UMAP visualization of the neural tube trajectory, highlighting cells from either Ttc21b-KO (left) or wild-type (middle) C57BL/6, or other C57BL/6 mutants (right), after downsampling a uniform number of cells per plot. c, Haematoxylin and eosin staining of the developing eye of homozygous mutant (Hom), heterozygous mutant (Het) and wild-type Ttc21b E13.5 embryos (Methods). Structures are lost exclusively in homozygous mutants. Scale bars, 500 μm. d, UMAP visualization of co-embedded cells of the limb mesenchyme trajectory from ZRS limb enhancer KO and wild-type FVB embryos. The same UMAP is shown twice, highlighting FVB wild-type (left) or ZRS limb enhancer KO (right) cells. A subpopulation exhibiting extreme loss in the ZRS limb enhancer KO is circled. e, The same UMAP as in d, coloured by expression of marker genes that seem specific to proximal (Meis1 and Meis2) or distal (Hoxa13, Hoxd13, Lhx9 and Msx1) limb development (Supplementary Table 4). MHB, midbrain–hindbrain boundary; Di, diencephalon.

First, Atp6v0a2 KO and Atp6v0a2^R755Q, distinct mutants of the same gene²⁶, exhibit highly consistent patterns of change, both with respect to which sub-trajectories are nominally significant as well as the direction and magnitude of changes. The consistency supports the validity of this analytical approach.

Second, mutants varied considerably with respect to the number of sub-trajectories nominally significant for compositional differences. At the higher extreme, 30 of 54 sub-trajectories were nominally altered by the Sox9 regulatory INV mutation, consistent with the wide-ranging roles of SOX9 in development^31,32. At the lower extreme, TAD boundary KOs exhibited very few changes, consistent with the paucity of gross phenotypes in such mutants¹². Nonetheless, all TAD boundary KOs did show some nominal changes, including specific ones (for example, the lung epithelial and liver hepatocyte trajectories were specifically decreased in Dmrt1 and Tbx3 TAD boundary KOs, respectively).

There were a few extreme examples (for example, in which a sub-trajectory seemed to be fully lost). For example, Ttc21b, which encodes a cilial protein and whose KO is associated with brain, bone and eye phenotypes^20,33, exhibited a marked reduction in retinal neuron (log₂[ratio] = −7.16; Fig. 2b), lens (log₂[ratio] = −2.40) and retina epithelium (log₂[ratio] = −1.65) trajectories (Extended Data Fig. 5a–c). Validations through haematoxylin and eosin staining support these patterns, as the homozygous Ttc21b mutant exhibits a visible collapse in structures that are detectable within the wild-type eye at E13.5. Specifically, the retinal neurons, lens and optic nerve were missing in the homozygous mutant (Fig. 2c). The retinal epithelium was delocalized and reduced as well (Fig. 2c and Extended Data Fig. 5c).

We next examined the positive control, the ZRS limb enhancer KO, a well-studied mutant that shows a loss of the distal limb structure at birth²⁸. Eight sub-trajectories were nominally altered in this mutant, mostly mesenchymal. Although the reduction in limb mesenchymal cells was modest (24% or log₂[ratio] = −0.39), co-embedding of limb mesenchyme cells from ZRS limb enhancer KO and wild-type FVB embryos identified a subpopulation that specifically expressed markers of the distal mesenchyme of the early embryonic limb bud, such as Hoxa13 and Hoxd13, that was markedly affected (Fig. 2d,e and Extended Data Fig. 5d). Such heterogeneity was not observed for the seven other nominally altered sub-trajectories (Extended Data Fig. 5e), consistent with the specificity of this phenotype.

Transcriptional heterogeneity within cell types

We next sought to develop a more sensitive approach for detecting deviations in transcriptional programs within cell-type trajectories. For this, we developed the local cellular heuristic neighbourhood enrichment specificity score (lochNESS), a score calculated on the basis of the ‘neighbourhood’ of each cell in a sub-trajectory co-embedding of a given mutant (all replicates) versus a pooled wild type (all replicates of all backgrounds). Briefly, lochNESS takes aligned PC features of each sub-trajectory and finds k nearest neighbours for each cell from other embryos. For each mutant cell, we compute the fold change of the observed versus expected number of mutant cells in its neighbourhood (Fig. 3a and Methods; similar methods developed independently in ref. ³⁴).

Fig. 3: LochNESS highlights mutant-related changes.

a, A schematic of lochNESS calculation and visualization. b, UMAP of the limb mesenchyme trajectory from ZRS limb enhancer KO and wild-type cells, coloured by lochNESS. The colour scale is centred at the median. Cells corresponding to a subset of ZRS limb enhancer KO cells with more extreme loss in Fig. 2d are outlined with a dashed line. c, Distribution of lochNESS in the neural tube sub-trajectories of Ttc21b-KO and Gli2-KO mutants. The dashed outlines highlight shifted distributions of the retinal neuron sub-trajectory of the Ttc21b-KO mutant and the floor and roof plate sub-trajectories of the Gli2-KO mutant. d, UMAP of co-embedded cells of floor plate and roof plate sub-trajectories from the Gli2-KO mutant and pooled wild type, coloured by lochNESS. e, The same as in d, but coloured by selected marker gene expression. f, Immunofluorescence staining of TTR (ChP marker) in brain regions (LV, lateral ventricle; 4V, fourth ventricle) in sections from the wild type and the Gli2-KO mutant (Methods). Scale bars, 50 μm.

Visualization of lochNESS in the embedded space highlights areas with enrichment or depletion of mutant cells. For example, returning to the ZRS limb enhancer KO embryos, we observe markedly low lochNESS in the distal mesenchyme of the early embryonic limb bud (Figs. 2d and 3b). This highlights the value of lochNESS, as within a sub-trajectory (limb mesenchyme), an effect is both detected and assigned to a subset of cells in a label-agnostic fashion.

Globally, the distribution of lochNESS is unremarkable for some mutants (for example, most TAD boundary KOs) but aberrant for others (for example, pleiotropic mutants such as Sox9 regulatory INV; Extended Data Fig. 6a). After carrying out additional quality control checks (Methods and Extended Data Fig. 6b–d), we examined lochNESS for each mutant in each sub-trajectory. Consistent with earlier analyses, our data show low lochNESS for the retinal neuron sub-trajectory in Ttc21b-KO mice (Fig. 3c and Extended Data Fig. 6e). We also observe a strong shift towards low lochNESS for the floor plate sub-trajectory in Gli2-KO mice, and a subtle change for the roof plate trajectory, which is forming opposite to the floor plate along the dorsal–ventral axis of the developing neural tube³⁵ (Fig. 3c and Extended Data Fig. 6e).

To explore this further, we extracted and reanalysed cells corresponding to the floor plate and roof plate (Extended Data Fig. 7a). Within the floor plate, Gli2-KO cells consistently exhibited low lochNESS (Fig. 3d). However, there were only a handful of differentially expressed genes between the wild type and the mutant, and no significantly enriched pathways. For example, genes such as Robo1 and Slit1, involved in neuronal axon guidance, are specifically expressed in the floor plate relative to the roof plate (Fig. 3e), but are not differentially expressed between wild-type and Gli2-KO cells of the floor plate. Alternatively, our failure to detect substantial differential expression may be due to power, as there were fewer floor plate cells in the Gli2 KO (about 60% reduction). Overall, these observations are consistent with the established role of Gli2 in floor plate induction, its role as an activator of SHH in dorso-ventral patterning of the neural tube and the previous demonstration that Gli2 KOs fail to properly induce a floor plate^35,36.

Less expectedly, we identified two subpopulations of roof-plate-derivative cell types, one depleted and the other enriched in Gli2-KO embryos (Fig. 3d and Extended Data Fig. 7a–c). To annotate these subpopulations, we examined genes whose expression was predicted by lochNESS (Methods). The mutant-enriched group of roof plate cells was marked by cilial genes and Ttr, a marker for the choroid plexus (ChP), whereas the mutant-depleted group was marked by WNT-related genes (for example, Rspo1/2/3 and Wnt3a/8b/9a) indicating it to be a region close to the ChP of the lateral ventricle, namely the cortical hem (Fig. 3e, Extended Data Fig. 7d and Supplementary Tables 4 and 5). We also mapped the three clusters shown in Extended Data Fig. 7a to spatial transcriptomic data from E13.5 mouse embryos³⁷ (Extended Data Fig. 7e). Supporting our annotations, cluster 1 mapped to the floor of the neural tube, cluster 2 next to the lateral ventricle ChP, and cluster 3 to the ChP (both in the lateral (anterior) and fourth (posterior) ventricles). We then examined marker genes that further separate lateral ventricle and fourth ventricle ChP and found that in addition to the roof plate marker Lmx1a, cluster 3 expresses the fourth ventricle marker Meis1 and cluster 2 expresses the lateral ventricle markers Otx1 and Emx2 (Extended Data Fig. 7f and Supplementary Table 4).

To experimentally validate these observations, we examined developmental progression of the neural tube and brain in E13.5 Gli2-KO mutant and wild-type embryos. In coronal sections of the mutant, we observed severe developmental defects including deformed forebrain lobes and delayed neural tube development (Extended Data Fig. 8a). Immunofluorescence imaging of Pax6 expression revealed a severely disturbed shape of the neural tube, confirming the well-described ‘dorsalization’ phenotype of the neural tube (Extended Data Fig. 8b), and consistent with marked reductions in the proportion of floor plate cells in the Gli2-KO mutant (Fig. 3d). Turning to the less expected observation of increased ChP, we found that the lateral ventricle as well as the fourth ventricle exhibited a disturbed pattern of staining of Ttr expression. Whereas the wild type shows inner and outer Ttr signal within the single cell layer, the mutant exhibited a ‘double DAPI’ layer, indicating a disordered tissue organization (Fig. 3f and Extended Data Fig. 8c,d). Adjusting for the overall smaller size of Gli2-KO mutants at E13.5, we quantified cells positive for Ttr expression in the lateral and fourth ventricle, and found a proportional increase in the mutant relative to the wild type (Supplementary Table 6), again consistent with the marked increase in the proportion of ChP cells in this mutant (Fig. 3d). In summary, we could confirm both the expected reduction in floor plate and the unexpected increase in roof-plate-derived ChP in the mutant. Of note, the relatively subtle and opposing effects on these roof plate subpopulations were missed by our original analysis of cell-type proportions, and uncovered only by the granularity of lochNESS.

LochNESS distributions can be systematically screened to identify sub-trajectories exhibiting mutant-specific shifts. For example, although all TAD boundary KO mutants have similarly unremarkable global lochNESS distributions, when we plot these distributions by sub-trajectory, a handful of shifted distributions are evident (Extended Data Fig. 9a,b). For example, multiple epithelial sub-trajectories, including pre-epidermal keratinocyte, epidermis, branchial arch and lung epithelial trajectories, are most shifted in Tbx3 TAD boundary KO cells, with further analyses preliminarily supporting a role for Tbx3 in epidermal and lung development³⁸ (Methods, Extended Data Fig. 9c,d and Supplementary Table 7).

Mutant-specific and mutant-shared effects

Pleiotropy, wherein a single gene influences multiple, unrelated traits, is a pervasive phenomenon in developmental genetics, and yet remains poorly understood³⁹. Although here we have ‘whole-embryo’ molecular profiling of just 22 mutants, we sought to investigate whether we could distinguish between mutant-specific and mutant-shared effects within each major trajectory. In brief, within a co-embedding of cells from all embryos from a given background, we computed k nearest neighbours as in Fig. 3a, and then calculated the observed versus expected ratio of each genotype among a cell’s k nearest neighbours. The ‘similarity score’ between one genotype versus all others is defined as the mean of these ratios across cells of the genotype (Methods). To assess whether any observed similarities or dissimilarities are robust, we can also calculate similarity scores between individual embryos. For example, for the mesenchymal trajectory of C57BL/6 mutants, similarity scores are generally higher for pairwise comparisons of individuals with the same genotype (Fig. 4a and Extended Data Fig. 10a–c). Pairs of individuals with the Scn11a-GOF mutation exhibited the most extreme similarity scores, consistent with our earlier observation that they clustered with E12.5 rather than E13.5 embryos (Fig. 1c). Following further analysis, we believe that the most parsimonious explanation is incorrect staging of these litters, rather than mutation-specific, global developmental delay (Extended Data Fig. 10d–g and Supplementary Note 1).

Fig. 4: Similarity scores identify mutant-shared and mutant-specific effects.

a, A heatmap showing similarity scores between individual C57BL/6 embryos in the mesenchymal trajectory. The rows and columns are grouped by genotype and labelled by embryo id and genotype. b, UMAP of the lateral plate and intermediate mesoderm sub-trajectory for mutants from the C57BL/6 background strain, coloured and labelled by subcluster and detailed cell type (marker genes in Supplementary Table 4). c, Boxplots showing the composition of the top six subclusters for individual Atp6v0a2-KO, Atp6v0a2^R755Q and Gorab-KO embryos (blue, n = 4 each genotype) and other C57BL/6 embryos (red, n = 23). Centre lines show medians; box limits indicate 25th and 75th percentiles; replicates are represented by dots. d, The same as in b, but coloured by log-transformed expression of selected marker genes.

We also observed that the similarity scores between three mutants (Atp6v0a2 KO, Atp6v0a2^R755Q and Gorab KO) were consistent with shared effects, in the mesenchymal, epithelial, endothelial, hepatocyte and neural crest (peripheral nervous system glia) trajectories in particular; in other major trajectories, such as neural tube and haematopoiesis, the Atp6v0a2 KO and Atp6v0a2^R755Q exhibited high similarity scores with one another, but not with the Gorab KO (Fig. 4a and Extended Data Fig. 10a,c,f). In human patients, mutations in ATP6V0A2 and GORAB cause overlapping connective tissue disorders, which is reflected in the misregulation of the mesenchymal trajectory of Atp6v0a2 and Gorab mutants^25,26. However, only the ATP6V0A2-related disorder exhibits a prominent central nervous system phenotype, consistent with the changes in the neural tube trajectory seen only in Atp6v0a2 mutants (Extended Data Fig. 10a,c,f).

To further explore phenotypic sharing between these mutants, we co-embedded cells of the lateral plate and intermediate mesoderm sub-trajectory from C57BL/6 strains. We resolved the identity of most subclusters using marker genes and spatial mapping, identifying multiple subsets for which Atp6v0a2-KO, Atp6v0a2^R755Q and Gorab-KO mice are similarly distributed compared to other C57BL/6 genotypes (Fig. 4b and Extended Data Fig. 11). Some subsets are enriched for cells from these mutants (for example, proepicardium, hepatic mesenchyme and lung mesenchyme) whereas others are depleted (for example, gastrointestinal smooth muscle; Fig. 4c,d and Supplementary Table 4). Although individually subtle, the consistent shifts in cell-type proportions between the two Atp6v0a2- and Gorab-KO mutants across these subsets of mesenchyme derived from lateral plate mesoderm presumably underlie their high mesenchymal similarity scores (Fig. 4c).

Altogether, these analyses illustrate how the joint analysis of mutants subjected to whole-embryo scRNA-seq can reveal sharing of molecular and cellular phenotypes. This includes global similarity (Atp6v0a2 KO versus Atp6v0a2^R755Q) as well as instances in which specific aspects of phenotypes are shared between previously unrelated mutants (Atp6v0a2 mutants versus Gorab KO).

Mesenchymal stalling in a Sox9 regulatory mutant

About half of the mutants profiled here model disruptions of regulatory, rather than coding, sequences. Among these, the Sox9 regulatory INV mutant stands out in having a markedly shifted lochNESS distribution, particularly in mesenchyme (Fig. 5a and Extended Data Fig. 6a). Sox9 encodes a pleiotropic transcription factor crucial for development of the skeleton, the brain, sex determination and other systems, orchestrated by a complex regulatory landscape^40,41,42. This particular mutant features an inversion of a 1-megabase upstream region bearing several distal enhancers and a TAD boundary, essentially relocating these elements into a TAD with Kcnj2, which encodes a potassium channel⁵ (Fig. 5b). Like the Sox9 KO, the homozygous Sox9 regulatory INV is perinatally lethal, with extensive skeletal phenotypes including digit malformation, a cleft palate, bowing of bones and delayed ossification. In addition to the loss of 50% of Sox9 expression, the inversion causes pronounced misexpression of Kcnj2 in the digit anlagen in a wild-type Sox9 pattern⁵. However, the extent to which Kcnj2 and Sox9 are misexpressed elsewhere, as well as the molecular and cellular correlates of the widespread skeletal phenotype, have yet to be deeply investigated.

Fig. 5: Apparent stalling and redirection of mesenchyme differentiation in the Sox9 regulatory INV mutant.

a, LochNESS distributions for all G4 mutants in the mesenchymal trajectory (left) and the Sox9 regulatory INV mutant in mesenchymal sub-trajectories (right). b, Model of the Sox9 regulatory INV mutation depicting ectopic Kcnj2 regulation through enhancer adoption. c, RNA velocity of mesenchymal G4 wild-type and Sox9 regulatory INV cells coloured by sub-trajectories (top) or genotype (bottom). d, Sox9 regulatory INV heterozygous mutant and littermate wild-type RNAscope images (red: Kcnj2; green: Sox9), with insets below highlighting a region corresponding to developing bone (area outlined with white dots line). Scale bars, 500 μm. e, RNA velocity of G4 wild-type and Sox9 regulatory INV cells in the limb mesenchymal sub-trajectory labelled by annotation (top) or genotype (bottom). f, The same as in e, but coloured by branch number. g, Dot plot of the top differentially expressed genes in the four branches shown in f.

At the level of mesenchymal sub-trajectories, shifts in the lochNESS distribution for Sox9 regulatory INV were consistently observed, but limb mesenchyme and connective tissue were particularly enriched for cells with extremely high lochNESS (Fig. 5a, right). Notably, two of three major enhancers (E250 and E195) known to drive Sox9-mediated chondrogenesis in mesenchymal stem cells are located within the inverted region⁴⁰ (Fig. 5b). Cell-type composition analysis (Fig. 2a) showed that Sox9 regulatory INV mutants harbour considerably larger numbers of cells classified as limb mesenchyme, at the expense of osteoblasts, lateral plate and intermediate mesoderm, chondrocytes and connective tissue trajectories. This shift can also be seen in a uniform manifold approximation and projection (UMAP) embedding (Fig. 5c), a topic we revisit further below.

These changes in cell-type composition were accompanied by reduced expression of Sox9 and increased expression of Kcnj2 in bone (aggregate of chondrocyte, osteoblast and limb mesenchyme; Extended Data Fig. 12a), although the number of cells expressing Kcnj2 was generally low. This suggests that the Sox9 regulatory inversion is resulting in increased Kcnj2 expression (through Sox9 enhancer adoption) and Sox9 reduction (through boundary repositioning) not only in the digit anlagen, but in skeletal mesenchyme more generally. To validate this, we carried out RNA in situ hybridization (RNAscope) on sections of developing bones of the rib cage at E13.5, comparing a heterozygous Sox9 regulatory INV mouse with a wild-type littermate. Consistent with our scRNA-seq data derived from homozygous mutants, our data show a Sox9-patterned increase in Kcnj2 levels, together with losses in Sox9 expression, in the developing bone (Fig. 5d and Extended Data Fig. 12b).

As the inverted Sox9 regulatory region also hosts multiple enhancers active in other tissues (for example, E161 in lung and E239 in cerebral cortex)⁴⁰, we wondered whether these patterns were also seen in other tissues. Indeed, both scRNA-seq and RNAscope quantification show increased Kcnj2 levels in all other tissues examined. Whereas reductions in Sox9 expression, clear in bone, were not observed in most other tissues by scRNA-seq, RNAscope showed Sox9 reductions in telencephalon and lung as well (Extended Data Fig. 12a,b). Taken together, these data suggest marked changes in mesenchyme due to reduced Sox9, together with broader increases in Kcnj2 expression. As expected on the basis of the role of Sox9 in chondrogenesis, hallmark pathways related to chondrocyte proliferation and differentiation^43,44,45,46 were downregulated; less expectedly, several immune-related pathways were upregulated (Extended Data Fig. 12c).

To explore the apparent accumulation of limb mesenchyme in the Sox9 regulatory INV (Fig. 5c and Extended Data Fig. 12d) in more detail, we reanalysed mutant and wild-type cells from the limb mesenchyme sub-trajectory, which revealed subpopulations of condensing mesenchyme, perichondrium and undifferentiated mesenchyme (Extended Data Fig. 12e,f). RNA velocity analyses suggested that most limb mesenchyme ‘accumulation’ in mutant embryos is due to cells that are delayed or stalled in an undifferentiated or stem-like state (Fig. 5c,e and Extended Data Fig. 12e). This accumulation is even more apparent in integrated views of the limb mesenchyme sub-trajectory, for which we observe branches that are highly enriched for Sox9 regulatory INV mutant cells, within undifferentiated mesenchyme (Fig. 5e and Extended Data Fig. 12g,h).

To investigate these branches further, we subclustered undifferentiated mesenchyme cells from the mutant and wild type (Fig. 5f,g). Notably, the most differentially expressed genes in ‘branch 2’ were largely neuronal (for example, several neurexins and neuregulin 3), an observation supported by gene set enrichment analysis (Extended Data Fig. 12i,j). A cellular composition analysis revealed that these neuronal-like cells were not restricted to the Sox9 regulatory INV mutant, but also found in wild-type embryos, albeit much less frequently (Extended Data Fig. 12g,h). To validate this unexpected ‘neural-like’ branch of mesenchymal cells as well as to assess their anatomical distribution, we mapped these cells to spatial transcriptomic data from E13.5 mouse embryos³⁷. Strikingly, this analysis placed branch 2 cells along the neural tube and the brain regions (Extended Data Fig. 13a). To address concerns that artefacts might arise from mapping single-cell data onto non-single-cell spatial maps, we also integrated our data with sci-space⁴⁷ spatial transcriptomic data (E14.5), as these retain single-nucleus resolution. The results are consistent, in that branch 2 mesenchymal cells are enriched in brain regions, branch 0 cells are enriched in limb bud regions, and branch 1 and 3 cells are diffusely distributed but largely excluded from brain regions (Extended Data Fig. 13b).

Taken together, these analyses support the validity of this neural-like subset of mesenchyme (present in the wild type and increased in Sox9 regulatory INV mutants). The observation is consistent with the reports that mesenchymal stem cells can be differentiated to neuronal states in vitro⁴⁸.

Discussion

Here we set out to establish whole-embryo scRNA-seq as a new paradigm for the systematic, scalable phenotyping of mouse developmental mutants. On data obtained for 22 mutants in a single experiment, we developed analytical approaches to identify deviations in cell-type composition, subtle differences in gene expression within cell types (lochNESS), and sharing of sub-phenotypes between mutants (similarity scores). Overall, the results are encouraging, and show how systematic, outcome-agnostic computational analyses of data obtained at the whole-embryo scale may in some cases reveal molecular and cellular phenotypes that are missed by conventional phenotyping.

We emphasize that the concurrent analysis of many mutants proved essential to the contextualization of particular observations (that is, to understand how specific or nonspecific any apparent deviation really was) against a background of dozens of genotypes and more than 100 embryos. This also enabled us to discover shared aspects of phenotypes between previously unrelated genotypes (for example, between Gorab and Atp6v0a2 mutants). Looking forward, profiling of additional mouse mutants might enable the further ‘decomposition’ of developmental pleiotropy, a poorly understood phenomenon, into ‘basis vectors’.

The diverse mutants analysed yielded a variety of results that speak to the utility of whole-embryo scRNA-seq for phenotyping. For example, an abnormal eye phenotype in Ttc21b mutants was previously described, but considered probably to be secondary to a more general craniofacial defect^20,33. The scRNA-seq analysis of E13.5 Ttc21b mutants demonstrated that multiple retinal cell trajectories were essentially absent. Detailed histological analysis confirmed this, suggesting that the eye abnormality is probably not a secondary effect, but rather that the overactive SHH signalling has a primary effect on retinal development in this mutant.

The utility of pursuing whole-embryo scRNA-seq was also demonstrated by an unexpected finding of both a depleted and an enriched cell population of roof plate cell derivatives in the Gli2-KO mutant. The ‘dorsalization’ of the neural tube in the absence of SHH signalling is well described^20,35,36 and was confirmed in our histological analysis of this line (Extended Data Fig. 8). However, there have been no described changes in the roof plate or its derivatives so far in Gli2-KO mice³⁶. By contrast, whole-embryo scRNA-seq uncovered that derivatives of the roof plate depict changes in composition (a primary finding) and tissue development (a finding based on secondary validation) in the mutant, illustrating how this approach can potentially yield new insight into even well-studied developmental pathways. However, owing to our dataset capturing only one time point, whether Gli2 misexpression causes the structural change directly in the derivative tissue or earlier during roof plate formation remains elusive.

Our MMCA has limitations. First, we profiled only four replicates per mutant at a single developmental time point. On the basis of a simulation analysis of the analytical approach that considers only cell proportions, four replicates of each mutant is probably sufficient to detect modest changes in abundant cell types (for example, a 10% change for cell types at 10% abundance) but only large changes in rarer cell types (for example, a 25% change in cell types at 1% abundance; Extended Data Fig. 4b). As such, to detect more subtle changes in model organisms such as mice for which very large numbers of replicates are not feasible, more sophisticated strategies such as lochNESS, which is not based on counts of cell types but rather directly considers the distribution of cells derived from different genotypes in a complex embedding, may be essential. It is important to note that our cell composition analysis, which includes both wild-type and mutant cells from the same strain to generate a pooled reference, assumes that the cell-type proportions of non-wild-type genotypes are roughly consistent, at least on the whole, with those of wild-type cells. This assumption may be more problematic in studies of biologically related mutants. Of note, in concurrently published studies in this issue, a similar approach is taken for genetic and environmental perturbations in zebrafish (ref. ⁴⁹), such that dozens to hundreds of replicate embryos of each genotype can be profiled and phenotypic variability quantified.

Second, profiling only a small fraction of cells present in E13.5 embryos potentially limits sensitivity. However, for any given mutant, we had more than 1.5 million cells from other genotypes (wild type or other mutants), which facilitated the detection of mutant-specific phenotypes for rare cell types (for example, in the retina (Ttc21b KO) and roof plate (Gli2 KO)).

Third, we were not able to explore all mutants in detail, nor to thoroughly investigate other aspects of the data (for example, the differences between wild-type strains). In the future, we anticipate that community input and domain expertise will be essential to extract full value from these data. To facilitate this, we created an interactive browser that allows exploration of mutant-specific effects on gene expression in trajectories and sub-trajectories, together with the underlying data (https://atlas.gs.washington.edu/mmca_v2/). Additionally, some of the phenotypes identified here have probably not been described before owing to the lack of resolution of conventional phenotyping. New secondary validation strategies need to be developed to confirm subtle defects in molecular programs or subtle changes in the relative proportions of specific cell types. A promising approach would be to complement whole-embryo scRNA-seq with rapidly advancing methods for whole-mouse-body antibody labelling and three-dimensional imaging⁵⁰.

Fourth, our results emphasize the importance of a well-matched control; although data from our wild-type embryos could be reused as control data for future studies of additional mutants, that risks batch effects, and a safer strategy would be to always include a well-matched, ‘in-line’ wild-type control while profiling mutant embryos.

In 2011, the International Mouse Phenotyping Consortium set out to drive towards the ‘functionalization’ of every protein-coding gene in the mouse, by generating thousands of KO mouse lines⁵¹. In principle, the whole-embryo scRNA-seq phenotyping approach presented here could be extended to all Mendelian genes or even to all 20,000 mouse gene KOs.

Methods

Data reporting

No statistical methods were used to predetermine sample size. Embryos used in experiments were randomized before sample preparation. Investigators were blinded to group allocation during data collection and analysis. Embryo collection and sci-RNA-seq3 analysis were carried out by different researchers in different locations.

Embryo collection

Mutants were generated through conventional gene-editing tools and breeding or tetraploid aggregation and collected at the embryonic stage E13.5, calculated from the day of vaginal plug (noon = E0.5). Collection and whole-embryo dissection were carried out as previously described⁵². The embryos were immediately snap-frozen in liquid nitrogen and shipped to the Shendure Lab (University of Washington) in dry ice. Sets of animals with the same genotype were either all male or half male–half female. All animal procedures were carried out in accordance with institutional, state and government regulations.

Isolation and fixation of nuclei

Snap-frozen embryos were processed as previously described³. Briefly, the frozen embryos were cut into small pieces with a blade and further dissected by resuspension in 1 ml ice-cold cell lysis buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl₂, 0.1% IGEPAL CA-630, 1% SUPERase In and 1% BSA) in a 6-cm dish. After adding another 3 ml of cell lysis buffer, the sample was strained (40 µm) into a 15-ml Falcon tube and centrifuged to a pellet (500g, 5 min). By resuspending the sample with another 1 ml of cell lysis buffer, the isolation of nuclei was ensured. The nuclei were pelleted again (500g, 5 min) and then washed and fixed in 10 ml 4% paraformaldehyde (PFA) for 15 min on ice. The fixed nuclei were pelleted (500g, 3 min) and washed twice in the nucleus suspension buffer (500g, 5 min). The nuclei finally were resuspended in 500 µl nucleus suspension buffer and split into two tubes, each containing 250 µl of sample. The tubes were flash frozen in liquid nitrogen and stored in a −80 °C freezer, until further use for library preparation. The embryo preparation was carried out randomly for nuclei isolation to avoid batch effects.

sci-RNA-seq3 library preparation and sequencing

The library preparation was carried out as previously described⁵³. In short, the fixed nuclei were permeabilized, sonicated and washed. Nuclei from each mouse embryo were then distributed into several individual wells in four 96-well plates. We split samples into four batches (about 25 samples randomly selected in each batch) for sci-RNA-seq3 processing. The ID of the reverse transcription well was linked to the respective embryo for downstream analysis. In a first step, the nuclei were then mixed with oligo-dT primers and dNTP mix, denatured and placed on ice; afterwards, they were processed for reverse transcription including a gradient incubation step. After reverse transcription, the nuclei from all wells were pooled with the nuclei dilution buffer (10 mM Tris-HCl, pH 7.4, 10 mM NaCl, 3 mM MgCl₂, 1% SUPERase In and 1% BSA), spun down and redistributed into 96-well plates containing the reaction mix for ligation. The ligation proceeded for 10 min at 25 °C. Afterwards, nuclei again were pooled with nuclei suspension buffer, spun down and washed and filtered. Next, the nuclei were counted and redistributed for second strand synthesis, which was carried out at 16 °C for 3 h. Afterwards, tagmentation mix was added to each well, and tagmentation was carried out for 5 min at 55 °C. To stop the reaction, DNA binding buffer was added and the sample was incubated for another 5 min. Following an elution step using AMPure XP beads and elution mix, the samples were subjected to PCR amplification to generate sequencing libraries.

Finally after PCR amplification, the resulting amplicons were pooled and purified using AMPure XP beads. The library was analysed by electrophoresis and the concentration was calculated using Qubit (Invitrogen). The library was sequenced on the NovaSeq platform (Illumina; read 1: 34 cycles, read 2: 100 cycles, index 1: 10 cycles, index 2: 10 cycles).

Processing of sequencing reads

Read alignment and cell × gene expression count matrix generation was carried out on the basis of the pipeline that we developed for sci-RNA-seq3 (ref. ³) with the following minor modifications: base calls were converted to fastq format using Illumina’s bcl2fastq v2.20 and demultiplexed on the basis of PCR i5 and i7 barcodes using the maximum-likelihood demultiplexing package deML⁵⁴ with default settings. Downstream sequence processing and cell × gene expression count matrix generation were similar to sci-RNA-seq⁵⁵ except that the reverse transcription (RT) index was combined with the hairpin adaptor index, and thus the mapped reads were split into constituent cellular indices by demultiplexing reads using both the RT index and the ligation index (Levenshtein edit distance (ED) < 2, including insertions and deletions). Briefly, demultiplexed reads were filtered on the basis of the RT index and ligation index (ED < 2, including insertions and deletions) and adaptor-clipped using trim_galore v0.6.5 with default settings. Trimmed reads were mapped to the mouse reference genome (mm10), using STAR v2.6.1d⁵⁶ with default settings and gene annotations (GENCODE VM12 for mouse). Uniquely mapping reads were extracted, and duplicates were removed using the unique molecular identifier (UMI) sequence (ED < 2, including insertions and deletions), reverse transcription (RT) index, hairpin ligation adaptor index and read 2 end-coordinate (that is, reads with UMI sequence less than 2 ED, RT index, ligation adaptor index and tagmentation site were considered duplicates). Finally, mapped reads were split into constituent cellular indices by further demultiplexing reads using the RT index and ligation hairpin (ED < 2, including insertions and deletions). To generate the cell-x-gene expression count matrix, we calculated the number of strand-specific UMIs for each cell mapping to the exonic and intronic regions of each gene with Python v2.7.13 HTseq package⁵⁷. For multi-mapped reads, reads were assigned to the closest gene, except in cases in which another intersected gene fell within 100 base pairs of the end of the closest gene, in which case the read was discarded. For most analyses, we included both expected-strand intronic and exonic UMIs in the cell-x-gene expression count matrix.

The single-cell gene count matrix included 1,941,605 cells after cells with low quality (UMI ≤ 250 or detected gene ≤ 100) were filtered out. Each cell was assigned to its original mouse embryo on the basis of the reverse transcription barcode. We applied three strategies to detect potential doublet cells. As the first strategy, we split the dataset into subsets for each individual, and then applied the scrublet v0.1 pipeline⁵⁸ to each subset with parameters (min_count = 3, min_cells = 3, vscore_percentile = 85, n_pc = 30, expected_doublet_rate = 0.06, sim_doublet_ratio = 2, n_neighbors = 30, scaling_method = ‘log’) for doublet score calculation. Cells with doublet scores above 0.2 were annotated as detected doublets (5.5% in the whole dataset).

As the second strategy, we used an iterative clustering strategy based on Seurat v3 (ref. ⁵⁹) to detect the doublet-derived subclusters for cells. Briefly, gene count mapping to sex chromosomes was removed before clustering and dimensionality reduction, and then genes with no count were filtered out and each cell was normalized by the total UMI count per cell. The top 1,000 genes with the highest variance were selected. The data was log-transformed after adding a pseudocount, and scaled to unit variance and zero mean. The dimensionality of the data was reduced by PCA (30 components) first and then with UMAP, followed by Louvain clustering carried out on the 10 PCs (resolution = 1.2). For Louvain clustering, we first fitted the top 10 PCs to compute a neighbourhood graph of observations (k.param = 50) followed by clustering the cells into subgroups using the Louvain algorithm. For UMAP visualization, we directly fitted the PCA matrix with min_distance = 0.1. For subcluster identification, we selected cells in each major cell type and applied PCA, UMAP and Louvain clustering similarly to the major cluster analysis. Subclusters with a detected doublet ratio (by Scrublet) above 15% were annotated as doublet-derived subclusters.

We found that the above Scrublet and iterative clustering-based approach is limited in marking cell doublets between abundant cell clusters and rare cell clusters (for example, less than 1% of the total cell population); thus, we applied a third strategy to further detect such doublet cells. Briefly, cells labelled as doublets (by Scrublet) or from doublet-derived subclusters were filtered out. For each cell, we retained only protein-coding genes, long intergenic noncoding RNA genes and pseudogenes. Genes expressed in fewer than 10 cells and cells expressing fewer than 100 genes were further filtered out. The downstream dimension reduction and clustering analysis were carried out with Monocle v3 (ref. ³). The dimensionality of the data was reduced by PCA (50 components) first on the top 5,000 most highly variable genes and then with UMAP (max_components = 2, n_neighbors = 50, min_dist = 0.1, metric = ‘cosine’). Cell clusters were identified using the Leiden algorithm implemented in Monocle v3 (resolution = 1 × 10⁻⁶). Next, we took the cell clusters identified by Monocle v3 and first computed differentially expressed genes across cell clusters with the top_markers function of Monocle v3 (reference_cells = 1,000). We then selected a gene set combining the top 10 gene markers for each cell cluster (filtering out genes with fraction_expressing <0.1 and then ordering by pseudo_R2). Cells from each main cell cluster were selected for dimension reduction by PCA (10 components) first on the selected gene set of top cluster-specific gene markers, and then by UMAP (max_components = 2, n_neighbors = 50, min_dist = 0.1, metric = ‘cosine’), followed by clustering identification using the Leiden algorithm implemented in Monocle v3 (resolution = 1 × 10⁻⁴). Subclusters showing low expression levels of markers specific for target cell clusters and enriched expression levels of markers specific for non-target cell clusters were annotated as doublet-derived subclusters and filtered out in visualization and downstream analysis. Finally, after removing the potential doublet cells detected by either of the above three strategies, 1,671,270 cells were retained for further analyses.

Whole-mouse-embryo analysis

As described previously³, each cell could be assigned to the mouse embryo from which it derived on the basis of its reverse transcription barcode. After removing doublet cells and another 25 cells that were poorly assigned to any mouse embryo, 1,671,245 cells from 103 individual mouse embryos were retained (a median of 13,468 cells per embryo). UMI counts mapping to each sample were aggregated to generate a pseudobulk RNA-seq profile for each sample. Each cell’s counts were normalized by dividing them by the estimated size factor, and then the data were log₂-transformed after adding a pseudocount followed by carrying out the PCA. The normalization and dimension reduction were carried out in Monocle v3.

We previously used sci-RNA-seq3 to generate the MOCA dataset, which profiled about 2 million cells derived from 61 wild-type B6 mouse embryos staged between stages E9.5 and E13.5. The cleaned dataset, including 1,331,984 high-quality cells, was generated by removing cells with <400 detected UMIs as well as doublets (http://atlas.gs.washington.edu/mouse-rna). UMI counts mapping to each sample were aggregated to generate a pseudobulk RNA-seq profile for each embryo. Each cell’s counts were normalized by dividing them by the estimated size factor, and then the data were log₂-transformed after adding a pseudocount, followed by PCA. The PCA space was retained, and then the embryos from the MMCA dataset were projected onto it.

Cell clustering and annotation

After removing doublet cells, genes expressed in fewer than 10 cells and cells expressing fewer than 100 genes were further filtered out. We also filtered out low-quality cells on the basis of the proportion of reads mapping to the mitochondrial genome (MT%) or ribosomal genome (Ribo%) (specifically, filtering cells with MT% > 10 or Ribo% > 5). We then removed cells from two embryos that were identified as outliers on the basis of the whole-mouse-embryo analysis (embryo 41 and embryo 104). This left 1,627,857 cells (median UMI count 845; median genes detected 539) from 101 individual embryos that were retained for all subsequent analyses.

To eliminate the potential heterogeneity between samples due to different mutant types and genotype backgrounds, we sought to carry out the dimensionality reduction on a subset of cells from the wild-type mice (including 15 embryos with 215,575 cells, 13.2% of all cells) followed by projecting all remaining cells, derived from the various mutant embryos, onto this same embedding. These procedures were carried out using Monocle v3. In brief, the dimensionality of the subset of data from the wild-type mice was reduced by PCA, retaining 50 components, and all remaining cells were projected onto that PCA embedding space. Next, to mitigate potential technical biases, we combined all cells from wild-type and mutant mice and applied the align_cds function implemented in Monocle v3, with MT%, Ribo% and log-transformed total UMI of each cell as covariates. We took the subset of cells from wild-type mice, using their ‘aligned’ PC features to carry out UMAP (max_components = 3, n_neighbors = 50, min_dist = 0.01, metric = ‘cosine’) by uwot v0.1.8, followed by saving the UMAP space. Cell clusters were identified using the Louvain algorithm implemented in Monocle v3 on three dimensions of UMAP features, resulting in 13 isolated major trajectories (Fig. 1e). We then projected all of the remaining cells from mutant mouse embryos onto the previously saved UMAP space and predicted their major-trajectory labels using a k-nearest-neighbour (k-NN) heuristic. Specifically, for each mutant-derived cell, we identified its 15 nearest-neighbour wild-type-derived cells in UMAP space and then assigned the major trajectory with the maximum frequency within that set of 15 neighbours as the annotation of the mutant cell. We calculated the ratio of the maximum frequency to the total as the assigned score. Of note, more than 99.9% of the cells from the mutant mice had an assigned score greater than 0.8. The cell-type annotation for each major trajectory was based on expression of the known marker genes (Supplementary Table 2).

Within each major trajectory, we repeated a similar strategy, but with slightly adjusted PCA and UMAP parameters. For the major trajectories with more than 50,000 cells, we reduced the dimensionality by PCA to 50 PCs; for the other major trajectories of more than 1,000 cells, we reduced the dimensionality by PCA to 30 PCs; for the remaining major trajectories, we reduced the dimensionality by PCA to 10 PCs. UMAP was carried out with max_components = 3, n_neighbors = 15, min_dist = 0.1, metric = ‘cosine’. For the mesenchymal trajectory, we observed a notable separation of cells by their cell-cycle phase in the UMAP embedding. We calculated a g2m index and an s index for individual cells by aggregating the log-transformed normalized expression for marker genes of the G2M phase and the S phase and then included them in the align_cds function along with the other factors. Applying these procedures to all of the major trajectories, we identified 64 sub-trajectories in total. Similarly, after assigning each cell from the mutant mice with a sub-trajectory label, we calculated the ratio of the maximum frequency to the total as the assigned score. Of note, more than 96.7% of the cells from the mutant mice had an assigned score greater than 0.8. The cell-type annotation for each sub-trajectory was also based on the expression of known marker genes (Supplementary Table 2).

Identification of correlated cell trajectories between datasets

To identify correlated cell trajectories between MOCA and MMCA datasets, we first calculated an aggregate expression value for each gene in each cell trajectory by summing the log-transformed normalized UMI counts of all cells of that trajectory. For consistency during the comparison to MOCA, we manually regrouped the cells from the MMCA dataset into 10 cell trajectories, by merging the olfactory sensory neuron trajectory into the neural crest (peripheral nervous system neuron) trajectory, merging the myotube trajectory, the myoblast trajectory and the cardiomyocyte trajectory into the mesenchymal trajectory, and splitting the hepatocyte trajectory into the lens epithelial trajectory and the liver hepatocyte trajectory. Next, for the two datasets, we applied non-negative least-squares regression to predict gene expression in a target trajectory (T_a) in dataset A based on the gene expression of all trajectories (M_b) in dataset B: T_a = β_0a + β_1aM_b, based on the union of the 3,000 most highly expressed genes and 3,000 most highly specific genes in the target trajectory. We then switched the roles of datasets A and B; that is, predicting the gene expression of the target trajectory (T_b) in dataset B from the gene expression of all trajectories (M_a) in dataset A: T_b = β_0b + β_1bM_a. Finally, for each trajectory a in dataset A and each trajectory b in dataset B, we combined the two correlation coefficients: β = 2(β_ab + 0.001)(β_ba + 0.001) to obtain a statistic, for which high values reflect reciprocal, specific predictivity. We repeated this analysis on sub-trajectories within each major trajectory.

Identification of significant cell composition changes in mutant mice using beta-binomial regression

A cell number matrix of all 64 developmental sub-trajectories (rows) and 101 embryos (columns) was created and the cell numbers were then normalized by the size factor of each column that was estimated by the estimate_size_factors function in Monocle v3. Ten sub-trajectories with a mean cell number across individual embryos <10 were filtered out. The beta-binomial regression was carried out using the VGAM package of R. The following code was used: vglm(cbind(n_celltype, n_total - n_celltype) ~ genotype, family = betabinomial), where n_celltype refers to the trajectory-specific cell number, and n_total refers to the total cell number of that embryo. Of note, embryos from the four different mouse strain backgrounds were analysed independently.

We reason that the power of our strategy to detect the cell proportion changes between different genotypes is affected by three factors: the abundance of a given cell type; the number of replicates in each genotype group; and the effect size. To evaluate power, we carried out a simulation analysis that varied these factors, implemented as follows.

• 1. We selected the 20 most abundant cell types in wild-type embryos. Their abundances ranged from about 1% to about 20%. The proportions of these cell types served as the basis for our simulations.

• 2. We simulated ten groups of ‘wild-type’ samples with 4, 8, 16, …, 40 replicates in each group, wherein each sample consisted of cells drawn from the 20 cell types. For each replicate, the simulated number of cells of each cell type was calculated as the product of: (a) the cell-type proportions, simulated by fitting a Dirichlet model based on the real proportions from step 1; and (b) the total number of cells recovered for that replicate, simulated on the basis of the mean (n ≈ 15,000) and standard deviation of the cell numbers across replicates in the real dataset.

• 3. We simulated ten groups of ‘mutant’ samples by repeating the above step except adding shifts to the numbers of cells within each cell type. The shifting scales were based on different effect sizes. For instance, effect size =  0.1 represents a 10% reduction in the number of cells.

• 4. We carried out beta-binomial regression (the same test used in Fig. 2a) to test whether the cell-type proportions were significantly changed between simulated ‘wild-type’ and ‘mutant’ samples, further checking the results as stratified by cell type (with different abundances), the number of replicates and the effect size.

The results are in line with our hypothesis that the detection power of our strategy varies among comparisons with different effect sizes, sample sizes or cell-type abundances (Extended Data Fig. 4). The main ‘take-home’ messages are summarized below.

• 1. Changes of 25% are robustly detectable, even for rare cell types (for example, <2%), with modest numbers of replicates.

• 2. Changes of 10% are possible to detect, but only for abundant cell types (for example, >5%). More replicates can help in this zone.

• 3. Changes of 1% are almost impossible to detect with a cell proportion approach, even with very large numbers of replicates.

In general, at the level of single-cell sampling carried out in our study, four samples (corresponding to the number of samples used in the manuscript) would be sufficient to detect a 25% effect size for those cell types present at a 1% proportion in wild-type embryos.

Defining and calculating lochNESS

To identify local enrichments or depletions of mutant cells, we aim to define a metric for each single cell to quantify the enrichments or depletions of mutant cells in its surrounding neighbourhood. For these analyses, we consider a mutant and a pooled wild type combining all four background strains in a major trajectory as a dataset. For each dataset, we define lochNESS as:

$${\rm{lochNESS}}=\frac{{\rm{number\; of\; mutant\; cells\; in}}\,k\, \mbox{-} \,{\rm{NNs}}}{k}/\frac{{\rm{number\; of\; mutant\; cells\; in\; dataset}}}{N}-1,$$

in which N is the total number of cells in the dataset, $k=\frac{\sqrt{N}}{2}$ scales with N, and the cells from the same embryo as the cell of interest are excluded from the k-NNs. Note that this value is equivalent to the fold change of mutant cell percentage in the neighbourhood of a cell relative to in the whole major trajectory. For implementation, we took the aligned PCs in each sub-trajectory as calculated above, and for each cell in an embryo we find the k-NNs in the remaining mutant embryo cells and wild-type cells. We plot the lochNESS in a red–white–blue scale, for which white corresponds to 0 or the median lochNESS, blue corresponds to high lochNESS or enrichments, and red corresponds to low lochNESS or depletions.

At present, we calculate lochNESS using a pooled wild type combining all four background strains to include larger numbers of cells in constructing the k-NN graph. If the numbers of cells are sufficient, a wild type from the matched background strain can be used. Additionally, if the numbers of cells are sufficient, one set of lochNESS can be calculated for each wild-type sample separately and the variability between samples can be considered.

Examining global distributions of lochNESS

Plotting the global distributions of lochNESS for each mutant across all sub-trajectories, we further observed that some mutants (for example, most TAD boundary KOs; Scn11a GOF) exhibit unremarkable distributions (Extended Data Fig. 6a). However, others (for example, Sox9 regulatory INV; Scn10a/11a double KO) are associated with a marked excess of high lochNESS, consistent with mutant-specific effects on transcriptional state across many developmental systems. For reference, we simultaneously create a null distribution of lochNESS using random permutation of the mutant and wild-type cell labels, simulating datasets in which the cells are randomly mixed. Of note, we confirmed that repeating the calculation of lochNESS after random permutation of mutant and wild-type labels resulted in bell-shaped distributions centred around zero (Extended Data Fig. 6b). As such, the deviance of lochNESS can be summarized as the average Euclidean distance between lochNESS versus lochNESS under permutation (Extended Data Fig. 6c). In addition, we computed lochNESS between wild types from different background strains and observed minimal variation in cell distribution between the wild type from G4, FVB and BALB/c strains and potential strain-specific distributions in wild-type C57BL/6 mice (Extended Data Fig. 6d).

Comparing lochNESS with the batch-mixing score the local inverse Simpson index

LochNESS shares conceptual similarities with batch-correcting measurement scores such as the local inverse Simpson index (LISI)⁶⁰, which quantifies the amount of mixing in a cell’s neighbourhood by counting the number of batches represented in the neighbourhood. As a direct comparison, we calculated LISI on each mutant with a pooled wild-type reference in PCA space. We calculated LISI with a dynamic perplexity based on the dataset size (${\rm{p}}{\rm{e}}{\rm{r}}{\rm{p}}{\rm{l}}{\rm{e}}{\rm{x}}{\rm{i}}{\rm{t}}{\rm{y}}={\rm{f}}{\rm{l}}{\rm{o}}{\rm{o}}{\rm{r}}(0.5\times \sqrt{(N)/3}$), K = 3 × perplexity), similar to our strategy for determining the neighbourhood size for lochNESS. Focusing on the G4 mutants as an example, the results show a correlation between LISI and lochNESS, for which LISI values close to 1 correspond to the more extreme positive or negative values of lochNESS as expected (Extended Data Fig. 6f). LochNESS has several conceptual advantages compared to LISI. First, lochNESS can easily determine whether the mutant sample is enriched or depleted in an area that is not well mixed using the sign of the value (positive = enrichment, negative = depletion), whereas LISI can separate only mixed (scores approaching 2) versus separated (scores approaching 1). Second, lochNESS can be easily extended to comparisons between multiple samples, whereas LISI is relatively restricted to pairwise comparisons. Third, lochNESS considers a dataset-specific neighbourhood size and baseline proportions.

Identifying lochNESS-associated gene expression changes

To identify gene expression changes associated with mutant-enriched or mutant-depleted areas, we find differentially expressed genes through fitting a regression model for each gene accounting for lochNESS. We use the fit_models() function implemented in Monocle v3 with lochNESS as the model_formula_str. This essentially fits a generalized linear model for each gene: $\log (\,{y}_{i})={\beta }_{0}+{\beta }_{i}\times {x}_{i}$, in which ${y}_{i}$ is the gene expression of gene i, ${\beta }_{n}$ captures the effect of the lochNESS ${x}_{n}$ on expression of gene i, and ${\beta }_{0}$ is the intercept. For each gene i, we test whether ${\beta }_{i}$ is significantly different from zero using a Wald test, and after testing all genes, we adjust the P values using the Benjamini–Hochberg procedure to account for multiple hypothesis testing. We identify the genes that have adjusted P value < 0.05 and large positive ${\beta }_{i}$ values as associated with mutant-enriched areas, and those with large negative ${\beta }_{i}$ values as associated with mutant-depleted areas.

Systematic screening of lochNESS distributions

LochNESS distributions can be systematically screened to identify sub-trajectories exhibiting substantial mutant-specific shifts. For example, although all TAD boundary KO mutants have similarly unremarkable global lochNESS distributions, when we plot these distributions by sub-trajectory, a handful of shifted distributions are evident (Extended Data Fig. 9a). Such deviations, summarized as the average Euclidean distances between lochNESS and lochNESS under permutation, are visualized in Extended Data Fig. 9b. For example, multiple epithelial sub-trajectories, including pre-epidermal keratinocyte, epidermis, branchial arch and lung epithelial trajectories, are most shifted in Tbx3 TAD boundary KO cells. Co-embeddings of mutant and wild-type cells of these sub-trajectories, together with regression analysis, identify multiple keratin genes as positively correlated with lochNESS, consistent with a role for Tbx3 in epidermal development³⁸ (Extended Data Fig. 9c,d and Supplementary Table 7). The lung epithelial cells were separated into two clusters, with the cluster more depleted in Tbx3 TAD boundary KO cells marked by expression of Etv5, which encodes a transcription factor associated with alveolar type II cell development, as well as Bmp signalling genes that regulate Tbx3 during lung development (Bmp1/4), and the distal airway markers Sox9 and Id2 (Supplementary Table 4). Of note, the shifts that we observed in Tbx3 TAD boundary KO cells remain preliminary and would need to be confirmed by further validation experiments.

Spatial mapping with Tangram

We computationally map our dataset onto a spatially resolved transcriptomics dataset, the mouse organogenesis spatiotemporal transcriptomics atlas (MOSTA) generated with Stereo-seq³⁷. The atlas has a total of 53 sagittal sections from C57BL/6 mouse embryos from E9.5 to E16.5 in 1-day intervals, and we obtained one section from the most relevant E13.5 data (E13.5_E1S1.MOSTA.h5ad) from the data-sharing website associated with the manuscript: https://db.cngb.org/stomics/mosta/download/. To map the cells for each single cell cluster on the spatially resolved transcriptomics dataset, we used a machine learning-based method called Tangram⁶¹. Briefly, Tangram is a computational tool that uses a Bayesian approach to infer the spatial locations of cells in a single-cell transcriptomics dataset on the basis of their transcriptomic profiles and the spatial patterns of gene expression in the spatially resolved dataset. The relevant subset of the MMCA data was preprocessed in Scanpy, but the metadata were inherited from the results generated in the section above entitled Cell clustering and annotation. We used Tangram with default parameters to estimate the spatial coordinates of cells from each cluster in the single-cell dataset and visualized results on the coordinates provided by MOSTA. We trained the Tangram model in gpu mode using an NVIDIA A100 GPU. Overall, Tangram provided a powerful method for mapping the cells from the scRNA-seq dataset onto MOSTA, enabling us to infer the spatial locations of different cell clusters of interest within the tissue.

Calculating mutant and embryo similarity scores

We can extend the lochNESS analysis, which is computed on each mutant and its corresponding wild-type mice, to compute ‘similarity scores’ between all pairs of individual embryos from the same background strain. We consider all embryos in the same background in a major trajectory as a dataset. For each dataset, we define a ‘similarity score’ between cell n and embryo j as:

$$\begin{array}{c}{{\rm{s}}{\rm{i}}{\rm{m}}{\rm{i}}{\rm{l}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{t}}{\rm{y}}{\rm{s}}{\rm{c}}{\rm{o}}{\rm{r}}{\rm{e}}}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}n,{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}j}=\frac{{\rm{n}}{\rm{o.}}\,{\rm{o}}{\rm{f}}\,{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}{\rm{s}}\,{\rm{f}}{\rm{r}}{\rm{o}}{\rm{m}}\,{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}\,j\,{\rm{i}}{\rm{n}}\,k\, \mbox{-} \,{\rm{N}}{\rm{N}}{\rm{s}}\,{\rm{o}}{\rm{f}}\,{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}\,n}{k}\\ \,\,\,\,\,\,\,/\frac{{\rm{n}}{\rm{o.}}\,{\rm{o}}{\rm{f}}\,{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}{\rm{s}}\,{\rm{f}}{\rm{r}}{\rm{o}}{\rm{m}}\,{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}\,j\,{\rm{i}}{\rm{n}}\,{\rm{d}}{\rm{a}}{\rm{t}}{\rm{a}}{\rm{s}}{\rm{e}}{\rm{t}}}{N}\end{array}$$

in which N is the total number of cells in the dataset and $k=\frac{\sqrt{N}}{2}$. We take the mean of the similarity scores across all cells in the same embryo, resulting in an embryo similarity score matrix for which entries are:

$${{\rm{s}}{\rm{i}}{\rm{m}}{\rm{i}}{\rm{l}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{t}}{\rm{y}}{\rm{s}}{\rm{c}}{\rm{o}}{\rm{r}}{\rm{e}}}_{{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}i,{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}j}=\frac{1}{{n}_{i}}\mathop{\sum }\limits_{n=1}^{{n}_{i}}{{\rm{s}}{\rm{i}}{\rm{m}}{\rm{i}}{\rm{l}}{\rm{a}}{\rm{r}}{\rm{i}}{\rm{t}}{\rm{y}}{\rm{s}}{\rm{c}}{\rm{o}}{\rm{r}}{\rm{e}}}_{{\rm{c}}{\rm{e}}{\rm{l}}{\rm{l}}n,{\rm{e}}{\rm{m}}{\rm{b}}{\rm{r}}{\rm{y}}{\rm{o}}j}$$

in which ${n}_{i}$ is the number of cells in embryo i.

Identifying and quantifying developmental delay

To identify potential mutant-related developmental delay, we integrate MMCA with MOCA. We consider a mutant and its corresponding wild type in a sub-trajectory as a dataset. We take the cells from E11.5 to E13.5 with similar annotations from MOCA and co-embed with the MMCA cells. We take the raw counts from both datasets, normalize and process the data together without explicit batch correction as both datasets were generated with sci-RNA-seq3 and were similar in dataset quality. We visualize the co-embedded data in three-dimensional UMAP space and check for developmental delay in the mutant cells (that is, mutant cells embedded closer to early MOCA cells compared to wild-type cells). To quantify the amount of developmental delay, we find k-NNs in MOCA for each cell in MMCA and calculate ${\rm{time\; score}}=\frac{{\sum }_{n=1}^{k}{T}_{n}}{k}$, in which ${T}_{n}$ is the developmental time of MOCA cell n in the k-NNs of the MMCA cell. Afterwards, we test whether the average time scores of mutant cells are significantly different from that of wild-type cells using a Student’s t-test.

RNAscope in situ hybridization

For RNAscope, embryos were collected at stage E13.5 and fixed for 4 h in 4% PFA in PBS at room temperature. The embryos were washed twice in PBS before incubation in a sucrose series (5%, 10% and finally 15% sucrose (Roth) in PBS) each for 1 h or until the embryos sank to the bottom of the tube. Finally, the embryos were incubated in 15% sucrose in PBS and O.C.T. (Sakura) in a 1:1 solution before embedding the embryos in O.C.T. in a chilled ethanol bath and storing them at −80 °C until sectioning. The embryos were cut into 5-μm-thick sections on slides for RNAscope.

Simultaneous RNA in situ hybridization was carried out using the RNAscope technology (Advanced Cell Diagnostics (ACD)) and the following probes specific for Mm-K (catalogue number 476261, ACD) and Mm-Sox9-C2 (catalogue number 401051-C2, ACD) on 5-μm sections of the mouse embryos. RNAscope probes were purchased from ACD and designed as described⁶². The RNAscope assay was run on a HybEZ II Hybridization System (catalogue number 321720, ACD) using the RNAscope Multiplex Fluorescent Reagent Kit v2 (catalogue number 323100, ACD) and the manufacturer’s protocol for fixed-frozen tissue samples with target retrieval on a hotplate for 5 min. Fluorescent labelling of the RNAscope probes was achieved by using OPAL 520 and OPAL 570 dyes (catalogue numbers FP1487001KT and FP1488001KT, Akoya Biosciences), and stained sections were scanned at ×25 magnification using an LSM 980 with Airyscan 2 (Zeiss).

Image analysis

For quantitative analysis of the RNAscope images, representative fields of view for each stained section were analysed using the image processing software Fiji⁶³. The mRNA signal for each organ of interest was counted in a defined area (1 × 1 mm²), with n = 6 per condition. Statistics were calculated using Student’ t-test and evaluated (not significant, P > 0.05; *P < 0.05 to ≥ 0.01; **P < 0.01 to ≥ 0.001; ***P < 0.001).

Ttc21b- and Gli2-mutant fixation for haematoxylin and eosin staining and immunofluorescence

Homozygous and heterozygous Ttc21b mutants and wild-type E13.5 mouse embryos were fixed overnight in 4% PFA at 4 °C. To stop fixation, the samples were transferred into 70% ethanol, washed twice and dehydrated. In the following, the embryos were embedded in paraffin, and cut into 2.5-μm-thick sections.

Ttc21b-mutant haematoxylin and eosin staining

Histochemical staining was carried out on the eyes of the embryos using haematoxylin and eosin. Slides were scanned with a digital slide scanner (NanoZoomer 2.0HT) and analysed using NDP.view2 software (Hamamatsu Photonics). The following numbers of embryos were processed: 2 wild type; 2 heterozygous Ttc21b ; 4 homozygous Ttc21b.

Gli2-mutant haematoxylin and eosin staining and immunofluorescence

For the histological analysis, haematoxylin and eosin staining of E13.5 Gli2-KO mouse embryos, and respective wild-type littermates (n = 4 and n = 2, respectively), was carried out on 4% paraformaldehyde-fixed paraffin-embedded sections (3 μm). Stained paraffin sections were scanned using a digital slide scanner (NanoZoomer 2.0HT) and examined using NDP.view2 software. The cut regions and positions were annotated according to ref. ⁶⁴.

The spatial abundance patterns of prealbumin as a marker for ChP and PAX6 as a marker for neural tube development were analysed by immunofluorescence, using specific antibodies (rabbit monoclonal (EPR20971) to prealbumin (1:1,000, Abcam) and rabbit polyclonal antibody to PAX6 (1:200, AB2237 Merck Sigma) in an automated BOND Research Detection system. Antibody binding was detected by goat anti-rabbit Alexa Fluor 488-conjugated secondary antibody (Leica, A-11008). Nuclear counterstaining was achieved using 4′,6-diamino-2-phenylindole (DAPI). In negative-control sections, the primary antibodies were omitted and antibody diluent was applied.

Stained embryo sections were scanned with an AxioScan 7 digital slide scanner (Zeiss).

Fluorescence quantification

Quantification of prealbumin expression cells was carried out using the image analysis software Definiens Developer XD2 (Definiens). The regions of interest (1–4) within the fourth and lateral ventricle ChP were annotated manually in serial sections. The calculated parameter was the ratio of the total number of prealbumin-positive cells over the embryo section area (in micrometres).

Statistics and reproducibility

Haematoxylin and eosin staining of the developing eye (Fig. 2c) was carried out on homozygous Ttc21b mutants (n = 4), heterozygous Ttc21b mutants (n = 2) and wild-type E13.5 embryos (n = 2). Experiments on the sections were carried out in parallel to ensure consistency.

Haematoxylin and eosin staining of Gli2-mutant and wild-type embryo sections (Extended Data Fig. 8a) was carried out on homozygous Gli2-KO (n = 4) and wild-type (n = 2) samples. Experiments on the sections were carried out in parallel to ensure consistency.

Immunofluorescence staining of the ChP marker TTR and neural tube marker PAX6 (Fig. 3f and Extended Data Fig. 8b–d) was carried out on sections of homozygous Gli2-KO (n = 4) and wild-type (n = 2) samples. Immunofluorescence of the same antibody was carried out on all mutants in parallel to ensure consistency.

Sox9 and Kcnj2 expression of heterozygous E13.5 wild-type and Sox9 regulatory INV mutant embryos (n = 6 embryos for each condition) was measured by RNAscope image quantification in a defined area (1 ×1 mm²). Statistics were calculated using a two-sided Student’ t-test and evaluated as follows: not significant, P > 0.05; *P < 0.05 to ≥ 0.01; **P < 0.01 to ≥ 0.001; ***P < 0.001. RNAscope of the tissue was carried out on all samples in parallel to ensure consistency.

Clustering and annotation limb mesenchyme trajectory

Seurat v4.0.6 was used for the analysis. Wild-type cells in the limb mesenchyme trajectory from all wild-type mice (n = 15 mice, n = 25,211 cells) were used to first annotate the cells. The raw counts were log-normalized, after which PCA was carried out with default parameters on the top 2,000 highly variable genes selected using the vst method. Nearest neighbours were computed on the PCA space, with default parameters, except that all of the PCs computed earlier were used. Clustering was carried out using the Louvain community detection algorithm with a resolution of 0.1, resulting in three clusters. Positive marker genes for these clusters were identified using the Wilcoxon rank-sum test, for which only the genes expressed in at least 20% of the cells in either cell group were considered. The clusters were annotated on the basis of biologically relevant markers (Extended Data Fig. 12f). The newly assigned cell annotations for the limb mesenchyme trajectory cells in the wild-type dataset were transferred to the corresponding cells in the Sox9 regulatory INV mutant using the FindTransferAnchors and TransferData functions using default parameters, except that all of the computed PCs were used. A total of 92.3% of the transferred annotations had a score (prediction.score.max) greater than or equal to 0.8.

Density visualization and RNA velocity analysis

Using Seurat v4.0.6, the raw counts were log-normalized, and PCA was carried out with default parameters on the top highly variable genes 2,000 genes, selected using the vst method. Dimensionality reduction was carried out using PCA with default parameters, after which the UMAP embedding was carried out on all computed PC components. Density plots were created using the stat_2d_density_filled function in ggplot2 v3.3.5. For RNA velocity analysis using scVelo v0.2.4, the total, spliced and unspliced count matrices, along with the UMAP embeddings, were exported as an h5ad file using anndata v0.7.5.2 for R. The count matrices were filtered and normalized using scv.pp.filter_and_normalize, with min_shared_counts = 20 and n_top_genes = 2,000. Means and variances between 30 nearest neighbours were calculated in the PCA space (n_pcs = 50, to be consistent with the default value in Seurat). The velocities were calculated using default parameters and projected onto the UMAP embedding exported from Seurat.

Single-sample gene set enrichment analysis

Single-sample gene set enrichment analysis was applied to scRNA-seq data using the escape package in R⁶⁵. The msigdbr and getGeneSets functions were used to fetch and filter the entire hallmark (H; 50 sets) or the signature cell type (C8; 700 sets) M. musculus gene sets from MSigDB⁶⁶. enrichIt with default parameters, except for using 10,000 groups and variable number of cores, was carried out on the Seurat object containing data corresponding to the undifferentiated mesenchyme cells from the Sox9 regulatory INV mutant, after converting the feature names to gene symbols as necessitated by the escape package. The obtained enrichment scores for each gene set were compared between the two branches (Fig. 5f) using the two-sample Wilcoxon test (wilcox_test) with default parameters and adjusted for multiple comparisons using Bonferroni correction.

Integration and spatial mapping with sci-space data

We integrated our dataset with a spatial transcriptomics dataset on mid-gestational mice (E14.5), based on the sci-space method⁴⁷, in which a subset of transcriptionally profiled nuclei have known physical locations in sagittal sections within which they were mapped before scRNA-seq. We used anchor-based integration as implemented by Seurat for a co-embedding of a subset of MMCA and sci-space. For cells in the subset of MMCA, we find the nearest neighbour in sci-space data in the integrated co-embedding, and plot the location of the neighbouring sci-space cell if it is known.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.