MainHuman cardiogenesis requires coordinated interactions among diverse cell types that originate from the cardiogenic mesoderm and neural crest. These populations orchestrate key processes such as formation of the sinoatrial node (SAN), the primary pacemaker of the heart, and maturation of cardiomyocytes1,2. These mechanisms depend on precise spatial and temporal organization of cell states and their integration into multicellular niches that guide morphogenesis and functional maturation. Although genetic, epigenetic and environmental perturbations can disrupt these processes, which can lead to congenital heart disease (CHD), the underlying cellular and spatial mechanisms remain incompletely understood. Notably, the markedly increased incidence of CHD in trisomy 21 (T21) highlights the sensitivity of cardiac development to altered gene dosage3. Despite the identification of numerous CHD-associated genes, a comprehensive view of how early molecular and cellular perturbations emerge and are organized during human heart development remains lacking.Although cardiomyocytes derived from human pluripotent stem cells have provided insights into aspects of cardiomyocyte identity and maturation4,5,6, they do not fully recapitulate in vivo development, particularly with respect to tissue organization and cellular diversity. Recent efforts to profile the fetal heart using transcriptomic and spatial omics approaches have revealed spatial gene expression patterns and cell-type heterogeneity7,8,9,10. However, these studies have generally been limited in developmental coverage, spatial resolution or modality of data types. Furthermore, robust methodologies for the identification and annotation of cellular niches in the developing heart, using various resolutions of spatial omics data, have yet to be established. Consequently, there remains a clear need for comprehensive, multimodal and spatially resolved atlases that can be used to analyse the cellular and molecular dynamics of human cardiac development.Here to address these knowledge gaps, we construct a spatially resolved multimodal atlas of the developing human heart from 4 to 20 post-conceptional weeks (PCW), integrating single-cell transcriptomic data11 with paired single-nucleus ATAC with sequencing (snATAC–seq) and spatial transcriptomics. We delineate 21 tissue niches and introduce TissueTypist, a hierarchical classifier for the annotation of niches across spatial resolutions, benchmarked against histology and independent datasets. Common coordinate systems capture ventricular transmural organization using OrganAxis12 and the spatial progression of SAN pacemaker subtypes and parasympathetic inputs. Finally, comparisons of age-matched euploid and T21 hearts reveal depletion of compact ventricular cardiomyocytes (vCMs) and increased apoptosis in T21 hearts. This result is supported by experiments using an isogenic model of T21 cardiocytes derived from induced pluripotent stem (iPS) cells.Overall, this study offers a foundational resource and analytical framework for decoding human cardiac development and its perturbations in congenital disorders.Multiomic map of human heart developmentWe profiled euploid fetal hearts spanning 4–20 PCW by single-cell RNA sequencing (scRNA-seq)11,13, paired single-nucleus RNA sequencing (snRNA-seq) and ATAC–seq (Multiome, 10x Genomics) and by spatial transcriptomics (Fig. 1a, Extended Data Fig. 1a, Supplementary Tables 1 and 2 and Supplementary Fig. 1). After quality control, 297,473 cells and nuclei were retained. Spatial data comprised 32 capture areas from 10 hearts, including Visium standard definition (SD; 55 μm), Visium high definition (HD; 2 μm) and Xenium datasets (Fig. 1a).Fig. 1: Multiomic atlas of the developing human heart.Full size imagea, Schematic of the study. scRNA-seq, snRNA-seq11,13, single-nucleus Multiome and spatial transcriptomic datasets were generated from 26 euploid human fetal heart samples spanning 4–20 PCW and integrated with T21 heart datasets from 11–14 PCW for comparative analyses. An unnumbered black circle denotes one donor, whereas a numbered circle indicates the donor count at that time point. b, Concordance between transcriptomic and chromatin accessibility profiles across fine-grained cell types. The heatmap shows pairwise Euclidean distances between cell-type median profiles in the gene expression space (top triangle, red) and the gene accessibility score space (bottom triangle, blue). Edge colours indicate parent mid-grained cell-type categories. Concordance was assessed using a two-sided Mantel test with Spearman’s correlation and 999 permutations (r = 0.76, P = 0.01). c, Schematic of the spatial analysis workflow for Visium SD and Visium HD datasets. Visium SD spots were deconvoluted, whereas Visium HD 2-µm bins were aggregated into putative cells. d, UMAP of Visium-HD-derived cells, coloured by major cell type. The dashed region indicates aCMs used for downstream pacemaker cell analyses. e, UMAP of SHOX2+ pacemaker cells from the sinoatrial (SA) region, which revealed sinus horn, SAN head and SAN tail subpopulations with representative marker genes (Extended Data Fig. 2d,e). f, Spatial localization of pacemaker subpopulations. Insets show SA regions and organization along the inflow tract and superior vena cava (SVC). g, Estimated abundance of compact and trabeculated cardiomyocytes in the left ventricle of a representative Visium SD section. Comparable compact and trabeculated cardiomyocyte organization was observed in Visium SD datasets from four independent euploid donors. h, Schematic of developmental left ventricular (LV) structures mapped onto a common transmural coordinate framework (CCF). i, Normalized cell abundance along the left ventricular transmural axis across developmental age. Scale bars, 1 mm (c,f (bottom left),g,h), 0.5 mm (f, top left and middle left) or 0.1 mm (f, top right, middle right and bottom right). CMs, cardiomyocytes; EC, endothelial cell. For all the figures, abbreviations for the fine-grained cell types are listed in Supplementary Table 3. Illustrations in a, c andh created in BioRender; Kanemaru, K. https://biorender.com/j7m5q1g (2026).Cell-type annotations were transferred from our companion study11, which defines a three-level hierarchical system: 6 coarse-grained labels, 14 mid-grained labels and 63 fine-grained labels (Supplementary Table 3). Of the nuclei subjected to multiome sequencing, 167,022 from 11 hearts (spanning 4–20 PCW) also passed ATAC–seq quality control (Extended Data Fig. 1b) with cell-type annotations inherited from their RNA-derived labels. Chromatin-accessible regions (peaks) were identified for each fine-grained cell type, and their union resulted in a total of 508,040 peaks. Overall, 66% of these peaks overlapped—defined as more than 100 base pairs—with the ENCODE candidate cis-regulatory elements (registry v.3)14. We observed strong concordance between ATAC and RNA profiles. In detail, the transcriptionally defined cell types were distinctly separated in the ATAC uniform manifold approximation and projection (UMAP) embedding (Extended Data Fig. 1b). Moreover, a Mantel test comparing inter-cell-type distance matrices between RNA and ATAC modalities revealed a strong correlation (r = 0.76; Fig. 1b). These data provide support for the concordance between transcriptional and chromatin accessibility landscapes.Spatial mapping of cardiac cell typesTo map the spatial distribution of cells, we adopted distinct analytical strategies for Visium SD and high-resolution spatial datasets (Fig. 1c). Cell2location15 was used to map cell-type abundances to Visium SD spots across sections spanning 4–20 PCW. Visium HD data from 5, 6 and 9 PCW heart samples were aggregated from 2-μm bins into cells using bin2cell16. To address mixed-cell signals arising from segmentation or local transcript sharing, we applied reference-based compositional annotation using TACCO17, using matched single-cell RNA-seq and snRNA-seq data and high-confidence spatial singlets (Extended Data Fig. 2a,b), which produced 213,873 Visium HD cells (Fig. 1d). A multimodal morphology-based cell-segmentation algorithm (10x Genomics) and TACCO similarly produced 362,277 Xenium cells from 9 and 16 PCW heart samples (Extended Data Fig. 2c). We annotated fine-grained cell types through marker-driven manual annotation. The identified cell types were consistent with the single-cell and single-nucleus data (Extended Data Fig. 2b).Single-cell analysis of the Visium HD dataset identified approximately 3,300 pacemaker cells in the SAN regions. Subclustering these cells revealed three distinct subpopulations that express markers corresponding to pacemaker cells of the sinus horn (SHPCs), SAN head (SANPC-Hd cells) and SAN tail (SANPC-Tl cells)18 (Fig. 1e and Extended Data Fig. 2d,e). Although human fetal pacemaker cell subtypes have recently been profiled using snRNA-seq19, our dataset is one of the first to map them to their corresponding structures and confirm their annotation (Fig. 1f). In the Xenium-5K dataset, although pacemaker cell subpopulations could not be resolved with the gene panel, SAN pacemaker cells were in close spatial proximity to multiple cell types, including neural cells, coronary cells, fibroblasts and epicardial cells. This result demonstrates the ability of high-resolution spatial transcriptomics to capture fine-scale cellular organization (Extended Data Fig. 2f).Our Visium SD dataset covers a broader range of development than the Visium HD data (Fig. 1a). Spatial mapping of the compact and trabeculated populations of vCMs validated these annotations. That is, they showed consistent overlap with their respective histological structures in the corresponding haematoxylin and eosin (H&E) images (Fig. 1g). To analyse the spatiotemporal transition of compact and trabeculated cardiomyocytes across the ventricular wall, we defined a transmural axis between the manually annotated epicardial and endocardial layers of the left ventricle using OrganAxis12 (Fig. 1h and Extended Data Fig. 2g). This process provided a common spatial coordinate framework to compare features, such as genes and cell-type abundances, across developmental time points and hearts of varying sizes. During development, the abundance distribution of compact cardiomyocytes along this axis shifted from a tight distribution at the epicardium to a more dispersed distribution closer towards the endocardium, whereas an inverse shift was observed for trabeculated cardiomyocytes (Fig. 1i). Moreover, nuclear density was the highest in older hearts, particularly towards the epicardial boundary (Extended Data Fig. 2h). Collectively, these findings indicate that the compact cardiomyocyte population expands in hearts during the second trimester.Cellular niches during developmentHierarchical cellular niche prediction analysis of Visium HD data using CellCharter20 identified 19 spatially coherent cardiac niches (Fig. 2a and Extended Data Fig. 3a). Notably, in the SAN region, fine-grained clustering delineated distinct subdomains corresponding to previously described pacemaker cell populations, including sinus horn-associated, SAN head and SAN tail regions. Beyond the SAN, the framework resolved anatomically distinct niches across multiple regions. In the outflow tract, the ductus arteriosus segregated from adjacent great vessel niches and was characterized by a distinct smooth muscle cell population. In valvular regions, a dedicated valve niche was identified that was marked by colocalization of valve interstitial cells with endocardial cells. In the ventricular myocardium, a layered organization was resolved, which distinguished compact and trabeculated cardiomyocytes and ventricular conduction system components, with patterns consistent across developmental stages. Together, these results demonstrate that the niche-based framework captures both regional specialization and hierarchical tissue organization in the developing human heart (Fig. 2a). To assess the biological validity of these computationally defined niches, we compared their spatial distribution with matched histological annotations derived from paired H&E images (Extended Data Fig. 3b,c). This comparison demonstrated strong concordance between CellCharter-defined niches and histologically defined anatomical regions, thereby supporting its ability to capture physiologically relevant tissue architecture.Fig. 2: Spatial transcriptomics reveals multicellular niches.Full size imagea, Top, schematic of high-resolution niche annotation using Visium HD. Bin2cell-segmented cells were analysed using CellCharter with spatial neighbourhood information, followed by clustering and within-cluster re-clustering to identify 19 cellular niches. Bottom, representative regions from 6 and 9 PCW heart samples are shown with matched H&E images, inferred niches and cell-type compositions. The number of independent donors in which each niche was identified is shown in c. b, Niche annotation in Visium SD. Cell-type abundance was estimated per spot using deconvolution, followed by clustering to define 18 cellular niches. Representative anatomical regions are shown with H&E images, niche annotations and inferred cell-type distributions. The number of independent donors in which each niche was identified is shown in c. c, Hierarchical organization of 7 coarse-grained and 21 fine-grained cardiac niches. Colours indicate fine-grained niche labels and are used consistently across panels. The absence of a modality circle indicates that the corresponding niche was not identified in that modality. Numbers in circles indicate the number of independent donors in which the corresponding niche was identified. Numbers are shown at the finest annotation level resolved for each modality. d, Schematic of the TissueTypist framework. Logistic regression models predict coarse-grained and fine-grained anatomical niches using gene expression from each spatial unit, neighbouring spatial units and distance to the tissue edge. e, TissueTypist prediction in an independent MERFISH dataset9. Predicted niches are shown with expression of representative SA region markers. Gene expression is shown as normalized counts. f, TissueTypist prediction in Xenium data. Predicted niches are compared with independently annotated cell-type localizations, which showed consistency between niche predictions and cellular composition across platforms. Scale bars, 1 mm (a (first column, top), b, e (top left), f (top left)), 0.5 mm (a, first column, bottom) or 0.1 mm (a, second column), e (top right, middle left), f (bottom left and right panels)). AV junction, atrioventricular junction; AVN, atrioventricular node; DA-SMC, ductus arteriosus smooth muscle cell; IFT, inflow tract; ParaN, parasympathetic neuron; PC, pacemaker cell; VCS, ventricular conduction system. Illustrations in b created in BioRender; Kanemaru, K. https://biorender.com/j7m5q1g (2026).For the Visium SD dataset, we performed Leiden clustering on an integrated latent space derived from the spot-by-cell type abundance matrix (Fig. 2b and Extended Data Fig. 4a). Clusters were annotated on the basis of their cell-type composition and histological features identified through expert annotation of the corresponding H&E images. This approach revealed spatially coherent gene-expression signatures that were conserved across developmental stages, thereby supporting the overall consistency of the dataset (Extended Data Fig. 4a). Each annotated tissue structure in the Visium SD dataset showed a distinct combination of cell-type abundances that reflected unique multicellular niches (Extended Data Fig. 4b). LYVE1+ tissue-resident macrophages were enriched in vessel-associated structures and the epicardium, whereas the ventricular compact myocardium showed high abundance of coronary vascular cells, including capillary endothelial cells and pericytes. Notably, the adventitia of the great vessels was enriched for neuron progenitors, consistent with the established increase in neuroattractant signatures in vascular smooth muscle cells11. In line with known developmental timing, parasympathetic neurons were more abundant than sympathetic neurons across cardiac regions, particularly at earlier stages, with neuronal populations predominantly localized to atrial-associated tissues (Extended Data Fig. 4c,d).Together, these observations demonstrate that distinct cellular compositions are organized into spatially coherent multicellular niches across cardiac regions during development. Visium HD resolves fine-grained cellular architecture at near single-cell resolution, whereas Visium SD provides complementary coverage across broader tissue contexts, which enables consistent characterization of cardiac organization across scales.TissueTypist for niche annotationWe developed TissueTypist (https://github.com/Teichlab/TissueTypist), a logistic regression classifier that transfers niche annotations to new spatial datasets. To train a classifier that generalizes across Visium platforms, we unified niche annotations across Visium SD 3 prime, Visium SD formalin-fixed paraffin-embedded (FFPE) and Visium HD FFPE datasets using a hierarchy of 7 coarse-grained and 21 fine-grained niches (Fig. 2c). Samples lacking resolution for fine-grained niche annotation retained an appropriate intermediate label. To bring HD data to the same spatial grain as SD, we aggregated HD cells into pseudobulk windows using a sliding-window approach. The TissueTypist model was trained not only on the intrinsic transcriptomic signatures but also on neighbouring expression and distance to the tissue edge (Fig. 2d). We designed the prediction workflow to be compatible with datasets of various spatial resolution.To assess the robustness of TissueTypist, we performed a systematic cross-validation analysis using a leave-one-section-out framework across all spatial modalities. Across 14 sections spanning Visium SD 3 prime, Visium SD FFPE and Visium HD datasets, the model achieved consistently high performance at the coarse level (weighted F1 = 0.88–0.98), with stable performance at the fine-grained niche level (Extended Data Fig. 5a). Performance varied across individual niches, with per-niche accuracy positively associated with representation in the training data (Pearson’s r = 0.58; Extended Data Fig. 5b). Reduced performance was primarily observed for sparsely represented niches, whereas most well-represented niches were classified with high accuracy. Analysis of the confusion structure revealed that misclassifications were largely confined in the correct coarse anatomical categories, which indicated that errors predominantly reflect ambiguity between closely related substructures rather than incorrect assignment to unrelated tissues (Extended Data Fig. 5c). Notably, increased confusion was observed for boundary-associated niches, including ‘SA node–tail’ and ‘epicardial region’, a result consistent with the observed transcriptional and spatial continuity at tissue interfaces.To evaluate generalizability across spatial platforms and datasets, we applied the trained models to independent high-resolution datasets, including a previously published MERFISH dataset9 and our in-house Xenium-5K data. We observed concordance between TissueTypist predictions and the original MERFISH cellular niche annotations9 (Cohen’s κ = 0.51, accuracy = 0.62, macro F1 = 0.49), with significance confirmed by permutation testing (P 0.1) or region (log fold change > 0, proportion of accessibility > 0.01) being SAN pacemaker cell-specific (compared with other cardiomyocytes) were further selected. For the vCMLs, GRNs that target ‘compact’ and ‘common’ gene signatures were selected similarly to the GRNs of SAN pacemaker cells (compact cardiomyocyte GRNs). In silico TF perturbation analysis was performed using CellOracle32 (v.0.20.0) for the aCM and vCML analysis. A CellOracle object was made using the genes included in the SCENIC+ regulons of each population. Input baseGRN was also prepared using the same regulon (only activator regulons). The cluster labels were generated by combining fine-grained cell types and binned chronological age (n bin = 3). TFs that govern the pruned GRNs were tested for in silico perturbation.Transmural axis analysis using OrganAxisThe spots that belong to left ventricle regions were manually selected (Extended Data Fig. 2g). The inner (endocardial side) and outer (epicardial side) most layer were manually annotated. Using the OrganAxis algorithm12, we calculated the relative distances of each spot from the two manually annotated layers. We calculated the L2 distance of all spots to the closest corresponding points in each annotation by number of k-nearest neighbour (kNN) points (dist2cluster, kNN parameter = 4). To construct an axis anchored to the interface between two structures (‘transmural axis’), we fitted a normalized signoidal curve to the relative distances of each spot and placed spots in relative positions in respect of the two structures (axis_2p_norm). The gene expressions and cell-type abundances were aggregated per bin along the transmural axis and clustered along the transmural axis and visualized (seaborn.clustermap). The cell-type abundance was normalized for each individual sample (age) to account for potential technical variability related to different heart sizes. This normalization enabled us to analyse the transitions of cell-type localizations across the transmural axis. For gene selection, transmural axis was binned into three parts (outer, middle and inner) and DEGs for each part (compared with the rest) were selected (log fold change > 0.5, mean expression > 0.01). Pathway enrichment analysis was performed for the genes upregulated in the outer or inner part. Nucleus segmentation of the H&E images of Visium sections was performed using Cellpose67 implemented in Squidpy68 (squidpy.im.segment, method=cellpose_he, flow_threshold=0.8, channel_cellpose=1).Cell–cell interaction analysisFor the analysis of scRNA-seq and snRNA-seq data, we performed CellPhoneDB25 (v.5) analysis using all cells in the dataset to infer cell–cell interactions among fine-grained cell types: ‘statistical_analysis_method’ with the default parameters and selected ligands and receptors, all of which were expressed in at least 5% of the cells. The ligand–receptor interactions were further selected on the basis of mean expression levels and the biological questions as indicated in the main text and the figure legends. For the analysis using Visium HD data (Fig. 3j), we performed CellPhoneDB (v.5) analysis using cell types and ligand–receptor genes, constrained by their spatial localizations in the identified SAN niche, as determined by TissueTypist. The ‘statistical_analysis_method’ with the default parameters and selected ligands and receptors, all of which were expressed in at least 5% of the cells, were used. For plotting, the ligand–receptor interactions and the cell types were further selected on the biological questions as indicated in the main text and the figure legends.Single-molecule FISHFFPE fetal heart tissue sections (5 μm) were placed onto SuperFrost Plus slides. Staining was performed using a RNAscope Multiplex Fluorescent Reagent kit v.2 assay (Advanced Cell Diagnostics, Bio-Techne), automated with a Leica BOND RX, according to the manufacturers’ instructions. Automated pretreatment included baking at 60 °C for 30 min, dewaxing, heat-induced epitope retrieval at 95 °C for 15 min in buffer ER2 and digestion with Protease III for 15 min. Following probe hybridization, tyramide signal amplification with Opal 520 and Opal 650 (Akoya Biosciences) and TSA-biotin (TSA Plus Biotin kit, Akoya) with streptavidin-conjugated Atto 425 (Sigma Aldrich) was used to develop RNAscope channels. All nuclei were DAPI stained. Stained sections were imaged using a Perkin Elmer Opera Phenix High-Content screening system with a ×20 water-immersion objective (NA of 0.16, 0.299 μm per pixel). The following channels were used: DAPI (excitation, 375 nm; emission, 435-480 nm); Opal 520 (excitation, 488 nm; emission, 500–550 nm); Opal 650 (excitation, 640 nm; emission, 650–760 nm); and Atto 425 (excitation, 425 nm; emission, 463–501 nm).Differential abundance test using MiloTo identify condition-specific shifts in cell populations between euploid and T21 heart samples (or WT and Dp1Tyb mice), we performed neighbourhood-level differential abundance testing using Milo42. Milo detects differential abundance of transcriptionally similar cell populations by first constructing a kNN graph from the single-cell latent space (here derived from scVI integration). Rather than testing predefined clusters, Milo defines a large number of overlapping local neighbourhoods on the kNN graph, each representing a small, transcriptionally coherent region of cell-state space. These neighbourhoods therefore capture subtle and continuous cell-state variation without requiring hard clustering. For each neighbourhood, Milo counts the number of cells originating from each biological sample and performs a negative binomial generalized linear model to test whether neighbourhood abundance differs between conditions while accounting for sample-level replication. P values were corrected for multiple testing using spatial FDR, which accounts for overlap between neighbouring graph regions. This neighbourhood-based framework increases sensitivity to detect condition-specific shifts in cell states, particularly when differences are gradual or confined to subregions of a broader cell type.T21 single-nucleus multiome analysisPost-quality-control T21 snRNA-seq data from four hearts (11–14 PCW) were further curated by removing non-cardiac cells and non-reproducible (across samples) cells. Coarse-grained and fine-grained cell-type annotations were assigned using a label-transfer approach (CellTypist), referencing euploid snRNA-seq data. The T21 dataset was then integrated with age-matched euploid samples (12–14 PCW, 3 hearts), referred to as ‘T21 + Euploid’ in the associated web portal. Using the integrated dataset, the proportions of compact cardiomyocytes in the left and vCM and right vCM populations were quantified. To generate a shared latent space for comparative analysis of T21 and euploid samples, and to assess T21-associated changes in cell population using Milo42, we performed donor-balanced subsampling across cell types (referred to as ‘T21 + Euploid, subsampled for Milo analysis’ in the web portal). Integration of the subsampled data was performed using scVI, accounting for donor identity and continuous covariates, including total transcript counts, %mito and %ribo. The shared latent space generated by scVI was used to generate two-dimensional UMAP embedding (Fig. 5a) and to construct a KNN graph (with n neighbours = 20), which formed the basis for identifying overlapping cellular neighbourhoods using milo.make_nhood. Differential abundance testing was then conducted between T21 and euploid samples, and significantly altered neighbourhoods were identified (spatial FDR