A single-cell resolution 3D molecular atlas

Single-cell RNA-sequencing (scRNA-seq) approaches have been used to profile the molecular features of individual cells during early development^12,13. However, most single-cell studies have profiled dissociated populations of cells, where spatial information of the single cells in the embryo is lost. Previous works of mapping the location of cells in biological structures on the basis of the concordance of the gene expression profile^15,19 have been confounded by mathematical uncertainties and false positives^20,21. A high-value attribute of the spatio-temporal transcriptome is the amenability of mining the dataset to identify population-specific signature transcripts as zipcodes (Fig. 1a) that could be applied for imputing the position of single cells in the germ layers of E6.5–E7.5 embryos^15,22. Leveraging our spatio-temporal molecular atlas and single-cell RNA-sequencing (scRNA-seq) dataset, we developed an analytics methodology that comprises tiered algorithms to infer the spatial distribution of single cells and to reconstruct a single-cell resolution 3D molecular atlas (Fig. 1e).

In our mapping study, the 3D corn plot template was used for mapping the spatial coordinates of single cells. By employing a multi-dimension single-cell mapping (MDSC Mapping) algorithm (Fig. 2a and “Methods”), single cells were mapped by the zipcodes embedded in their transcriptome to the inferred position in the model. To evaluate the level of precision of positional mapping, we tested the mapping of single cells isolated from known positions at five developmental time points between E6.5 and E7.5. The results showed that the single cells could be mapped to their best-fit site of origin at a significant confidence level (PCC values at 0.74–0.97, Fig. 2b, Supplementary Fig. ^4b and “Methods”). Other technologies, such as MERFISH and MERSCOPE, are amenable for capturing spatial transcriptome at single-cell resolution but have limited gene detection rate²⁰ compared with Geo-seq and therefore are less suitable for the spatial mapping of single cells.

Open in a separate window

Fig. 2

A single-cell resolution 4D molecular atlas of mouse gastrulation.

a The pipeline of multi-dimension single-cell (MDSC) mapping. The SRCCs of the expression values of the zipcodes of each single cell against all reference samples of the reference embryo were calculated, followed by the application of a spatial smoothing algorithm to impute the high-confidence (closest) location (see “Methods”). The mapping of cells to position 8P was shown as an example. b Verification of the results of MDSC Mapping of single cells isolated from a known position in an E7.0 embryo. The number on each corn indicates the number of cells mapped to the Geo-seq position in the germ layers. PCC values and confidence intervals are shown for the simulation. c Uniform manifold approximation and projection (UMAP) plot showing the data structure of the “Gastrulation Atlas” comprising 32,940 cells from E6.5 to E7.5 embryos, with the exclusion of extraembryonic cells. Twenty-five cell types are annotated (see color legend). Def. endoderm definitive endoderm, PGC primordial germ cells. d Fraction of cell type per time point, displaying a progressive increase in cell-type complexity during gastrulation. e MDSC Mapping results of exemplar cell types for E6.5, E6.75, E7.0, E7.25, and E7.5 embryos. The number in each corn indicates the number of cells mapped to the specific Geo-seq position. f The spatiotemporal distribution of all the single cells identified in the Gastrulation Atlas of E6.5–E7.5 mouse embryos. Cell types are annotated as in (c). g The spatiotemporal distribution of single cells annotated as “nascent mesoderm” in the epiblast and mesoderm of E6.5–E7.5 embryos. h The spatio-temporal distribution of Pou3f1-expressing cells in E6.5–E7.5 embryos. The color legend indicates the level of expression determined by the transcript counts.

Drawing from the single-cell transcriptome data of mouse embryos at gastrulation to early organogenesis (“Gastrulation atlas,” Fig. 2c, d)¹², we first applied the MDSC Mapping algorithm and mapped the 25 embryonic cell types to the germ layers of the E6.5–E7.5 embryos (Fig. 2e and Supplementary Fig. ^4c–l). To visualize the spatial distribution of the cells in the germ layers, an Annulus Model was applied to display the “Geo-seq position” of single cells (Fig. 1e, Supplementary Fig. ^5a–e and Methods) in a series of 3D-rendered spatiotemporal maps. After mapping single cells to the best-fit position in the Annulus Model (Supplementary Fig. ^{5d, e}), we applied mathematical modeling based on Bubble Sort Algorithm to re-position single cells within each Geo-seq position by incorporating information of vectorially graded molecular activity (e.g., gradient of signaling activity and gene expression levels) (Fig. 1e, Supplementary Fig. ^5f and “Methods”) that is known to associate with the regionalization of cell fate. This mathematical model enables an effective imputation of the coordinates of single cells in the Geo-seq zipcode position (Fig. 2f–h). By taking into consideration the pattern of distribution across the developmental stage, the data could be rendered digitally into the spatio-temporal (4D) Atlas at single-cell resolution (Fig. 2f and Supplementary Fig. ^6a). This 4D Atlas provides unique insights into the spatial distribution of, for example, specific group of single cells in the germ layers of E6.5–E7.5 embryos (Fig. 2g and Supplementary Fig. ^6a) and, in specific cases, the Pou3f1-expressing cells in the epiblast/ectoderm (Fig. 2h), and T-expressing cells in the three germ layers (Supplementary Fig. ^6b) during gastrulation.

The single-cell mapping results were generally consistent with prior knowledge of the identity of cell types annotated based on expression markers and lineage inference^11,12. For example, the anterior PS cells were mapped to the distal posterior epiblast (positions 1P–4P) at E7.0 (Supplementary Fig. ^4h), and rostral neuroectoderm cells were mapped to anterior regions of the ectoderm of E7.5 embryo (Supplementary Fig. ^4l). However, some cell types were mapped to areas of the germ layer that are apparently inconsistent with their cell identity. For example, cells of the PS were mapped widely in the epiblast outside the PS (Supplementary Fig. ^{4d, f, h, j, l}), and cells that are annotated as nascent mesoderm were mapped predominantly to the posterior epiblast (6P–11P), besides in the mesoderm layer at E7.0 (Supplementary Fig. ^4h). It might be that these cells in the epiblast represent a transitional cell state and not the fully specified cell type^23,24. This observation raised the possibility that the PS cells and nascent mesoderm mapped to the epiblast may represent a population of progenitor cells at the transition of their allocation from the epiblast to the emerging germ layers. Drawing on the knowledge of the prospective fate of cells in the epiblast of gastrulating embryos^1,11, we re-annotated the single cells mapped to the epiblast as progenitor/precursor cells (the intermediate cell types) of the germ layer derivatives based on the findings of heterogeneity of cell populations at a coordinate position in the germ layers (see next section).

Heterogeneity of space-registered cell population

To visualize the composition of the single-cell population at each Geo-seq position, an optimization algorithm based on Euclidean distance was applied (Fig. 3a and “Methods”). Using the gene-expression matrix of single cells mapped to a specific position as an input, computing a normalized Euclidean distance matrix followed by applying the least square method, an optimized coordinate was assigned to each single cell. Applying this optimization algorithm, we simulated the spatial pattern of cell types within each Geo-seq position. For example, for position-6P at E6.75, this refined mapping revealed the presence of a heterogeneous cell population of the epiblast, PS, and nascent mesoderm cells (Fig. 3b, c). This finding prompted us to examine the heterogeneity of cell types in all the Geo-seq positions in the germ layers of the embryo at the five developmental time points.

Open in a separate window

Fig. 3

Mathematical modeling for single-cell spatial distribution and the collation of spatio-temporal heterogeneity map.

a The mathematical model for re-ordering the spatial distribution of the single cells within a Geo-seq position by Euclidean-distance derived optimal coordinates (see Methods). b t-distributed stochastic neighbor embedding (t-SNE) plot showing the single cells mapped to position-6P at E6.75. Cell types are annotated (see legend). c The imputed spatial distribution of single cell types within position-6P at E6.75. Cell types are annotated as in (b). d Heterogeneity Map. Corn plots (top row) show the composition of cell types (shown as a pie chart in each corn) at different Geo-seq positions in the germ layers across the five-time points of gastrulation. Corn plots (bottom row) show the heterogeneity of cell populations descending from the primitive steak-like cells in different Geo-seq positions in the germ layers of the gastrulation stage embryo. e The molecular trajectory of the descendants of primitive streak-like cells of the E6.5 posterior epiblast in the germ layers of E6.75, E7.0, E7.25, and E7.5 embryo, imputed using the Population Tracing algorithm. Cell types are listed in Supplementary Fig. ⁸. The rectangle represents the spatialized cell type (color indicating the cell types), with the size indicating the propensity of branching trajectory, and the width of the edge indicates the strength of correlation between the connected cell types.

To uncover the heterogeneity of single-cell population in the Geo-seq position, we applied t-SNE to re-cluster cells at each position in the germ layers of embryos at the five developmental stages (Supplementary Fig. ⁷). This analysis revealed that some cell types that were previously annotated differently were transcriptionally similar, illustrating the challenges of grouping of cells into discrete entities by primarily transcriptome. For these cells with different annotations that mapped to an unexpected germ-layer domain but displayed proximity of transcriptome, we re-annotated the cell identity based on both marker gene expression and prior knowledge of the prospective cell fate gleaned from fate-mapping and lineage tracing studies^25,26 (Supplementary Figs. ⁷ and ⁸). For example, analysis of the heterogeneity of cell types identified two clusters of single cells at Geo-seq position-6P in the posterior epiblast of the E6.75 embryo (Supplementary Fig. ^9a–c). One cluster displayed upregulation of mesoderm genes and enriched gene ontology (GO) terms of mesoderm development, and the other cluster (Cluster 2) displayed a transcription signature of ectoderm. These two clusters were therefore annotated as ‘PS→Mes’ and ‘Epi→Ect’ respectively (Supplementary Fig. ^9c, see methods for nomenclature). RNAscope analyses validated the expression of Mesp1 (Supplementary Fig. ^{9d, e}) and co-existence of Pou3f1-positive cells as well as T-positive cells in the proximal-posterior epiblast of E6.75 embryo (Supplementary Fig. ^{9d, f}). Based on the fractions of re-annotated cell types of every Geo-seq position, a series of Heterogeneity Map was constructed for the gastrula embryo (Fig. 3d), and the presence of all cell types in the three germ layers during gastrulation was identified (Fig. 3d and Supplementary Fig. ⁸).

Inferring the developmental trajectory

The information on cellular heterogeneity further revealed the progressive diversification of cell types in the germ layers (Fig. 3d), which might mirror the events of lineage development during gastrulation. The maps showed that a homogeneous cell population is found in the epiblast and ectoderm domains initially, and the heterogeneous cell types emerge during gastrulation in cell populations in the posterior epiblast/ectoderm, the PS, the mesoderm, and the endoderm. To reveal the developmental connectivity of the cell populations with increasing heterogeneity, the Population Tracing algorithm (Supplementary Fig. ^2f and “Methods”) was applied to infer the putative molecular trajectory of the PS-like cells (PSLCs) in the posterior epiblast of E6.5 embryo to cells in the germ layers at advancing stages of gastrulation (Fig. 3d). Embedded in the transcriptome of the inferred trajectories is the information of the molecular regulatory networks that drive lineage development, for example, the derivation of PS-Mes cells from E6.5 PSLCs, and subsequent allocation of PS-Mes cells to diverse mesodermal cell types at E7.5 (Fig. 3e).

Cells that were originally annotated as “PS” were highly correlated with the PSLCs in our study. It is worth noting that the PSLCs diverge to progenitors of three germ layers at E6.75. A portion of cells in the proximal-posterior epiblast contributes to the epiblast/ectoderm lineage, while the bulk of PSLCs is partitioned into proximal and distal cell groups that were allocated to mesoderm and mesendoderm lineages. The “PS→Mes” cells from proximal PS contribute different mesoderm derivatives at later stages, including the blood and cardiac mesoderm that emerge from proximal PS-mesoderm populations at E7.0. Uniquely, the PSLCs in the E7.25 distal PS contain a cell population that may be the putative MEP, potentially contributing to both mesoderm and endoderm lineages (Fig. 3e and Supplementary Fig. ^{10a, b}). Overall, the inferred molecular trajectories are broadly in line with the prospective fate of the epiblast cells, which provide an entry point for inferring the molecular control of the specification and differentiation of cell types derived from the PS-like cells at successive time points of germ layer development.

Chirality of cell populations in proximal epiblast

An unexpected discovery of the cellular heterogeneity of the epiblast was the finding of different cellular composition of the proximal-lateral epiblast located on contralateral sides (positions 8R2 and 9R2 versus 8L2 and 9L2) of the E7.5 embryo. The populations on the right side are an admixture of cell types, while those on the left side are homogeneous (Fig. 3d, e). To delineate the diversity of cell types in the population, we computationally re-clustered the single cells mapped to position-8R2 and -9R2 and the adjacent positions: 8R1, 9R1 and contralateral positions 8L1, 8L2, 9L1 and 9L2 in the proximal-lateral epiblast, as well as position-9P in the PS and position-9MA and -9MP in the mesoderm. The clustering analysis affirmed the heterogeneity of cell types in positions 8R2 and 9R2 (Fig. 4a–d and Supplementary Fig. ^11a–d)—in particular, a large fraction of the cells at these two positions were allocated to a distinct cluster that was enriched for mesenchyme characteristics¹², in contrast to the remaining cells, which grouped with cells from other positions and were enriched for ectoderm characteristics. The mesenchyme cell type is transcriptionally similar to cells in the mesoderm and is tightly clustered with cells in position-9MA and -9MP (Fig. 4e, f and Supplementary Fig. ^{11e, f}). Applying the Euclidean distance-based Optimization Algorithm, we recapitulated the spatial distribution of single cells within ectoderm section-8/9 (Fig. 4g), where mesenchyme-like cells, distinct from other cell types, were enriched in position-8/9R2.

Open in a separate window

Fig. 4

Heterogeneity of cell types and developmental trajectories of single cells in the proximal-lateral ectoderm of E7.5 embryo.

a–f t-SNE plots showing the heterogeneous cell clusters in position-9R2 and -8R2. Cells are annotated by position (a, c, e) and cell type (b, d, f). Dashed arrows (in a, c, e) denote the heterogeneous cell clusters in position-9R2 and -8R2 versus other lateral positions. Panel a–d showed the single cells of proximal-lateral ectoderm positions in sections 9 (a, b) and 8 (c, d). Panel e, f showed the single cells of proximal-lateral ectoderm, primitive streak, and mesoderm positions in section 9. g The spatial distribution of single cells in the Annulus Model of ectoderm at Geo-seq section 8/9. h–k t-SNE plots (h, j) and the molecular trajectories (i, k) of single cells (imputed using the population tracing algorithm) at position-8R2/9R2 (h, i) and position-8L2/9L2 (j, k) in E7.5–E8.5 embryos. Developmental time points (stage) and cell types (see legend) are indicated in the t-SNE plots. Cell types in E7.75–E8.5 embryos are annotated according to the ‘Gastrulation Atlas’.

To identify the descendants of these cell clusters in position-8R2 and -9R2, we applied the Population Tracing algorithm to infer the descendants of these cells in E7.75–E8.5 embryos (based on single-cell data in the “Gastrulation Atlas”)¹² (Supplementary Fig. ^11g). The cell clusters follow different trajectories, with the mesenchyme cluster giving rise to a multitude of mesodermal tissues, and the ectoderm cluster to the surface ectoderm and rostral neuroectoderm (Fig. 4h, i and Supplementary Fig. ^11h). These two cell clusters at position 8/9R2 are likely to be the progenitors of the lateral mesoderm, and the surface ectoderm and neuroectoderm respectively. In contrast, cells that populate the other six proximal-lateral positions (8/9L1, 8/9L2, and 8/9R1) contribute primarily to ectoderm derivatives and make a minor contribution to ectomesenchyme (presumptively the neural crest cells) (Fig. 4j, k and Supplementary Fig. ^11i–l). Consistent with previous findings¹¹, anterior lateral epiblast (L1, R1) contributes mainly to the neuroectoderm lineage, whereas posterior lateral epiblast (L2) preponderantly contributes descendant to the surface ectoderm. This preponderant contribution of the posterior proximal epiblast on the right side of the embryo to the mesoderm derivatives may implicate an alignment of this developmental event to the specification of left–right (L–R) body asymmetry of the embryo.

Whole-mount in situ hybridization

RNA whole-mount in situ hybridization of the gastrula was performed as previously reported⁴³. Briefly, embryos were fixed in 4% paraformaldehyde (PFA; sigma; #P6148) overnight, rehydrated through the series of 75%, 50%, and 25% methanol at room temperature, washed by 3 times in DPBS, treated with 10mg/mL proteinase K (Life Technologies, # AM2548) in PBS for 8min and post-fixed for 20min in 4% paraformaldehyde. Riboprobes were synthesized by DIG RNA labeling kit and in vitro transcription kit (Roche Applied science, 11277073910; mMESSAGE mMACHINE T7 ULTRA KIT, AM1345; MEGAclear™ Transcription Clean-Up Kit, AM1902) and the probes were amplified using primers listed in Supplementary Data ⁴. Post-fixed embryos were incubated with approximately 1μg/mL of digoxigenin-labeled RNA probe at 68°C overnight, washed and stained with anti-digoxigenin antibody to visualize the expression pattern of the gene transcripts hybridized with the riboprobes.

RNAscope

OCT-embedded E7.5 mouse embryos were cryo-sectioned at 20 μm thickness serially from the distal to proximal region and mounted on electrostatic glass slides. Analysis of gene expression was performed using RNAscope® Multiplex Fluorescent Reagent Kit v2 (Advanced Cell Diagnostics, 323100) following manufacturer’s instructions with minor modifications (proteinase plus treatment for 15min instead of proteinase IV for 30min) using probes supplied by Advanced Cell Diagnostics: mm-Pou3f1 (436421), mm-Mesp1-C3 (436281-C3), mm-T-C3 (423511-C3), mm-Lefty2-C2 (436291-C2), mm-Nodal-C3 (436321-C3), mm-Smad6-C4 (528041-C4), mm-Hand1 (429651). Images were acquired using the Leica TCS SP8 STED system. The fluorescence intensity of each region was calculated using Image J software following standard steps (8-bit format transformation>setting a threshold for background removal>signal measurement). In order to minimize signal variance between embryos and experimental batches, the integrated fluorescent intensity of each region was normalized by calculating the relative ratio to the corresponding region (e.g., relative ratio to MPL region for Lefty2, relative ratio to MPR for Hand1 and Smad6). For each RNAscope probe, at least three biological replicates were examined, and statistical significance was assessed by a student’s t-test.

Small interfering RNA (siRNA) perturbation of BMP signaling

siRNA oligos for Bmp4, Bmpr1a, and Smad1 tagged with Cy3 orange-fluorescent dye and fluorescein amidite (FAM)-labeled negative control siRNA were synthesized by GenePharm. The siRNAs were first dissolved using RNase-free water to make 50μM stock solutions following the manufacturer’s instructions. The efficiency of each oligo for knocking down the transcripts was pre-determined in mouse embryonic stem cells. The siRNA with the highest KD efficiency was selected for the experiments. Equal volumes of the oligos for Bmp4, Bmpr1a, and Smad1 were mixed into one siRNA preparation. For lipofectamine-mediated transfection, 4μL siRNA mixture, 1.5μL LipoFectamine 2000 (Invitrogen, 11668-027), and 4.5μL OptiMEM were pre-mixed and incubated at room temperature for 5min. E7.0 mid-streak stage mouse embryos were collected from the pregnant mice and transferred to drops of PB1 medium placed on the 37°C warm plate under the microscope. Approximate 0.05μL siRNA-lipid mix was micro-injected into the extracellular space between the ectoderm and endoderm of the right side or the left side of the recipient embryo using a flat-tip microinjection pipette⁴⁴. Into control embryos, the same volume of negative control siRNA was injected. The injected embryos were cultured ex vivo as described below.

Ex vivo culture of mouse embryo

Mid-primitive-streak (E7.0) embryos were cultured in the medium of 50% CMRL (Gibco, 11530037) and 50% rat serum, supplemented with 1×Glutamax (Gibco, 35050061), 1×NEAA (Hyclone, SH30238.01), and Glucose (4mg/mL) under 5% CO₂ incubator at 37°C. 1μM chemical inhibitor of BMP signaling (LDN193189, Selleck, S2618-2mg) was added to the culture medium in the respective experimental groups. The cultured embryos were collected at various time points (Figs. 6 and ​and77 and Supplementary Fig. ¹⁵) for whole-mount RNA in situ hybridization or Geo-seq. The whole-mount RNAscope experiments were performed following the published protocol⁴⁵. Images were captured using the LiTone XL system (Light Innovation Technology Limited). Embryos subjected to Geo-seq were fixed in paraformaldehyde for 30min at 4°C, dehydrated in an alcohol series, embedded in OCT compound, and cryo-sectioned serially. Cells were sampled by laser capture microdissection following the published protocol¹⁶. The acquired samples were treated with protease K (Invitrogen, AM2546) at 56°C for 30min, and the lysate was precipitated in ethanol solution and prepared for RNA-sequencing as previously described¹⁶.

Pre-processing of RNA-seq data

The sequencing quality of raw sequencing data was evaluated by FASTQC. Tophat2 v2.0.4 program⁴⁶ was used to map raw reads to mm10 version of the mouse genome with the following parameters: -g 1 -N 4 --read-edit-dist 4 --microexon-search -G annotation.GTF. The mapping ratio was calculated based on the number of mapped reads and total reads for each sample. We calculated fragment per kilobase per million (FPKM) as expression level using Cufflinks v2.0.2 with default parameters⁴⁷. For each embryo, genes with FPKM>1.0 in at least two samples across all samples were retained for further analysis. Finally, the expression levels were transformed to logarithmic space by using the log2 (FPKM+1).

Identification of spatial domain and zipcodes

For E6.5 and E6.75 embryos, zipcodes were identified as follows: (1) applied an adaptive clustering algorithm, Bayesian Information Criterion-Super K means (BIC-SKmeans) algorithm⁴⁸, to identify an optimum number of clusters that can best capture the variance in the data, these optimum clusters were defined as spatial domains, and (2) identify inter-domain DEGs using RankProd⁴⁹ with P value<0.05 and fold change>1.5. The top 50 DEGs of each spatial domain were denoted as zip codes for E6.5 and E6.75 embryos.

For E7.0, E7.25, and E7.5 embryos, zipcodes were identified as follows: (1) Use ComBat⁵⁰ to remove potential batch effects based on the expression of all genes (log₂-transformed) due to variation between samples of different sequencing conditions. (2) Use Python program v2.7.13 to calculate the variance of each expressed gene across all samples and select the top 6000 genes as highly variable genes, then perform principle component analysis (PCA) using FactoMineR⁵¹ package in R. The top 50 highest and 50 lowest PC loading genes from the most significant PCs (E7.0, PC1–4; E7.25, PC1–5; E7.5, PC1–5), which showed inter-domain specific expression patterns, were denoted as zip codes for E7.0, E7.25, and E7.5 embryos respectively. Heatmaps were generated using Cluster 3.0 and JavaTreeView⁵².

Population tracing algorithm

To trace the developmental trajectory of cell populations in different spatial domains, we developed a digital tracing algorithm with the following computational procedures: (1) Use the union of genes in zipcodes of pair of stages of interest as the input gene set, (2) calculate the Euclidean distance of any two domains from two embryos of adjacent stages (for example, the distance between D_A1 and D_B2, D-A1 denotes domain 1 at time point A, D-B2 denotes domain 2 at time point B) (formula-1),

(3) Calculate the mutual nearest neighbors (minimum matrix distance and variation less than 10%, the formula-2 below) for each domain from one developmental timepoint against the domains from the next developmental timepoint and connect the spatial domains. This procedure was repeated across the timepoints from the beginning to the end of the developmental series.

and (4) Convert the distance of any two connected spatial domains to logarithmic space by using log₂ transformation and visualize the genealogy of cell populations by Sankey plot using Google Charts¹¹.

3D Modeling for embryo structure and visualization of Geo-seq data

To mimic the spatial structure and visualize the spatial pattern of gene expression in the embryo, a geometric concentric-oval model that mirrored the architecture of the cup-shaped gastrula-stage mouse embryo was developed. The RNA-seq data of cell samples were assigned to the positions defined by spatial coordinates on the ‘3D corn plot’ model to depict the spatial pattern of expression of the gene or gene group of interest, with the expression levels indicated by a color scale computed from the transcript counts in the RNA-seq dataset.

Multi-dimension single-cell mapping (MDSC mapping)

Based on the 3D model and the zipcode signature, a mathematical algorithm was developed for imputing the location of single cells in the germ layers of the mouse embryo. The mathematical operation included: (1) Calculate Spearman’s rank correlation coefficient (SRCC): The SRCCs between the expression values of the zip codes of each single cell and all samples of the reference embryo were computed to generate, for example, 74 SRCC values for each single cell against 74 Geo-seq samples of E7.0 embryo, (2) Apply a spatial smoothing algorithm to determine the high-confidence location of each cell. This mapping method contrasts with the previous imputation method¹⁵, in which single cells are mapped to the position of the maximum SRCC value. While the higher SRCC may indicate a strong probability of matching to a position, there were cases where a single cell could match to several adjoining positions. Therefore, for mapping the single cells, the top 3 matching positions with maximum SRCCs were extracted, then applied the formula-3 below to calculate a 3D coordinate:

The distance between this location and every Geo-seq sample position was then calculated, and the sample with minimum distance was determined as the best-mapped position of the cell.

Verification of the mapping efficiency

The efficiency of the MDSC Mapping pipeline was evaluated by mapping single cells manually isolated from known positions in E6.5, E6.75, E7.0, E7.25, and E7.5 C57BL/6J embryos. First, embryos were dissected from the decidua in 10% FBS-DMEM medium, and cells were isolated from a specific position of the embryo by mouth pipetting guided by microscopy. Altogether, 19 single cells from E6.5, 32 single cells from E6.75, 37 single cells from E7.0, 24 single cells from E7.25, and 6 single cells from E7.5 were subjected to automatic Smart-seq2 amplification and library construction with the Agilent Bravo automatic liquid-handling platform⁵³. Data preprocessing encompassed mapping, quality control, and normalization using the same criteria previously established¹¹. Based on the transcriptome data of these single cells, their position was mapped to the reference embryo by MDSC Mapping.

To access the accuracy of the mapping methodology, we performed the following analysis: (1) Apply image processing to identify the spatial positions of single cells isolated by pipetting, (2) Apply Gaussian Distribution to mathematically simulate the Confidence Intervals. For E6.5–E7.25 embryos, we employ the standard normal distribution (formula-4) to simulate the confidence intervals:

At E7.5, since the cell number is fewer, the standard normal distribution is adjusted to a more concentrated distribution: f(x)≈1. (3) Calculate the Pearson Correlation Coefficients (PCC) between the Confidence Intervals and MDSC Mapping results. Through this modeling protocol, we showed that the single cells could be mapped at significant fidelity to their site of origin; the PCC of MDSC mapping per stage is 0.7399 (E6.5), 0.8893 (E6.75), 0.7943 (E7.0), 0.7602 (E7.25), and 0.9738 (E7.5).

Single-cell datasets for positional mapping

The 10X Genomics single-cell data were downloaded as raw files from the Gastrulation Atlas: https://github.com/MarioniLab/EmbryoTimecourse2018. Steps of quality control, normalization, batch correction, and clustering were performed using the same criteria as previously described¹².

3D modeling for single-cell resolution map

To visualize the spatial distribution pattern of single cells, an Annulus Model was developed to reconstruct the embryo spatial structure in single-cell resolution. The model comprises (1) Concentric annuli in the anatomical section of the embryo, from the inside outward, representing the epiblast/ectoderm, mesoderm, and endoderm germ layer, respectively. (2) Division of the annulus into interior spaces matching the Geo-seq defined position. Single cells that mapped to a Geo-seq position were distributed uniformly across the interior space of the position in each annulus section. (3) Within each interior space, apply the Bubble Sort Algorithm to align the cells along the known gradient of gene expression level or signaling intensity. For example, on the basis that Bmp4 expression in the posterior epiblast streak decreases along the proximal-distal axis of the embryo, the algorithm code (see below) was applied to rearrange cells in domain 8P in the E6.5 epiblast (x-cd, y-cd, and z-cd are abbreviations of 3D spatial coordinate (x, y, z)), based on the gradient of Bmp4 expression:

for i in range len(8P mapped SCs):

for j in range len(8P mapped SCs)-i-1:

if z-cd[j] <z-cd[j+1] and Bmp4[j] > Bmp4[j+1]:

x-cd[j], x-cd[j+1] = x-cd[j+1], x-cd[j]

y-cd[j], y-cd[j+1] = y-cd[j+1], y-cd[j]

z-cd[j], z-cd[j+1] = z-cd[j+1], z-cd[j]

This imputation enabled the assignment of spatial coordinates specific for each single cell and (4) the construct of the 3D positional map for all single cells for the visualization of the spatial pattern of the single cells of different cell/tissue lineages that are individually defined by transcriptome features. Further information can be drawn from the 3D embryo map for the position of single cells displaying different levels of expression (indicated by the color scale computed from the transcript counts in the “Gastrulation Atlas”) of the gene/gene group of interest.

The migrating mesoderm model

To model the mesoderm migration during gastrulation, we refined the Annulus Model by devising a Gradually Extended Annulus Model for the mesoderm of E7.0 and E7.25 embryos. From distal to proximal, the incomplete annuli that represent the anatomical section of the mesoderm progressively envelop the epiblast from posterior to anterior. At E7.25, the annuli of the proximal mesoderm (7–12M) completely encircles the epiblast. The bubble Sort Algorithm was then applied to assign coordinates to each single cell within the interior domain space of each annulus of the mesoderm.

Euclidean-distance derived optimal coordinates

To refine the spatial distribution pattern of single cells within each Geo-seq position, an optimization algorithm based on Euclidean distance was developed. The mathematical operations include: (1) Compute the gene-expression matrix (GEM) of single cells mapped to a specific position, and transform the matrix to logarithmic space using log₂((normalized count)+1). (2) Based on the logarithmic transformed matrix, calculate the Euclidean distance between every two cells to generate the Euclidean distance matrix (EDM). In order to adapt EDM to the interior space of the Annulus Model, the EDM was normalized through the formula-5 below to make the maximum matrix distance (d_max) equivalent to the maximum length (length_max) within the interior space.

(4) Apply the least square method (formula-6) under the constraint conditions of the spatial mathematical model of the interior space to derive the optimized coordinate for each single cell:

(This example shows the constraint conditions of the mathematical model for position 8P at the E6.5 stage) and (5) visualize the spatial pattern of single cells of different types in each position. All the spatial algorithms and 3D modeling were operated using MATLAB.

Heterogeneity analysis

To reveal the composition of the single-cell population mapped to a Geo-seq position, clustering analysis by t-SNE in Seurat¹⁹ was performed on the single cells for every position in the E6.5–-E7.5 embryos. The cell clusters were annotated on the basis of marker gene expression and the knowledge of the prospective fate of cells in specific regions of the embryo, gleaned from lineage tracing and fate mapping studies^25,26. Based on the composition of annotated cell types in the single-cell population, a Heterogeneity Map was constructed for all Geo-seq positions in the gastrula-stage embryos, with the cell-type composition displayed as pie charts in the corn plots. To construct the molecular trajectory of specific cell types across the developmental stages, single cells of the same type or aligned with a specific lineage were grouped. The Population Tracing algorithm was applied using the averaged gene expression level data to infer the trajectories, which was visualized by the Sankey plot using Google Charts¹¹.

Nomenclature for cell type annotation

The nomenclature ‘X→Y’ and ‘X(Y)’ represent different cell states. X represents the germ layer information of the cell population. ‘X→Y’ indicates these cells are representing a transitional cell state from X to Y. And ‘X(Y)’ represents the precursor of a specified cell type (Y) in the germ layer X.

Clustering and deconvolution of E7.5 mesoderm single cells

Transcriptome data of E7.5 mesoderm single cells were extracted from the ‘Gastrulation Atlas.’ The Seurat package¹⁹ was used to perform single-cell clustering analysis, and t-SNE was applied to visualize the results. To evaluate the representation of identified cell clusters in Geo-seq samples, we prepared cell-cluster labels for each single cell, the single-cell gene-expression matrix, and all Geo-seq samples of mesoderm from the left and right side of the E7.5 embryo. CIBERSORT⁵⁴ was applied to perform cell-cluster deconvolution analysis with default parameters (without quartile normalization). The proportion of each cell cluster in the population was visualized on a t-SNE plot.

Functional enrichment analysis

Functional enrichment of gene sets with different expression patterns was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID)⁵⁵ version 6.8.

Signaling pathway enrichment analysis

In addition to BMP, FGF, Nodal, WNT, Notch, and Hippo-Yap signaling pathways¹¹, Hedgehog, and TGFβ pathways were included in the analysis⁵⁶. Potential signaling-target genes of each pathway were identified by comparing control samples with treatment samples using RankProd (P<0.01 for Hippo-Yap and P<0.001 for others) from published perturbation data (Gene Expression Omnibus accession numbers GSE48092, GSE41260, GSE17879, GSE69669, GSE15268, GSE31544, GSE58664, and GSE90567). Fisher’s exact test, followed by Benjamin–Hochberg correction, was applied to determine the significance of the overlap of the target genes of signaling pathways in different DEG groups.

Generation of Pkd1l1, Dand5 (Cerl2), and Dnah11 (iv) mutant embryos

For the Pkd1l1 rks mutant, substitution mutation in the Pkd1l1 gene³³ was generated by adenosine base editing. Cerl2 mutation was generated by introducing a stop codon in the first exon, which phenocopied the Exon 2 deletion mutations³⁴. For the iv mutant, the first P domain in the Dnah11 gene was deleted³⁶.

Small guided RNAs targeted to the respective genomic regions were designed by using the online tool Chop-chop (http://chopchop.cbu.uib.no/). The DNA fragments containing the T7 promoter and sgRNA sequence and scaffold were transcribed in vitro using the MEGAshortscript Kit (Invitrogen, AM1354), and the products were purified using the MEGAclear kit (Invitrogen, AM1908). DNA fragments for T7-CRISPR-Cas9⁵⁷, T7-ABE 7.10⁵⁸, and T7-YE1-BE4max⁵⁹ were transcribed in vitro using MMESSAGE MMACHINE T7 Ultra Kit (Invitrogen, AM1345) and MEGAclear kit (Invitrogen, AM1908).

C57BL/6 female mice (4 weeks old) were superovulated and mated with the male C57BL/6 mice. Twenty-four hours later, fertilized embryos were collected from oviducts. T7-CRISPR-Cas9 (for Dnah11 mutant), T7-ABE 7.10 (for Pkd1l1 rks mutant), T7-YE1-BE4max (for Cerl2 mutant) mRNAs (100ng/µl), and corresponding sgRNA (100ng/µl) were mixed in HEPES-CZB medium containing 5μg/ml cytochalasin B (CB) and injected into the cytoplasm of fertilized eggs using a FemtoJet microinjector (Eppendorf) with constant flow settings. The injected embryos were cultured in KSOM with amino acids at 37°C under 5% CO₂ in the air to reach the 2-cell stage after 24h in vitro. Two-cell embryos were transferred into pseudo-pregnant ICR female mice, and embryos were collected at E7.5 and E8.5 for analyses.

The authors would like to thank F. Tang for critical discussions. This work was supported in part by the National Key Basic Research and Development Program of China (2019YFA0801402, 2018YFA0107200, 2018YFA0801402, 2018YFA0800100, 2018YFA0108000, and 2017YFA0102700), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA16020501 and XDA16020404), National Natural Science Foundation of China (31630043, 31900573, 31900454, 31871456, and 32130030), and China Postdoctoral Science Foundation Grant (2018M642106). P.P.L.T. was supported by the National Health and Medical Research Council of Australia (Research Fellowship grant 1110751).