Introduction

Embryogenesis involves not only the creation but, crucially, the organization of new cells and new cell identities. 3D gastruloids model an important stage of embryogenesis: gastrulation. Gastruloids break symmetry, elongate along an anterior-posterior (AP) axis, and form many of the cell types found in embryos at a similar developmental stage, as shown by single-cell RNA-seq analyses [1–6]. These characteristics mean they could hold enormous potential to enable us to deeply characterize, perturb, and test developmental hypotheses, circumventing many practical and ethical constraints imposed by using in vivo models. Yet despite extensive single-cell characterization of gastruloids, we still lack detailed information about the arrangement of cell types within individual gastruloids.

Sequencing-based gene expression analyses of gastruloids have at most been able to generate coarse-grained maps of the spatial arrangement of the constituent cells via sectioning [2,7]. These methods cannot resolve the position of individual cells, nor can they measure how tissues are organized other than along the AP axis. Computational methods that measure gene expression and identity of individual cells can only estimate AP axis location through inference methods and do not measure organization in any other dimension [8,9]. Additionally, detailed comparative spatial analyses of individual gastruloids have not been performed, as the experiments either require pooling gastruloids or are low-throughput, such that there are insufficient samples to make meaningful comparisons. While some studies point to stringent reproducibility of gene expression along the AP axis in gastruloids [10,11], others show that there is great variation, gastruloid-to-gastruloid, in the organization of specific tissues [12]. Without high-resolution spatial characterization across multiple gastruloids, it’s unclear what aspects of tissue development can be effectively modeled with this system.

The cell types present in gastruloids have largely been inferred from single-cell sequencing analyses. As is standard for these experiments, cells are clustered on the basis of their overall gene expression profiles, and types are assigned based on a qualitative assessment of the genes differentially expressed by the cluster [2,7,13–15]. However, it is unclear whether discrete binning cells into categories is the only, or even the most useful, approach to single-cell studies of gastruloids, a system nearly entirely composed of cells in the process of differentiating and transitioning between identities. A more gene-centric analysis, especially one coupled to a spatial measurement, could more naturally represent these transitions and how they are coupled to the arrangement and physical location of cells. Whether the diversity of cell types observed in sequencing experiments truly maps onto spatially organized collections of cells, and if this organization is consistent between gastruloids, remains an open question.

To address these gaps, we developed a spatially resolved, single-cell molecular map of the location, identity, and gene expression of cells within 26 individual gastruloids. We found that despite some morphological variability, elongated, polarized gastruloids had largely reproducible cell type composition. The order of types and expression of individual genes along the AP axis was also quite consistent. Leveraging our spatial single-cell expression data, we found that the posterior of the gastruloids was characteristically patterned with distinct clusters of progenitors and differentiated cell types. We present novel computational techniques useful for identifying genes with mutually exclusive expression patterns, which illuminated the spatial component of a bipotent progenitor differentiation decision and revealed a spatial hierarchy of gene expression patterns. Using evidence generated from this metric, we show evidence of spatial and functional control of endothelial cells, suggesting that 3D gastruloids may be used in the future to model early stages of vasculogenesis and organization.

Results

Cell types’ locations and relative proportions are consistent across morphologically normal gastruloids

To measure the spatial distribution of gene expression, we prepared gastruloids using mouse E14TG2a cells and a standard protocol (see Methods). We harvested mature gastruloids after 120 hours of growth and assayed with seqFISH [16]. There was substantial variation in the length, width, and relative amounts of anterior and posterior tissues in the final dataset, although all were within the range of what would be considered a ‘morphologically normal’ gastruloid [1,10]. We assigned an AP axis to each gastruloid using the expression of T, a canonical marker for the posterior (Figure 1a). When we compared how gene expression varies along the AP axis, we saw good agreement at a coarse-grained level with a previous study that sectioned gastruloids along the axis and analyzed gene expression in each section [2] (Figure S1.3a).

Figure 1:

We assigned a cell type to each nucleus using a cell type scoring method with known marker genes. Unsupervised clustering produced similar results (Figure S4.4a), but we achieved greater cell type distinction when we scored using a marker gene panel; see Methods for additional details. A representative gallery of typed gastruloids is shown in Figure 1b; the full dataset is in Figure S1.2. Because each cell received a cell type score for each type, we could use the entropy of the cell type score probability distribution to assess confidence of our cell type assignment: a cell that received a similar score for multiple cell types would have a high entropy distribution, while one which scored highly for one type and low for the rest would have low entropy. On average, the cell type entropies for most cells in a given cell type were low, with the exception of paraxial mesoderm and cardiac mesoderm cells, which had intermediate values (see additional discussion below). The generally low entropies indicate that most of our cell type assignments were high-confidence (Figure S1.1a).

Once we had the cell type identity and spatial location of each cell in all the gastruloids, we characterized the organization of each by mapping where each cell type was found relative to other types and overall morphology. The approximate location of each cell type relative to gastruloid morphology was quite consistent: the posterior region, although highly variable in size (Figure S1.2a,b,c), mainly consisted of neuromesodermal precursors (NMP, turquoise), a bipotent cell type that contributes to both neural and mesodermal tissues [17–19]. This localization is consistent with previous observations of the expression of NMP genes in the gastruloid posterior [2,7,20]. However, unlike these previous studies, which inferred NMPs from the expression of single genes and could not resolve other cell types, we were able to see that NMPs were intermixed in a patchy fashion with mesodermally-fated cells that are on the path to forming somitic tissues (presomitic mesoderm, PSM, red). In some instances, we saw PSM cells forming symmetric patches on either side of a central neurally-fated region (Figure 1d, center gastruloid, Figure S1.2a v., xi., xii, b iv), mimicking the organization found in embryos [21–23] and thought to occur in gastruloids [8,24]. Generally, we saw 1-2 distinct clusters of neurally fated cells (spinal cord precursors, black), often adjacent to or opposite the region associated with somitic fate acquisition (“differentiation front”, light blue), although occasionally found in the center of the gastruloid (Figure S1.2a ii., vi., xvii). Although the distinct organization of neurally fated cells has been observed in other embryo models [3], the organization of spinal cord precursors in gastruloids in the context of other cell types has not yet been systematically characterized; staining for spinal cord-associated genes showed compelling patterning, but the patterns varied by gene, complicating interpretation of the organization of the tissue [4]. Consistent with [2], the differentiation front was nearly always found in one distinct cluster marking the boundary between the anterior and posterior (Figure S1.2a i., iii., vi., b iii., vi). The anterior regions, although clearly not organized into somites, consisted mainly of somitic cells, presomitic cells, and cardiac mesoderm (yellow, green, and purple), with highly variable mixing (Figure S1.2 a iii., v., x., xv., xvi) . Also present in the posterior of almost all the gastruloids were small clusters of endothelial precursors (pink) (Figure S1.2 a iii., vi., vii., xvii., c), in some cases associated with a cluster of endodermal cells (blue) (Figure S1.2 a x., xi., xiv., xv). In general, if an endodermal region was present, it was located in the center of the gastruloid (Figure S1.2 a v., viii., xii., xvii., b ii., iv., c), consistent with observations in [12]. Although there was clearly variation in the details of cell type arrangement, there were obvious organizational principles that were conserved between gastruloids independent of size, consistent with the fidelity of gene expression along the AP axis reported in reference [10]. This simultaneous mapping of all cell types is consistent with many observations in the literature about the expression of individual genes and integrates them into one consistent picture of cell type organization.

We sought to quantify variability in cell type composition between gastruloids. Early single-cell datasets relied on pooling multiple gastruloids, thus obscuring the degree to which the overall cell type distribution was reflected in each individual gastruloid. However, recent single-cell measurements of individual gastruloids have suggested substantial gastruloid-to-gastruloid variation in cell type proportions [13]. Figure 1c shows distributions of cell type proportions across samples. Figure S1.1b shows the per-gastruloid distributions. In general, we found that cell type proportions did not vary significantly between gastruloids, with the exception of endoderm, paraxial mesoderm, and somite. Unlike other cell types that were consistently present but varied in proportion, endoderm variation largely resulted from the absence of this tissue in some gastruloids, reflecting the developmental variability of endodermal tissues previously documented in gastruloids [12]. The variation in the proportion of somite cells could be at least partially explained by variability in the amount of anterior tissue relative to posterior tissue (Figure 1d, left), and given that somitic tissue forms progressively over the course of gastruloid development, could potentially be explained by desynchronization of somite formation such that some gastruloids had, at 120 hours, formed more of this cell type which emerges later in gastruloid development [13]. We did not observe gastruloids that were as strongly neurally-biased as those reported in [13], but we did see some gastruloids with a relatively high proportion of spinal cord precursor cells (Figure S1.2a ii., xv., b vii), and overall the proportion of spinal cord had a negative covariation with the mesodermally-derived cell types, consistent with the anticorrelation also reported in [13] (Figure S1.1c).

The proportion of somite cells was strongly positively correlated with the proportion of presomitic mesoderm cells (covariation = 0.63, Figure S1.1c). Presomitic mesoderm cells progressively take on somitic fate at and after the differentiation front; a positive correlation between these two would suggest that the amount of presomitic mesoderm is linked to how much somite is ultimately formed. How much presomitic mesoderm is formed depends on fate decisions made by NMPs: NMP progeny can become either mesodermal or neural, depending on their local signaling environment. This chain of correlation suggests that fate decisions made by NMPs may drive the overall composition of individual gastruloids and that, as in embryos, a balance between the amount of somitic and mesodermal tissue is achieved, although likely to a less stringent degree [25]. Further experiments could delineate whether, as in embryos, Wnt3a drives this balance [25], and the degree to which growth of differentiated progeny affects compositional variation [20].

Finally, paraxial mesoderm variation could only be weakly explained by morphological characteristics (Figure 1d, right), and couldn’t be explained by covariation with any other cell type (Figure S1.1e). Given that paraxial mesoderm also has a relatively higher average cell type entropy score, it suggests either that our marker genes for this cell type were not sufficient to confidently assign a score, that their gene expression is more variable than other types, or that this type itself is transient in gastruloids, as has been suggested in the literature [26].

Cell types’ arrangement is consistent across gastruloids, but varies by type

Given that the proportions of cell types within each gastruloid were fairly consistent, to what degree did their spatial organization vary? We first projected each cell’s expression onto the AP axis and looked at the distribution of AP axis locations across gastruloids (Figure 1e). The most posterior cell types (NMP and presomitic mesoderm) showed consistent, wide distributions from 0-30% of the axis. Centered around 30%, spinal cord precursors and endoderm had distinctive peaks (clusters), the exact location of which varied between gastruloids. The differentiation front was consistently located in one peak between 30-40% of the AP axis. The anterior cell types tended to be more broadly distributed without distinct peaks; however, endothelial precursors in particular had wide, multimodal distributions that varied substantially between gastruloids despite a comparably small proportional variation (Figure 1c). Generally speaking, variation in the proportion of a cell type across gastruloids was not directly predictive of the variation in location of that cell type relative to the AP axis. From these observations, we concluded that the amount of any given cell type is not the sole determining factor in where it is located along the AP axis.

Several studies of gene expression in gastruloids have used pooled measurements to infer the AP axis-location of genes and cell types [3,7,13,24,27] and our data are largely consistent with these lower-resolution findings (Figure S1.3a). Yet it is obvious from individual gene staining [1,4,8,28] and our detailed 2D maps of cell identity and location that gastruloid organization is much more complex than the average order of cells along the AP axis. We thus needed an analytical method for quantifying spatial organization beyond distributions along the AP axis. To characterize spatial organization, we quantified the overall mixing of cell types within each gastruloid using a mixing coefficient that quantifies how often each cell is found next to a cell of a different type or the same type. These neighbor relationships are summarized into a single score that varies from 0 (completely mixed) to 1 (completely segregated) (Figure 1f). The mixing coefficient values range from 0.29 to 0.58 (Figure 1g). With the exception of cardiac mesoderm, the amount of any given cell type did not significantly correlate with the mixing coefficient (Figure S1.1d). The proportion of cardiac mesoderm negatively correlates with mixing, implying that this cell type is particularly likely to be found as a singlet.

Our spatial transcriptomics is imaging-based, and the fact that we image transcripts in a single plane admits the possibility that, in any individual gastruloid, we may collect data from a different part of the gastruloid. We controlled for this to the extent possible within experimental limitations by imaging multiple gastruloids across several experiments and keeping our imaging parameters, particularly the instrument z-depth relative to the coverslip, nearly identical across experiments. The relatively wide distribution of mixing coefficients demonstrates that even within gastruloids with broadly similar morphologies and cell type proportions, the underlying organization of cell types can vary substantially.

Pairwise cell type interactions suggest more mixing of cell types in the gastruloid anterior than posterior

Our simultaneous observation of variability in mixing coefficients and consistent cell type organization along the AP axis led us to ask whether mixing coefficients were driven largely by local topology and whether this effect was disproportionately due to the arrangement of a subset of cell types. The interaction between long-range signaling gradients and short-range cell-cell interactions is thought to drive gastruloid symmetry breaking and patterning [29–31], and our spatial data allowed us to query the cell-cell level organization and interaction between cells. First, we asked to what degree each cell type is exposed to all other cell types (including those that could not be typed). We calculated an “exposure index”, which quantifies the degree to which each cell type (source) is exposed to another cell type (neighbors). It is similar to the mixing coefficient used to score the overall organization of a gastruloid, but differs in that the mixing coefficient only measures differences in self vs. non-self interactions between cell types, whereas the exposure index tracks interactions between all pairs of cell types. The measure is asymmetric in that the exposure of cell type A to B may not be the same as the exposure of cell type B to A; for example, a cell type found in singlets within a larger tissue will be more frequently exposed to cells of that common tissue type, while the common type will see itself much more frequently than the intermixed cell. Thus, this asymmetry accounts for the fact that local patterning can influence the degree to which any one cell type is found next to another. Asymmetry in off-diagonal entries in the exposure matrix (Figure 2a), therefore, can reflect local patterning. For instance, somite and endothelial precursors have an asymmetric pattern: although endothelial cells are often exposed to somite cells (exposure index 0.207), somite cells are comparatively infrequently exposed to endothelial cells (0.07). This observation, coupled with a fairly strong endothelial-endothelial exposure index (Figure S2a), suggests small clusters of endothelial cells are found among somite cells, reflecting our visual impressions. In contrast, cardiac mesoderm cells were often exposed to all other anterior cell types, but their self-exposure index was low, matching our observation that they are found in inconsistent and disorganized clumps or individually. Like paraxial mesoderm, this cell type also has a high average cell type score entropy (Figure S1a).

Figure 2:

The overall interaction matrix (Figure 2a) revealed that there were more inter-type exposures between the anterior cell types of the gastruloid than between the posterior types. This quantification matches the impression from visual inspection that the gastruloids contain both somitic cells and paraxial mesoderm cells that properly localize to the anterior (as per [2,8]) but are arranged in a disorganized fashion. Cells of the posterior, by contrast, were much less mixed, with individual cell types forming large, distinct clusters. These results demonstrate that our analytical methods were able to capture and quantify the degree of spatial disorder in these two compartments. The contrast between organized posterior clustering and disorganized anterior mixing matches expectations based on literature that shows that self-organization mechanisms in gastruloids in the anterior vs. posterior are differentially sensitive to culture conditions, with somitic patterning requiring external matrix support [1,4,8,24], distinguishing it from the seemingly more autonomous organization observed in posterior cell types.

Local topology of cell type interactions is preserved across gastruloids

Our analyses have been focused on pairwise interactions, but these data also have the potential to reveal higher-order interactions between cell types. To measure such interactions, we first determined all the possible combinations of three cells (triplet motif) from the 9 cell types we had identified (with replacement). Then, for each cell in the gastruloid, we identified the two nearest neighbors and determined which triplet type it was. We calculated the enrichment of each triplet motif relative to its expected frequency based on the amount of each of the constituent cell types. Figure 2b shows this distribution of triplet motifs that were found in over half the gastruloids and in each were significantly enriched or depleted based on overall cell type frequency. Figure S2b shows the same distribution but for triplets that were found in any gastruloid. We found that motifs formed from 3 of the same cell type were overwhelmingly the most significantly enriched, and all were found in over half of the gastruloids. The enriched 2+1 (two same) motifs were composed of all posterior cell types or all anterior cell types. 3/4 of the most common significantly enriched motifs that were of three different cell types consisted of cardiac mesoderm with two other anterior cell types or differentiation front, again suggesting that the cardiac mesoderm cells are uniquely spread out as individuals.

We also identified motifs that were highly underrepresented compared to what would be expected given the frequency of occurrence of each of the individual cell types: nearly all of the significantly underrepresented motifs were of 3 different cell types, with the rarest being somite+NMP+presomitic mesoderm. This observation was initially surprising given that these cells are all thought (in gastruloids) to derive from the same lineage; however, the early mesodermal precursors are spatially segregated from the anterior by the differentiation front, which is where somite determination happens [2,3,24] The relative paucity of this triplet indicates that the spatial organization of the differentiation process occurs with very high fidelity, and that less differentiated precursors rarely make it into the anterior, where they would be exposed to somite and paraxial mesoderm cells.

Overall, these data speak to the consistency of gastruloid organization, although broad-scale morphological differences can explain some differences in cell type proportion and rare individual gastruloids may have unique cell type arrangements. In general, morphologically normal gastruloids consistently patterned posterior cell types and moved more differentiated cell types to the anterior. Along the way, the granular interactions between individual cells were reproduced with remarkable fidelity. It is notable that this arrangement across length scales was achieved with only minimal external chemical cues.

Co-expression of cell type marker genes delineates progressive NMP differentiation in gastruloids

One of the central goals of studying developmental processes is to determine the means by which cells make fate decisions and when these decisions are driven by the physical location of the cell as opposed to cell-intrinsic programs. A large part of gastruloid development is thought to be driven by the differentiation, growth, and migration of NMPs and their progeny, which are fated towards either a mesodermal or posterior neural (spinal cord) fate [3,28,32]. In embryos, this process is thought to happen in a continuous fashion; NMP progeny simultaneously move anteriorly and acquire either a mesodermal or spinal cord fate [18,25,33,34]. The differently fated cells are driven to different parts of the developing spinal column; lineage tracing evidence suggests that fate commitment is at least partially driven by physical location and the local signaling environment [25], and is progressive in nature [33]. We wondered how this coupling between space and differentiation was resolved in gastruloids, and whether the continuous transition from progenitor to more differentiated cell type was reflected in spatially continuous patterns of gene expression.

Our data represent a single slice of time during which NMPs, presomitic mesoderm, and spinal cord precursors exist simultaneously; if this process happened continuously, we would also expect to see a subset of cells somewhere on an identity continuum between these types (Figure 3a). For simplicity, when initially assigning cell types, we used differential expression of panels of marker genes to assign identity. However, many of the cells in the posterior had non-zero cell type scores in all three categories (Figure S3.1a). This observation motivated us to leverage the advantages of our dataset to delineate where individual genes were expressed in space with single-cell resolution and determine whether that additional level of resolution gave additional information about NMP differentiation trajectories.

Figure 3:

Since both presomitic mesoderm and spinal cord cells derive from NMPs, we reasoned that it was possible, especially if the cells were still in the process of differentiating, for cells typed as presomitic mesoderm or spinal cord to residually express NMP genes. Indeed, there were genes in our panel that were shared between each type and the precursor (referred to here as NMP+presomitic mesoderm genes and NMP+spinal cord genes or “shared genes”). Figure 3c shows the total expression of each gene group versus the NMP genes for all cells in all gastruloids that we typed as NMP, presomitic mesoderm, or spinal cord. As expected, there was a clear correlation between the mixed categories and NMP genes, supporting the notion of a continuous differentiation process.

In an effort to locate cells that had fully committed to one lineage or the other, we looked at the expression of genes thought to be exclusively associated with presomitic mesoderm or spinal cord. Our expectation was that the expression of these genes would be mutually exclusive, and several genes did display this pattern (Figure 3b, right). However, we were surprised to find many examples where these supposed marker genes for different cell types were expressed in the same cell (Figure 3b, left). Spot assignment in imaging data has the potential to mis-assign transcripts, which could, in principle, lead to artifactual mixed states. However, we don’t believe misassignment explained our observations. We used stringent nuclear dilation (small enough that ∼30% of transcripts were not assigned to a nucleus). Moreover, when we calculated the exposure index, we found very little exposure between presomitic mesoderm and spinal cord cells (0.07 and 0.04, Figure 2a), meaning the majority of the time these cells were not nearest neighbors, and thus had little potential for overlap and transcript mis-assignment.

Was this trend of co-expression true for the two groups of marker genes as a whole? When we compared the total per-cell expression of the genes in both categories (spinal cord and presomitic mesoderm), we noted that while they did show more mutually exclusive (i.e. L-shaped) expression than when compared to the NMP genes, there were clearly many cells that had substantial counts from both PSM and spinal cord gene categories (Figure 3d). One possible explanation for the co-expression of genes associated with both arms of the differentiation trajectory is that these cells had not yet fully committed to either lineage. It is also possible that a subset of the genes were better (or at least more exclusive) markers than others (as our observation of the heterogeneity in expression patterns suggested), and that there may in fact be a continuum of markers, some of which may distinguish certain cell types and not others. To determine if any of these possibilities were true, we needed a way to identify genes that were expressed in a mutually exclusive manner.

The L-metric quantifies mutually exclusive gene expression

Marker genes for distinct lineages are often assumed to be mutually exclusive in their expression. However, the above results show that such assumptions may not hold. We were curious whether there were genes that were truly exclusively associated with each cell type, and if such genes in general would prove to be more effective marker genes than those that were merely associated with the cell type.

To this end, we developed a pairwise metric between genes that reported the degree of mutually exclusive expression. It is calculated by rank ordering cells by the expression of one gene and measuring the degree to which the expression of the other gene is anti-rank-ordered (see Methods for details). We call this metric the “single-cell L-metric” (scL-metric) because when the per-cell expression of mutually exclusive genes was plotted against one another, the data made an L shape (Figure 3e, right). A value of -1 represents perfectly mutually exclusive expression, which only happens when the genes are never found in the same cell. Higher values indicate more co-expression. The maximum possible value is 1, which is obtained when both genes express in exactly the same cells, independent of the total amount. The scL-metric is inherently asymmetric: if the two pair orderings produce very different L values, this indicates that there is a corresponding asymmetry in their expression. For example, all counts of the first gene may be found in cells expressing the second gene, yielding an scL-metric closer to 1. However, the second gene may be expressed in many cells with no expression of the first gene, and, thus, swapping the pair ordering results in an scL-metric value closer to 0. In most cases, the values were similar, so we averaged them for simplicity, thereby symmetrizing the metric. The scL-metric is parameter- and threshold-free, which is an advantage over methods that require thresholds to binarize data before calculating an odds ratio or mutual information.

Our expectation is that ubiquitously expressed genes like cell cycle and housekeeping genes will have relatively high scL-metric values no matter which genes they are compared with, whereas genes that are specifically associated with a single cell type will have an scL-metric value close to -1 when compared with genes specific to other types, but higher values when compared with genes associated with the same cell type. An example with two NMP genes (Cdx2 and Cdx4) with a high degree of co-expression (scL-metric=0.41) is in Figure S3.1e. In the case of NMP differentiation, we could use this metric to determine whether there are pairs of genes that definitively delineate whether a cell has spinal cord or mesodermal identity, and furthermore, if there are genes that distinguish these more differentiated cells from the parental cell type.

We compared the distribution of L-metric values between spinal cord and presomitic mesoderm genes; as our previous results suggested, there were many genes that had very strong mutually exclusive expression, while others were expressed in the same cells. For example, Pax6 (spinal cord) and Eogt (presomitic mesoderm) are often expressed in the same cells; they have a mild positive correlation as measured by Pearson’s R (0.24), and an scL-metric of -0.55 (Figure 3e, left). However, Rfx4 (spinal cord) and Eogt (presomitic mesoderm) were expressed mostly in distinct cells; their expression was nearly uncorrelated when measured by Pearson’s R (-0.02), but they have an scL-metric of -0.96 (Figure 3e, right). Thus, Rfx4 appears to be a good marker for differentiated spinal cord cells, while Pax6 might mark cells that are neurally fated but not yet committed. We also note that it is not possible to distinguish from these analyses whether all differentiated spinal cord cells express Rfx4, merely that it is exclusively associated with this type.

Overall, the distribution of spinal cord and presomitic mesoderm gene sc-L-metric values was strongly skewed towards -1, with at least 20 pairs of genes strongly distinguishing between the two types (Figure 3f, left). From the 54 possible pairs of spinal cord and presomitic mesoderm genes, we found that 24 (44.4%) had scL-metric values < -0.9, likely indicating that they mark different populations of cells.

We predicted that there would be more cell-wise co-expression (i.e., higher scL-metric) between NMP genes and the other two categories (presomitic mesoderm and spinal cord), since these two types of cells derive from NMPs. When we compared the distribution of scL-metric values of NMP and spinal cord or NMP and presomitic mesoderm gene pairs, we did in fact see a shift of the distribution to higher values compared to the presomitic mesoderm and spinal cord pairs, indicating more co-expression (Figure 3f, right).

Despite this shift to higher values, some pairs still seemed to delineate the cell types. For example, Fgf9 (presomitic mesoderm) and Wnt3a (NMP) have a very low (L = -0.84). There were even more low-L pairs when comparing NMP and spinal cord genes, such as Sox21 and Wnt3a (L=-0.98) and Scube2 (spinal cord) and Cdx1 (NMP) (L=-0.91). Notably, the average scL-metric of presomitic mesoderm and NMP gene pairs was much higher than spinal cord and NMP gene pairs, which suggests that perhaps commitment to a neural pathway involves a more systematic reorganization of gene expression.

Analysis using this metric allowed us to distinguish genes such as Eogt and Rfx4 that clearly mark distinct subsets of cells, from those like Eogt and Pax6 that partially overlap and may mark intermediates as well as more differentiated cells. We also observed a unique pattern when we compared the expression of Rfx4, which we identified as a good marker of spinal cord genes, with Nkx1-2, a shared gene annotated as being associated with both NMPs and spinal cord cells. The expression of Nkx1-2 was strongest in the posterior, but extended into the region marked by Rfx4. It seemed to be progressively downregulated in this region from posterior to anterior, producing a graded expression pattern with the anterior-most cells only expressing Rfx4 (Figure S3.1c, d). This observation suggests that shared genes can mark cells on their way to becoming spinal cord cells, and at least some are turned off as the cells differentiate, suggesting that the continuous differentiation pathway observed in embryos is also found in gastruloids.

These graded expression patterns led us to be more broadly interested in the patterns of spatial expression and scL-metric values of genes that were annotated as belonging to more than one cell type. Up until now, we had excluded such genes from our analyses, but we reasoned that their distribution and co-expression might reveal more about the path that NMPs take on their way to becoming presomitic mesoderm or spinal cord cells. Shared NMP+presomitic mesoderm genes included T, Fgf17, and Hoxaas3, while shared NMP+spinal cord genes included Sox2, Gbx2, and Nkx1-2. When we compared both shared categories to their terminal cell type, we found on average that the distributions of L-metric values were positively skewed compared to the NMP and terminal cell type distributions, indicating that, as expected, there was even more co-expression among these categories (Figure S3.1b). Interestingly, we found that both distributions were bimodal: a subset of genes were clustered around 0, indicating high co-expression with the terminal cell type genes (this included Dll1 (presomitic mesoderm) and T (NMP+presomitic mesoderm), and Cdh6 (spinal cord) and Nkx1-2 (spinal cord+NMP). Others were closer to -1, meaning the gene in the shared category is either more associated with the parental type (NMPs) or potentially is specific to intermediates. Examples of these pairs include Mesp2 (presomitic mesoderm) and Msgn1 (presomitic mesoderm+NMP), and Rfx4 (spinal cord) and Hoxb6 (spinal cord+NMP).

We wondered whether we could use the scL-metric scores to refine our marker gene panel. We reasoned that ‘good’ marker genes would have very low scL-metric scores with most other genes and have high scL-metric values with a small subset of genes (presumably those belonging to the same type). However, filtering on these criteria did not markedly improve our cell type scoring, likely because the genes in our panel were hand-selected to be good type representatives, so most fit this criterion (Figure S4.3a, b and Methods). However, we found that scL-metric analysis could be used to quantify how related cell types were in terms of marker gene expression (Figure S4.3c), and moreover was extremely effective at assessing the outcome of unsupervised clustering, including quantifying overlap in gene expression between clusters and identifying gene pairs that were mutually exclusively expressed between clusters, even when they were highly related (Figure S4.4).

This analysis demonstrates that we could use the scL-metric to delineate, in an automated fashion, whether genes reported in the literature as associated with multiple clusters in single-cell data were truly associated with multiple cell types, without parameter tuning and using only the cell-by-gene table.

The L-metric captures the spatial distribution of gene expression despite being calculated without spatial information

We noticed that sometimes, even genes that had been annotated as belonging to the same cell type showed mutually exclusive expression (scL-metric = -1) (Figure S3.1b). We noticed that the spatial distribution of expression of the genes in these pairs was often quite distinct, and often did not match the spatial cell type distribution as well. Thus, we wondered whether an overall spatial organization, compatible with but independent of cell typing, would emerge from clustering genes based on the degree of mutual exclusivity of their expression.

For each gene, we calculated a vector of scL-metric values with all other genes, effectively cataloging the extent of mutual expression across our dataset. We clustered these vectors; a heatmap of this clustering is shown in Figure 3g. Broadly speaking, the genes associated with each cell type clustered together, and the shared category genes generally clustered closely to the NMP genes. Interestingly, although both the NMP and spinal-cord exclusive genes were each found in one large type-specific cluster (Figure 3g ii., v.), the presomitic mesoderm exclusive genes were split into two groups (Figure 3g i., iv.). There was also a subset of NMP-exclusive and spinal cord-exclusive genes that clustered together and separately from the rest of their group (Figure 3g iii.).

To determine what could be driving these differences, we plotted where the genes were expressed in space. Since the clustering was hierarchical, we first made a tree representing this hierarchical relationship, and then at each node of the tree, we made a summed density plot of all the genes contained below that node. The density plots at the top of the tree show broad distributions of expression, but as we followed the branches down closer to the leaves, we saw fascinating spatial patterns begin to emerge. We were surprised to find that the L-metric, which doesn’t inherently contain any spatial information, was able to capture a spatial tree of gene expression.

Most of the NMP-exclusive genes were found very close to the posterior and were diffusely expressed. However, in the mixed cluster, which contained both NMP-exclusive and spinal cord-exclusive genes, we saw a different pattern — here the NMP genes were more patchily expressed and extended far more anteriorly (Figure 3h i..), suggesting that genes in this group may mark cells destined to become part of the spinal cord. Indeed, one such gene was Epha5, an ephrin receptor differentially expressed in a cluster labeled as NMPs by single-cell sequencing [2] but which has separately been annotated to be involved in neural development [35] We also saw a clear division between more anteriorly and more posteriorly expressed spinal cord genes (Figure 3h ii..) Thus, clustering L-metric values allowed us to delineate differences in behavior between cell type-associated genes: one cluster of spinal cord and NMP genes appears to mark neurally fated but uncommitted cells (Figure 3h i..), while another marks committed spinal-cord fated cells (Figure 3h ii..).

We also looked more deeply into the split presomitic mesoderm exclusive genes. Not only were they expressed in different parts of the gastruloid (compare Figure 3h iii., iv.), we realized that these genes were associated with presomitic mesoderm cells in different states: Mesp2 is turned on in a specific and spatially constrained location at the onset of somitogenesis [3], a process which has only just started in this gastruloid, while Dll1 more generally marks the presomitic mesoderm. Thus, using only differences in scL-metric values, we were able to identify both general presomitic mesoderm versus presomitic mesoderm that was starting to undergo somitogenesis.

These analyses demonstrate that information contained within the hierarchical relationships between genes determined by scL-metric reveals novel information about cell states within cell types and can identify distinct spatial locations of cells in this cell state, all without explicit encoding of spatial information, but rather quantifying and clustering the degree to which genes are mutually exclusively expressed with one another.

Clustering scL-metric vectors clearly resolves cell types and reveals novel genetic interactions

Given the amount of spatial and state information that was encoded in the scL-metric heatmap for a subset of our gene panel, we expanded our analyses to all genes, hoping to discover new genetic interactions or refine existing ones. We first calculated the scL-metric value for all genes in all gastruloids, then averaged across gastruloids and clustered the resulting interaction vectors (see Methods for details). The heatmap is shown in Figure 4a (heatmap including cell cycle genes is shown in Figure S4.1a). We noted that just as when we clustered genes associated with NMPs and their direct descendants, genes associated with cell types tended to cluster together. Specifically, NMP, spinal cord, endoderm, and endothelial genes clustered very strongly together, while presomitic mesoderm genes again were split into two groups, one of which was more closely associated with genes involved in early somitogenesis.

Figure 4:

We examined the largest cluster of presomitic mesoderm genes (Figure 4a, i..). Just as before, the clustering reflected distinct spatial distributions; a highlight of one diagonal block with two subclusters is shown in Figure 4b, along with the expression patterns in a representative gastruloid. The right subcluster has most of the presomitic mesoderm genes, and is particularly enriched in the left side of this gastruloid, despite the fact that there is a cluster of presomitic mesoderm genes on both sides of the central axis (Figure S1.2a). The other subcluster has genes associated with differentiation into somites, and contains Snai1, which is a gene associated with EMT and cells in a more migratory state, consistent with how somitogenesis is thought to occur in gastruloids [2,3,24]. These genes cluster together because both sets are strongly expressed in a central region of presomitic mesoderm, one is more posteriorly biased, while the other is more anteriorly biased.

The genes associated with more anterior cell types, like somite, paraxial mesoderm, and the differentiation front (which marks the boundary between posterior and anterior), also clustered together in one large block, although the genes associated with the different types tended to be mixed (Figure 4c). When we plotted their spatial distributions, the reason for this separation became clear — one block was genes more strongly expressed in the most anterior of the gastruloid, while the other, including almost all of the differentiation front genes, was expressed more posteriorly; these genes were associated with the onset of somitogenesis. We saw similar spatial patterning in the cluster of NMP genes. One sub-cluster, which contained genes like Wnt3a, Hoxc10, and Evx1, was very specifically far-posteriorly localized, while the other, which contained genes like Fgf8, Cdx4, and some of the shared NMP+presomitic mesoderm genes, extended further into the anterior (Figure S4.1b).

We also observed interesting differences between the subclusters of endothelial genes (Figure 4d). Cldn5 had the strongest average association with other endothelial genes and strong negative associations with genes associated with other cell types, indicating that it is likely the best marker in our panel for generically marking all endothelial precursors. The cluster that contained Cldn5 showed strong intra-cluster associations (high scL-metric values), and although the other cluster showed relatively high association with Cldn5, the average association within the cluster and with the genes in the other cluster was weaker. However, this cluster contained some interesting genes: Sox17, an endoderm marker, and Tgfb1 (Tgfβ), which is associated with reduction in inflammation, pro-growth signals, and differentiation [36]. Plotting the expression of these genes in representative gastruloids confirmed their co-expression with Cldn5, a consistent marker for endothelial cells (Figure 4e).

Tgfβ appeared to be exclusively associated with Cldn5-expressing cells, but only in some gastruloids (it was extremely lowly expressed in others and not clearly associated with any specific cell type). Tgfβ loss causes embryonic lethality due to hematopoiesis and vascular development failures [37], specifically in extraembryonic tissues. Extraembryonic endoderm has been inconsistently reported in gastruloids; some studies find evidence of this tissue [2], while others do not [13]. Moreover, Tgfβ can drive the transition from endothelial precursors to hematopoietic stem cells [38], so it is particularly exciting to see it specifically expressed in endothelial tissue. While we note that the regions expressing Tgfβ have unique morphology (Figure S1.2a, xv), further work is needed to determine the role this gene plays in the development of these tissues in gastruloids. These results demonstrate that the scL-metric can be used to identify novel cell states in subsets of cell types.

To validate the clustering produced by the L-metric, we compared our results to a state-of-the-art method for identifying gene programs in an unbiased fashion from single-cell data: consensus non-negative matrix factorization (cNMF) [39]. We pooled nuclei from all individual gastruloids and ran cNMF. We found that many of the resulting clusters (Figure S4.2a) corresponded to the clusters identified when the scL-metric tree was truncated to produce exactly the same number of clusters (Figure S4.2b). The similarities were even greater when the scL-metric clusters were hand-selected based on visual inspection of the tree and density of marker genes (Figure S4.2c). From these observations, we conclude that the two methods are capable of producing similar results at a high-level, but are different in their application. Individual cells receive component scores for cNMF gene programs, yielding more per-cell information, while the tree produced by L-metric clustering reveals hierarchical information about gene programs, which quantifies their similarity in expression on a more global scale.

Thus, by quantifying the similarity between genes based on the mutual exclusivity of their expression with all other genes, we were able to extract similarity in the expression of cell type genes without explicitly specifying them. Moreover, the spatial relationship between genes, including among genes associated with the same cell type, was encoded in the tree of relationships produced by clustering. Finally, we discovered novel functions of genes in the panel through their location in the scL-metric tree: although Tgfβ was initially included in our panel to generally detect inflammatory and growth signaling, clustering by expression patterns revealed its unique association with endothelial precursors. Together, these data demonstrate the power of mutually exclusive expression relationships to uniquely quantify relationships between genes and the cells that express them.

Calculating the L-metric on smoothed gene expression spatial profiles reveals higher-order tissue organization and unique cell type arrangement in certain gastruloids

Our analysis of mutually exclusive gene expression on a per-cell basis both yielded the tools to identify highly specific marker genes and revealed gene expression states that were characteristic of subsets of our annotated cell types. However, there are examples where graded gene expression in tissue, independent of the individual cells expressing the gene, is biologically meaningful: zonation in tissues such as intestine [40], liver [41] and kidney [42], and hot and cold areas of tumors [43]. Furthermore, many spatial analysis techniques are not resolved to a single cell. Therefore, we wondered whether we could adapt the scL-metric to identify patterns of gene expression at different length scales independent of cellular assignment.

To do so, we first noted that in a spatially-resolved dataset, the assignment of spots to nuclei is essentially a non-linear partitioning of genes in space — by being assigned to a nucleus, an individual expression spot loses its specific spatial location and is binned into a region defined by the cell boundary. An example of this binning is shown on the left side of Figure 5a. We reasoned that we could instead capture relationships between the purely spatial expression of genes by mapping each transcript onto an evenly spaced grid and then smoothing the resulting 2-dimensional distribution using kernel density estimation. The bandwidth of smoothing can be set in any units, but we found the average (per gastruloid) nearest-neighbor distance to be an intuitive and useful length scale to use. Figure 5a shows an example of how choosing different distances for smoothing captures different scales of patterning—any amount of smoothing that causes gene expression in adjacent cells to be averaged would increase the L-metric by increasing the degree of co-expression, and more smoothing yields a further increase in L. We chose to use density estimations with a distance value of 2 times the nearest-neighbor average distance, as this spacing seemed likely to capture meso-scale organizational patterns in gastroids while preserving finer structures of spatial variance (Figure 5a).

Figure 5:

The scL-metric is calculated by comparing gene expression on a cell-by-cell basis. To adapt this measurement to spatial data, we calculated the L-metric using the density estimation at each spatial bin across the entire gastruloid (i.e., each bin is treated the way a cell is treated when the scL-metric is calculated). We denote this value as the “spatial L-metric”.

We calculated the spatial L-metric for all the consistently-detected genes in our panel across all gastruloids (202/220 genes, see Methods). When we averaged across gastruloids and clustered the values, the resulting heatmap (Figure S5.1a) had many properties similar to the scL-metric heatmap shown in Figure 5a: For the most part, cell types clustered together; however, we noted a few unique patterns. The order of clustering seemed to more closely reflect overall position in the gastruloid: the more posterior cell types (NMP, presomitic mesoderm, spinal cord) were found close together, then endoderm, differentiation front, and a mixed cluster of paraxial mesoderm, somite, and endothelial genes. This ordering suggests that when the spatially smoothed relationships (spatial L-metrics) of all genes in the panel are clustered, the resulting relationships reflect the overall spatial organization of the underlying cell types within the larger tissue structure. Individual cell types appeared to cluster more strongly when using scL-metric (compare Figure 4a with Figure S5.1a), but the order and mixing of cell type genes with the spatial L-metric reflects tissue-level organization—NMP and presomitic mesoderm genes were close together and mixed, reflecting the mixing we see of the cell types associated with these genes.

The other difference of note was that, unlike in the scL-metric heatmap, the endothelial genes were much less strongly clustered, likely reflecting the diversity of locations they were found in within individual gastruloids (Figure S1.2). We looked at their patterns in the heatmaps for individual gastruloids and noticed that there was strong co-clustering of endothelial and endoderm genes in some of the gastruloids. An example is shown in Figure 5b, which is the clustered heatmap for an individual gastruloid (cell types and morphology shown in Figure 5c). A strong block of interaction on the diagonal (“endoderm/endothelium cluster”, outlined) indicates a strong association of these genes on a length scale beyond individual cells. When we looked at the cell type arrangement of this gastruloid, there were two clear populations of endothelial cells: those that were located in the anterior and those that were more posteriorly located and strongly intermixed with the endoderm tissue. Thus, the unique clustering behavior of genes in this gastruloid is reflective of a unique, tissue-level pattern that is only found in some of our samples. When we calculate the scL-metric heatmap for the same gastruloid, although endothelial and endoderm genes cluster close together, they form two distinct diagonal clusters, compared to the single, high-L block found in the spatial-L heatmap (compare the outlined clusters in Figure 5b and Figure S5.2a).

These comparisons suggested that the spatial L-metric could capture patterns on a tunable length scale. When we directly compared the intra-cell type L-metric values calculated using both methods (Figure S5.1b), we saw that while the values for some cell types, like NMP, were quite similar, others, such as spinal cord, showed a dramatic shift towards +1 when calculated using a smoothing of 2 nearest neighbors. This shift is consistent with our observation of consistent expression of NMP genes and a more “salt-and-pepper” expression of spinal cord genes, and suggests that by calculating the spatial L-metric over many length scales, one could learn the effective length scale of interaction for specific gene pairs.

Endothelial precursors show unique organization and distinct, spatially-dependent cell states

Throughout our analyses of both spatial organization and gene expression, endothelial cells emerged as displaying particularly heterogeneous patterning and gene expression that depended on their physical location and neighboring cells, which varied between gastruloids. The relative rarity of this type has precluded in-depth analysis of their location and function in other studies [2,13,14,44], and despite strong interest in modeling hematopoiesis and the development of vasculature [45], it is still not clear to what degree gastruloids produced by standard protocols create and organize endothelial tissue. Therefore, we more specifically examined the endothelial cells identified in our gastruloids, profiling their spatial location, organization, and gene expression.

We first used density-based spatial scanning and clustering [46] to, in an unbiased fashion, detect clusters of like cells in individual gastruloids that by eye appeared to contain small clusters of endothelial precursors in the anterior. These clusters are highlighted in Figure 6a. In some cases, we observed that the cells in the cluster formed a ring-like structure (Figure 6a highlight). We applied the same spatial clustering method to all other cell types and found that, on average, the clusters of endothelial cells were the most circular (Figure 6b, Figure S6.1a). This clustering behavior is consistent with the peri-gastrulation behavior of endothelial precursors in real embryos, which have specific migratory activity [47] and form aggregated clusters that go on to form endothelial and hematopoietic cells [48]. We detected endothelial clusters in all of the gastruloids (Figure S6.1c), and found that they formed more clusters than all other cell types except paraxial mesoderm (Figure S6.1b). Importantly, the relative location of these clusters along the AP axis varied substantially between gastruloids (Figure 6c, Figure S6.2a), highlighting the necessity of studying individual gastruloids to observe the unique behavior of this rare cell type.

Figure 6:

We also observed that in 5 out of our 26 samples, there was a large central patch of endoderm cells intermixed with endothelial precursors; these samples also had unique spatial L-metric clustering of endothelial and endoderm genes (Figure 5b). An example of one such gastruloid is shown in Figure 6d. Migration to and association with the endoderm is also a hallmark of endothelial development [47,48], and we were curious whether there were differences between these cells and the cells we observed forming anterior, somite-associated clusters. The endoderm-associated endothelial cells were generally in larger clusters, and when we computed the cell type exposure index for just this gastruloid, we found that, consistent with our visual observations, in this particular sample, endothelial and endoderm cells were much more frequently found next to one another than on average (Figure S6.1d). To determine whether these spatial and organizational differences reflected gene expression differences, we divided the gastruloid normal to the anterior-posterior axis to separate the endothelial cells into endoderm-associated and somite-associated and looked for differentially expressed genes between the two groups. The significantly differentially expressed genes we found are shown in Figure 6e.

We found that some endoderm genes were more abundant in endoderm-associated endothelial cells and that some somite genes were more abundant in somite-associated endothelial cells, likely reflecting the mis-association of genes from neighboring cells, which is common in nuclear-based spot assignment in highly packed tissues like gastruloids. However, we also found two endothelial genes that were differentially expressed in endothelial cells from the two different regions: Pecam1 and Cdh5, both of which are associated with angiogenesis, were more enriched in endoderm-associated endothelial cells, and Cxcr4, which is thought to aid haematopoiesis, angiogenesis, and cardiac cell organization by responding to homing signals. Pecam1 also clustered uniquely in our L-metric analysis (Figure 4c), suggesting this differential expression is conserved across gastruloids. Spatial expression of these genes is shown in the bottom row of Figure 6f. Although most endothelial cells are thought to be of mesodermal origin, some evidence suggests that, in the organogenesis of specific tissues like the liver, the endoderm can give rise to endothelial cells [49]. Our data demonstrate a molecularly driven organization distinct from the clustering we observed in the anterior and suggest that multiple mechanisms of endothelial specification could be modeled in gastruloids, even simultaneously within the same structure, although further characterization is needed to determine exactly what processes these unique endodermal/endothelial structures model.

The genes enriched in the anterior have less obvious functional distinctions; for example, it isn’t immediately clear how a gene like Apoe, which has been annotated to be associated with primordial germ cells, specifically functions in anterior endothelial precursors (Figure 6f, top row). Gadd45g, a stress-associated protein, is also differentially expressed, possibly reflecting the fact that cells in the anterior are exposed to more rapid growth, compression, and cell death.

Although this cell type has been consistently observed in single-cell measurements of gastruloids, its relative rarity has precluded in-depth analysis of subtypes or inference of spatial location. Our results strongly suggest that endothelial precursor formation, migration, and organization may all be modeled in 3D gastruloids; recent advances in 2D gastruloids have allowed modeling of cardiac and hepatic vascularization [45], and our data suggest that 3D gastruloids may similarly be adapted to model more specific aspects of hematopoiesis and vascularization. Early specification from a pool of mesodermal precursors is a hallmark of the endothelial lineage [47]; given the consistency with which we observe endothelial precursors, we speculate that this behavior is recapitulated in gastruloids, but further epigenetic measurements are required to validate this hypothesis.

Discussion

One of the central goals of development is to map how the identity and arrangement of cells are determined and coordinated. Although precise scaling of gene expression in gastruloids had been previously reported [10], the relationship between this scaling and the proportions of specific cell types or more complex tissue organization was not known. Towards that end, we have performed detailed measurements and analyses of the expression of hundreds of genes in tens of individual gastruloids, revealing how cell types are organized within and across gastruloids. Our results highlight that many aspects of embryonic development, such as the arrangement and differentiation trajectories of posterior cell types, are faithfully recreated on a per-gastruloid basis. Others, however, can vary, and conditions may need to be tuned to increase reproducibility if a specific aspect of their development needs to be modeled. For example, we reproduced observed variation in endodermal structures [12] and reported both the distinct organization of endothelial cells and variation in that organization between gastruloids. Although we do not fully reproduce the level of per-gastruloid cell type proportion variation reported in a previous single-cell sequencing dataset [13], we do find evidence of a tradeoff between neural and mesoderm tissues; differences in the number of cells analyzed per gastruloid could potentially explain this discrepancy. Moreover, we observe many of the same cell types as in this study, and show that their organization is quite consistent.

We introduce the L-metric, which quantifies mutually exclusive gene expression, and highlight two ways of using it that can be employed on different types of data. The cell-based L-metric, scL-metric, can be applied to any technique that generates a cell by gene table. Hierarchical clustering analysis of the resulting vectors can reveal how similar each gene is to others in terms of how it is co-expressed with all genes, and may reveal spatial information if the sample was collected from a spatially-organized tissue. We were able to verify the spatial information because our data were spatial, but the scL-metric is calculated only using the cell-by-gene table, a standard output for any single-cell analysis technique.

The spatial L-metric, which operates spatial bin-wise rather than cell-wise, can be used on spatial data even when cell assignment or nuclear segmentation isn’t possible — in principle, it can reveal which genes are patterned similarly at many different length scales, with the fundamental resolution limit set by the measurement technique. It was able to identify co-expression, even amongst lowly-expressed genes, such as the spinal cord markers, that had a more ‘salt-and-pepper’ expression pattern and were only strong markers when considered in aggregate. It is possible that some genes with bursty gene expression may appear not to express together, but spatial smoothing may allow for the identification of co-expression of these genes, which may be useful for cell type identification.

Although the identification of cell types has utility for high-level analyses of organization, our in-depth dissection using the L-metric revealed that, particularly for cells in the process of differentiating, a gene-centric view may yield a more holistic picture of tissue organization. It is particularly exciting to imagine these analyses performed in a time-resolved fashion, as dynamic changes in gene expression across a tissue in response to physical rearrangement and broad signaling patterns could reveal mechanisms or spatial aspects of differentiation not visible when cells are merely binned into predetermined categories.

Methods

Cell culture

All gastruloids were generated from E14TG2a cells (a generous gift from the Toettcher lab). Cells were cultured in 2i+LIF media at 37 °C, 5% CO₂ and grown on plates coated in 0.1% gelatin (Millipore/Sigma ES006B). 2i+LIF media consisted of GMEM (Gibco 11710-035), 10% stem-cell qualified FBS (R&D Systems S10250), 1x GlutaMAX (Gibco 35050–061), 1x MEM non-essential amino acids (Gibco 11140050), 1 mM sodium pyruvate (Gibco 11360–070), 100 μM 2-mercaptoethanol (Gibco 471 21985–023), and 100 units/mL penicillin/streptomycin (Gibco 15140–122), 1000 units/mL LIF (Millipore Sigma ESG1107), 2 μM PD0325901 (Tocris 4192), and 3 μM CHIR99021 (Tocris 4423). Cells were passaged every 2 days: cells were trypsinized (Life Technologies 12563011), counted using an automated cell counter (Invitrogen Countess II), and 5x10^-5 — 8x10^-5 cells were transferred to a 25cm² dish; cells were passaged at least once after thawing before gastruloid aggregation.

Gastruloid aggregation

The aggregation protocol was as described in [2] (protocol here). Briefly, cells were trypsinized, rinsed twice with 1X PBS to remove media, and then resuspended in N2B27 medium. N2B27 consisted of equal parts neurobasal medium (Gibco 21103049) and DMEM/F12 (Gibco 11320033) supplemented with 1X N2 (Gibco 17502048) and B27 (Gibco 17504-044), 100 units/mL penicillin/streptomycin (Gibco 15140–122), and 100 μM 2-mercaptoethanol. Cells were resuspended at a concentration of 2500-7500 cells/mL; 40 μL of this mixture were placed in each well of round-bottom ultra non-adherent U-bottom 96-well plates (Costar 7007), yielding 100-300 cells/well. The first two seqFISH runs used 300-cell gastruloids (n=8 total) and the last used 100-cell gastruloids (n=18). Cells were allowed to grow at 37 °C, 5% CO₂ for 48 hours and then supplemented with 150 μL fresh N2B27+3 μM CHIR99021 (Chiron, Wnt agonist). After 24 hours, 150 μL of media was removed and replaced with N2B27 without Chiron; this media change was repeated daily until the gastruloids reached maturity at 120 hours. All wells were scored for elongation at 120 hours just before fixation. The gastruloids selected for fixation were elongated and uniaxial; the experiments that generated them had an 80-94% success rate (see Figure S1.2d for full distribution).

Sample fixation

After gastruloids had grown for 120 hours, they were harvested by pooling together all successful (i.e., elongated and uniaxial) wells in a plate, washed with 1X PBS, and then fixed for 10 minutes in 4% PFA in PBS. After fixation, gastruloids were washed with 1X PBS and stored in methanol for up to a month before preparation for seqFISH. 3 separate seqFISH runs were performed on 6 different aggregations (one experiment pooled 4 plates of gastruloids), producing a total of 26 individual samples.

Embedding and sample preparation for seqFISH

∼30 gastruloids were placed in a Lo-Bind tube (Eppendorf 0030108523). After removing as much methanol as possible, the gastruloids were washed with 2XSSC, then incubated in 2XSSC + 1 μg/mL DAPI overnight at 4°C. Stained gastruloids were permeabilized with 1X PBS + 0.2% Tween20, arrayed on a gelatin-coated slide (Spatial Genomics 00200001) with an eyelash tool (Electron Microscopy Sciences 11852), covered in liquid acrylamide gel, and allowed to polymerize for 2 hours. Polymerized samples were washed with 8% SDS, then treated with 8x10^-3 units/mL Proteinase K (NEB P8107S). After this wash, the primary probes were hybridized for 16-40 hours, after which the samples were washed and loaded onto a GenePS machine (Spatial Genomics).

Detailed protocol available here.

seqFISH and raw data handling

Sequential rounds of secondary probe addition and imaging were performed in an automated fashion by the instrument [16]. The GenePS software that controlled the instrument was v1.6.0. The instrument was run with a 1 or 2 μm z-offset on 5-18 individual regions of interest, highlighting each gastruloid. The genes chosen for the seqFISH panel were based on differential expression in single-cell RNA sequencing experiments in references [2,13] or by consensus in the literature. 200 of the genes were detected via a sequential barcoding scheme (see [16] for details), and 20 were serially detected as in traditional smRNA-FISH [50]. The serial genes tended to vary in quality between gastruloids; genes with high background or in which expression was too high for the spot detection software to confidently assign spots were excluded. The list of excluded genes for each gastruloid is included with the raw data.

Gastruloid annotation and spot assignment

Detected spots were decoded and exported without nuclear assignment by SG Navigator software version 0.9.0. At this stage, gastruloids without even and consistent expression of the housekeeping gene Eef2, which indicated focus or technical issues or sample variability, were discarded. 5 out of 6 gastruloids were discarded from the first run, and 3 out of 11 in the second run. None of the samples from the third run were discarded. Raw TIFF files of the DAPI images were annotated using the NimbusImage platform, which is available here: https://github.com/arjunrajlaboratory/NimbusImage/ and hosted here: https://www.nimbusimage.com/. Expression of the gene T was used to locate the posterior end and corroborated by consistent differences in nuclear morphology between the posterior and anterior. An AP axis was inferred by drawing a line between that pole and the extreme opposite end of the gastruloid along the center of the gastruloid parallel to the edges and connecting with the furthest extended point to define the anterior pole. Segmentation was performed using cellposeSAM [51] on the high-resolution .tiff files of the DAPI-stained nuclei from the 3rd round of hybridization. Once the nuclear masks had been created, spots were assigned to nuclei using a custom software package called SGAnalysis (available here). Nuclei were dilated by 8 pixels (0.86 μm), and all spots that fell within this dilated nuclear area were assigned to that nucleus. In cases where the dilated regions overlapped, spots were assigned to the nucleus for which the edge of the dilated area was furthest away. On average, we were able to resolve 131,112 nuclei for each gastruloid, and of those nuclei, 79,607 had enough mRNA counts to be typed. Nuclei without enough mRNA were generally small, uniformly bright (in contrast to the usual mottled appearance of mouse nuclei), and lacked expression of the highly-expressed housekeeping gene Eef2, indicating the corresponding cell was likely dead or fragmented. We primarily found these nuclei in the anterior tissues and classified them as “none”-type, excluding them from our analyses unless noted. The average reads/nucleus per gastruloid ranged from 47-80; for typed nuclei the value was 105-265.

Cell type scoring

To assign cell type, we used an panel of marker genes for 9 consensus cell types thought to exist in gastruloids at this stage of development: Neuromesodermal precursor (NMP), presomitic mesoderm (PSM), spinal cord precursor (spinal cord), endoderm, differentiation front, cardiac mesoderm, paraxial mesoderm, somitic, and endothelial precursor (endothelium). Genes were chosen based on differential expression in references [2,13], or consensus in the literature. Each nucleus received a score for each cell type using the score_genes function in scanpy [52], and the highest-scoring cell type category was assigned to that cell. If none of the scores met a universally applied cutoff (0.1), then the cell was assigned ‘none’.

Cell type entropy score

The degree of uncertainty in cell typing was quantified by calculating the entropy of the cell type score distribution. The scores can be negative or zero, so they were first converted to probabilities using a softmax function with temperature scaling:

A temperature of 0.1 was used. Once the scores had been converted, the Shannon entropy of the distribution was calculated:

Mixing score calculation

The mixing score for a gastruloid was calculated by finding the nearest neighbor of each cell and determining how many cells had a nearest neighbor of the same type. This value was normalized by dividing by the expected number of same-type nearest neighbors based only on the frequency of each type in that gastruloid. Values range from 0 (perfectly mixed) to 1 (perfectly segregated).

Exposure index

The exposure of each cell type was calculated by finding the k nearest neighbors (“target”) of every cell (“source”) in every gastruloid (k=5 was used), counting the frequency of each cell type among these neighbors, and then normalizing based on the overall frequency of occurrence of the neighbor (“target”) cell type. This generated a pair-wise exposure score of each cell type to self (symmetric) and all other cell types (asymmetric) that ranged from 0 (never found next to one another) to 1 (always found next to one another).

Triplet identification and enrichment

Spatial enrichment of cell type triplets was assessed using a computational framework based on Delaunay triangulation. Each gastruloid was processed independently to account for potential inter-sample variability in cell density and spatial organization. Spatial neighbors were identified using Delaunay triangulation, which partitions the coordinate space into triangular tessellations connecting nearest neighbors. A distance cutoff of a maximum of 200 pixels was used; neighbors with greater than this distance were not considered. All groups of 3 nearest neighbors were used to evaluate triplet frequency. Statistical significance was evaluated through a permutation testing framework. For each sample, the observed triplet frequencies were compared against null distributions generated through random reassignment of cell type labels while preserving spatial coordinates. The number of permutations was set to 1,000 per sample. The data for all gastruloids were compiled together such that the significance of each triplet in each sample, and how many samples each triplet was found in, could be assessed. Only significantly enriched or depleted triplets were analyzed; significant triplets found in at least half of the gastruloids are shown in Figure 2b, and all are shown in Figure S2b.

Kernel density estimation for generating smoothed spatial profiles of gene expression

We performed kernel density estimation (KDE) to generate smoothed spatial profiles of gene expression patterns. Using the nonparametric.KDEmultivariate method from the statsmodels package in Python, a given gene’s KDE of gene expression in a given gastruloid was calculated using a two-dimensional Gaussian kernel for each RNA spot fitted to a bandwidth equal to the average nearest neighbor distance between cell centroids in that gastruloid multiplied by 2 (Table S1) in both dimensions. To produce an expression profile for a gene from its KDE on a given gastruloid, we first partitioned the gastruloid’s annotated hull to produce a grid containing square elements, each with side length equal to half of the average cell diameter, and then evaluated the KDE on the center of each grid element. Visualizations of multiple genes’ combined expression in a given gastruloid were obtained by summing individual genes’ normalized KDE values for each grid element.

Calculating the single-cell L-metric

The scL-metric for a given gene pair was calculated per-gastruloid using a cell-by-gene matrix filtered to exclude cells with less than 2 transcripts for both genes and sorted to place cells in descending order according to the first gene’s expression values. The scL-metric is asymmetric with respect to pair order in that the first gene provides the reference expression distribution to which observed and simulated expression distributions of the second gene are compared. On the ordered and filtered cell-by-gene matrix, the following three simulated expression distributions were created for the second gene: perfect coexpression, in which the second gene’s expression values were sorted in descending order, like the first gene; perfect anti-coexpression, in which the second gene’s expression values were sorted in ascending order; and, uniform coexpression, in which all cells were assigned the mean expression value of the second gene. Cumulative sums were calculated for the observed expression distributions of the first gene and second gene and the simulated expression distributions of the second gene. The cumulative sums of each expression distribution associated with the second gene were plotted against the cumulative sums of the first gene, and each of the resulting curves’ areas were obtained using trapezoidal integration and normalized by the product of the first gene’s total expression and the second gene’s total expression, producing A_observed, A_positive, A_negative, and A_uniform (Figure S3.2a,b). The scL-metric value for this gene pair, specific to its ordering, was then determined by comparing A_observed to each of the areas under the simulated distributions’ curves:

We computed the scL-metric across all gene pairs for each gastruloid in our dataset from 05/07/2025 (n=18 gastruloids). We did not integrate this dataset with our other two datasets (n=8 gastruloids), because relatively poorer detection quality in those datasets limited our ability to estimate the distributions necessary to calculate the scL-metric. Because some of the serial genes in the seqFISH panel were inconsistently detected between samples (see seqFISH section of Methods), the number of gene pairs varied slightly between samples. When an scL-metric matrix for an individual gastruloid is shown, the number of genes for that gastruloid is indicated. We averaged together the scL-metric matrices for the 18 individual gastruloids. For this analysis, a common set of 202 genes with good detection in all gastruloids was used.

All scL-metric analyses, with the exception of hierarchical clustering analysis, use a symmetrized form of the scL-metric obtained by averaging the two directional scL-metric values for a given gene pair. This symmetrization was motivated by the aim to summarize the coexpression relationship between two genes with a single value and the lack of significant differences between directional L-metric values.

Cell type analysis via unsupervised clustering

We used scanpy [52] to process the dataset from 05/07/2025. All nuclei were pooled together, nuclei with less than 40 total gene counts were filtered out, and then PCA (50 PCs) and UMAP embedding were performed using the standard pipeline. We performed Leiden clustering using scanpy.tl.leiden, and used visual inspection to choose a resolution at which to cluster. A resolution of 0.8 produced fairly well-separated clusters and roughly corresponded to the expected number of cell types. Clustering results and spatial projection are shown in Figure S4.4a.

Filtering marker genes

We re-scored nuclei using a refined marker gene panel. We retained only genes that met the following criteria: more than 10% of all scL-metric values < -0.8 (73% of all genes) and more than 12% of all scL-metric values > -0.3 (90% of all genes). 70% of the genes in our panel met both. Scores with this filtered gene panel are shown in Figure S4.3b.

Assessing marker gene panels and cluster-associated genes

To visually represent the relatedness of the genes either differentially expressed in clusters identified with unsupervised clustering, or the groups of genes associated with cell types in a marker gene panel, we calculated the average L-metric score of all possible combinations where the first gene in the L-metric score belonged to the first type or cluster, and the second gene belonged to the second group or cluster. In the case of the marker gene panel, we first filtered genes on expression level, considering only genes that had > 2.3 counts / nucleus in nuclei where they were expressed; this filtering retained 84% of the genes in the panel.

Using cNMF to find gene clusters

All nuclei were pooled together, and cNMF [39] was run on all nuclei and all non cycle-genes. The default parameters were used and 5-10 components were calculated. 7 components produced the highest stability, and the top 24 most positively differential (by ordered z-score) genes were used for visualization.

Calculating the spatial L-metric

The spatial L-metric was computed analogously to the scL-metric except that rather than rank ordering the expression of genes in individual cells, we instead fitted a KDE over the detected spots for a given gene and evaluated it on a square grid of ’spatial bins’ that spanned the entire gastruloid (as described above), then rank-ordered these spatial bins by the resulting probability densities of expression. We treated these bins the same way we had treated cells in the previous calculation, so that the resulting L value reflected mutually exclusive expression in space rather than per-nucleus. The length scale of the interactions considered was determined by the chosen size of the spatial bin. Throughout, we report values for bins at 2 times the nearest-neighbor distance (which was calculated on a per-gastruloid basis). We did not filter low-expressing bins the way we had for cells for 2 reasons: 1) we calculated the KDE on all spots detected, not just those that fell within a nuclear area, so there were fewer areas of low expression, and 2) smoothing resulted in a non-zero density throughout the imaging area. Just as for the sc-L-metric, we symmetrized the spatial L-metric by averaging the 2 values calculated for each gene pair, except when displaying the values as a heatmap.

Hierarchical clustering of L-metric matrices

We performed hierarchical clustering on L-metric matrices to group genes based on their pairwise coexpression relationships. Hierarchical clustering of a given L-metric matrix was done by applying the cluster.hierarchy.linkage method from the SciPy package in Python with parameters method and distance set to “ward” and “euclidean,” respectively, on the rows of the matrix, which represented interactions between a given gene (as the first gene in the pair) and all other genes. We performed this hierarchical clustering analysis on the following gene sets: 1) A subset of only NMP, presomitic mesoderm, and spinal cord genes, evaluated for 1 gastruloid (n=36 genes, Figure 3g); 2) a subset of all genes well-detected in all gastruloids, but not including cell cycle genes (n=166 genes, evaluated on the average value across all gastruloids in Figure 4); 3) a common set of 202 genes well-detected in all gastruloids (evaluated on the average L-metric values across all gastruloids in Figure S4.1a and Figure S5.1); and, 4) a subset of genes in one an example gastruloid excluding poorly-detected genes and cell cycle genes (n=171 genes, Figure 5b and Figure S5.2a).

Cell cluster identification

Cell clusters were identified using DBSCAN [46] implemented in sklearn.cluster. A maximum distance of 300 pixels was used as a cutoff for neighbor identification, and a minimum of 8 cells was required to seed a cluster.

Cluster circularity

Cluster circularity was calculated by taking the convex hull of the identified cluster, then comparing the shape of that convex hull to a circle. The area (A) and perimeter (P) of the convex hull were measured, and then the following formula was used to calculate circularity:

If C=1, the shape was a perfect circle. Values less than 1 indicated deviations from circularity.

Supplementary figures

Figure S1.1:

Figure S1.2:

Figure S1.3:

Figure S2:

Figure S3.1:

Figure S3.2:

Figure S4.1:

Figure S4.2:

Figure S4.3:

Figure S4.4:

Figure S5.1:

Figure S5.2:

Figure S6.1:

Figure S6.2:

Table S1.

Table S2.

Data availability

All data and code used to process the raw data and generate the figures can be found at the following link: https://www.dropbox.com/scl/fo/bchkqlbcjb8ub9m606did/AIudcWZaC566toXzb2L-jXc?rlkey=u0wgtkq8jerxoqb5ump599oip&dl=0. Additional custom scripts used to process the raw seqFISH data can be found on github; see Methods for specific usage and links to the relevant packages.

Acknowledgements

The authors acknowledge Evan Underhill, Harry McNamara, and Jared Toettcher for the generous gift of cell lines, protocols, and training, Sedona Murphy for sample prep and panel design help, Sedona Murphy, Harry McNamara, Grant Kinsler, and Vinay Ayyappan for helpful discussions, and all members of the Raj lab for emotional support and feedback on the manuscript.

C.T. is funded by a Damon Runyon Postdoctoral Fellowship (DRG-2465-22). P.S. acknowledges support from the Roy and Diana Vagelos Scholars Program in the Molecular Life Sciences and the Roy and Diana Vagelos Science Challenge Award. A.R. acknowledges support from a center grant from the Mark Foundation for Cancer Research, NIH Director’s Transformative Research Award R01 GM137425, NIH R01 CA238237, NIH R01 CA232256, and NIH 4DN U01 DK127405.

Additional information

Author contributions

C.T. performed experiments, C.T. and P.S. performed analysis and prepared figures, Y.H., P.S., and A.R. developed the L-metric and wrote associated code, C.T. and A.R. wrote the manuscript, and all authors read and approved the manuscript.

Declaration of generative AI and AI-assisted technologies

During the preparation of this work, the authors used Claude to generate and improve code for analyses and to suggest improvements to the clarity and flow of the text. The authors reviewed and edited the content and take full responsibility for the content of the manuscript.

Funding

Damon Runyon Cancer Research Foundation (DRG-2465-22)

The Mark Foundation for Cancer Research

National Institutes of Health (R01 GM137425)

National Institutes of Health (R01 CA238237)

National Institutes of Health (R01 CA232256)

National Institutes of Health (4DN U01 DK127405)