Modeling embryonic heart vascular plexus development and sympathetic innervation on a human heart organoid

TL;DR

Researchers created a human heart organoid that models heart vascular and nerve development, offering a new tool to study heart formation and function.

Contribution

A self-organized human heart organoid that models vascular plexus and sympathetic innervation development.

Findings

EMOs develop a functional coronary-like vascular plexus through VEGF and PDGFβ signaling.

Innervated EMOs show integrated neurovascular organization and respond to nicotine.

The model mirrors in vivo epicardial-to-mesenchymal transition and cell lineage development.

Abstract

Coronary vascularization and sympathetic innervation are tightly coordinated during heart development and are essential for normal cardiac function. Here, we present a self-organized human iPSC-derived epicardium-myocardium organoid (EMO) that mimics heart vascular plexus development and integration with myocardium sympathetic innervation. Through the modulation of VEGF and PDGFβ signaling, EMOs develop a functional, self-generated coronary-like vascular plexus (V-EMOs). and display active epicardial-to-mesenchymal transition trajectories into fibroblast, mural, and vascular cell lineages, mirroring in vivo processes. Assembly of these vascularized EMOs with human iPSC-derived sympathetic neuron spheroids yields innervated EMOs exhibiting integrated neurovascular organization and functional responses to nicotine stimulation (iV-EMOs). This modular, developmentally guided organoid system features coordinated human cardiac, vascular, and neural development, providing a physiologically relevant platform for studying heart morphogenesis, neurovascular interactions, and regenerative therapeutic strategies.

•hPSC-derived EMOs recreate key stages of human heart morphogenesis and architecture•EMOs show a functional self-generated and organized coronary-like vascular plexus•VEGF-PDGFβ regulates epicardial-derived cells specification and spatial patterning•Sympathetic innervated EMOs exhibit physiological responses to nicotine stimulation

hPSC-derived EMOs recreate key stages of human heart morphogenesis and architecture

EMOs show a functional self-generated and organized coronary-like vascular plexus

VEGF-PDGFβ regulates epicardial-derived cells specification and spatial patterning

Sympathetic innervated EMOs exhibit physiological responses to nicotine stimulation

Branco et al. describe a functional vascularized and innervated heart organoid that recreates, in a stepwise manner, the temporal and spatial development of the human heart. The self-generated cues present in this model support epicardial-derived cells specification and the establishment of a self-organized coronary-like vascular plexus, interconnected with a sympathetic neuronal network.

Figures

EMOs stimulation with VEGFA and PDGFBB improves vascularization and induces EMT of epicardial cells(A) Scheme illustrates EMOs generation from hiPSC-derived myocardium and pro-epicardium organoids (MO and PEO), with VEGFA and PDGFBB stimulation from day 6 to day 10 of co-culture.(B) Percentage of EMOs area that stains for the EC marker CD31, with and without VEGFA and PDGFBB supplementation. Data were obtained from representative IF images of at least 4 independent experiments. ∗*p* < 0.05, ∗∗∗∗*p* < 0.0001.(C) Representative IF staining of V-EMOs (3D and organoid slices), highlighting the dev

Single-cell RNA sequencing analysis identifying cell types in V-EMOs and P-EMOs(A) t-SNE plot of filtered cells (*n* = 19,024) from combined V-EMOs and P-EMOs conditions, indicating eight unbiased clusters. Dot plot shows the average expression and percentage of cells expressing canonical markers for each cluster.(B) t-SNE plot (*n* = 501) highlights the compact myocardial-like subpopulation within the cardiomyocyte cluster.(C) Violin plots show the expression of compact myocardium markers within cardiomyocyte cluster.(D) t-SNE plot (*n* = 11,229) of combined Epicardium/early EPDC and epicardi

Single-cell RNA sequencing analysis of epicardial and epicardial-derived cells in V-EMOs and P-EMOs(A and B) t-SNE plot of combined Epicardium/early EPDCs, Epicardium-derived mesenchyme/fibroblast-like cells, mural cells, endothelial cells, and proliferative fibroblast clusters (*n* = 16,490) (A), pseudotime trajectories inferred from transcriptional kinetics (B). Prediction data distinguished two different trajectories: the fibroblast trajectory (clusters 0, 1, and 2) and the vascular trajectory (clusters 0, 1, 3, and 4).(C) Dot plot shows the transcriptional kinetics along pseudotime of key

Impact of PDGFBB and VEGF stimuli on EMOs cell composition(A) t-SNE plot shows the subclustering analysis of epicardium-derived mesenchyme/fibroblast-like cells cluster.(B) Violin plots show the expression level of representative genes for each identified subcluster in epicardium-derived mesenchyme/fibroblast-like cells.(C) Relative frequency of each epicardium-derived mesenchyme/fibroblast-like cells subcluster in V-EMO and P-EMO populations.(D) t-SNE plot shows the subclustering analysis of EC cluster.(E and F) Violin plots show the expression level of representative coronary (E) and artery/

Assembly of V-EMOs and sympathetic neuron spheroids generates functional heart-innervated organoids(A) Representative BF images of SNSs and IF staining of replated SNSs (7 days after replating) after 40 days of differentiation, highlighting the expression of the sympathetic neuronal markers TH and PRPH. Scale bars, 100 μm.(B) Representative BF images of assembloids 2 days post-assembly, highlighting SNs migration from the SNS region toward the EMOs region of V-EMOs.(C and D) Representative IF z stack 3D reconstruction images of V-EMO-SNS assembloids 12-day post-assembly, highlighting SNs (Tuj1

Full text

Introduction

Coronary vascularization and innervation of the myocardium are crucial and interlinked processes in embryonic heart development, both essential for maintaining normal heart function. Sympathetic signals regulate heartbeat rate, relaxation time, contractile force, and conduction velocity, while the coronary vascularization ensures adequate oxygen supply to the working myocardium. Dysfunction of either system is the most prevalent cause of heart failure, arrhythmogenesis, and sudden cardiac death. Therefore, the development of physiologically relevant humanized heart models incorporating both coronary vascularization and sympathetic innervation offers a unique opportunity to dissect the disease mechanisms and explore novel therapies in a clinically relevant context.

To date, no in vitro model has successfully recapitulated the developmental process combining coronary vascularization and sympathetic innervation in hPSC-derived heart organoids. While several studies have described hPSC-derived vascularized cardiac models,1^,^2^,^3^,^4 none have reproduced the spatiotemporal sequence of early coronary vascular plexus development, nor its interactions with the sympathetic ganglionated network. Similarly, the establishment of hPSC-derived heart organoids that recreate myocardium innervation alone has not been achieved, with only a few studies describing 2D monolayer co-culture systems of primary or hPSC-derived ventricular cardiomyocytes (CMs) and sympathetic neurons (SNs).5^,^6^,^7^,^8^,^9^,^10^,^11^,^12^,^13^,^14

In this work, we adopt an embryonic development-inspired in-vivo approach, building on our previously established epicardium-myocardium organoid (EMO) model,15 as the starting point to develop a stepwise platform that mimics the coordinated development of coronary vascularization and SNs innervation in the heart. In the described platform, we assess the impact of VEGF and PDGFβ signaling modulation on epicardial-derived cells specification and spatial organization, and on the development of a functional self-generated and organized coronary-like vascular plexus. Furthermore, through assembly with hiPSC-derived sympathetic neuron spheroids, we established a vascularized-innervated EMO (iV-EMO), which exhibits a physiological response to nicotine stimulation.

Results

Previous work developed by our group has reported a hiPSC-derived epicardium-myocardium organoid (EMO) model (control EMOs), featuring a self-organized epicardium region surrounding the entire surface of a ventricle myocardium-like area15 (Figures S1A–S1C). Knowing that throughout embryonic heart development, epicardium formation is followed by the coronary vascularization and innervation of the myocardium, we hypothesized that EMOs could be the ideal scaffold to replicate and explore the synchronized developmental process of coronary vascularization and innervation of the human heart.

Impact of vascular endothelial growth factor and platelet-derived growth factor-β signal activation on epicardium-myocardium organoid development

To induce vascularization, we started by assessing EMO’s response to the modulation of the vascular endothelial growth factor (VEGF) and platelet-derived growth factor-β (PDGF-β) signaling pathways. PDGF-β signaling is required for efficient epicardium cell migration and derivation into coronary vascular smooth muscle cells in vivo16^,^17 and, accordingly, has been used to induce vascular smooth muscle cells specification from hiPSC-derived epicardial cells in vitro.18 VEGF is a potent angiogenic factor commonly used to induce endothelial cell (ECs) specification from hiPSCs and to support their maintenance. Therefore, we tested the impact of the supplementation of VEGFA and PDGFBB growth factors during the EMOs generation (Figure 1A).Figure 1EMOs stimulation with VEGFA and PDGFBB improves vascularization and induces EMT of epicardial cells(A) Scheme illustrates EMOs generation from hiPSC-derived myocardium and pro-epicardium organoids (MO and PEO), with VEGFA and PDGFBB stimulation from day 6 to day 10 of co-culture.(B) Percentage of EMOs area that stains for the EC marker CD31, with and without VEGFA and PDGFBB supplementation. Data were obtained from representative IF images of at least 4 independent experiments. ∗p < 0.05, ∗∗∗∗p < 0.0001.(C) Representative IF staining of V-EMOs (3D and organoid slices), highlighting the development of a CD31^+^ network of vascular cells surrounding the myocardium region.(D and E) Representative IF staining of V-EMOs (organoid slices), highlighting the sprouting of CD31^+^ cells toward the myocardium region (B), and the co-staining of the CD31 network with the ECM laminin and the pericyte marker NG2 (C).(F) Representative IF staining of the epicardium region in P-EMOs (organoid slices), highlighting the presence of two distinct layers in the epicardium region, an external mesothelium-like layer and a sub-epicardium space adjacent to the myocardial region (NKX2.5^+^).(G and H) Representative IF staining of P-EMOs, highlighting the vascular network of CD31^+^ cells and the enriched staining for VIM (E), NG2, fibronectin (FIBRO), collagen I (COL1), and laminin (LAM) (F), within the compact myocardium region. Scale bars, 100 μm. See also Figures S1 and S2.

Immunofluorescence (IF) staining of EMOs generated with VEGFA (V-EMOs) and PDGFBB (P-EMOs) supplementation revealed a significant increase in CD31^+^ ECs in both conditions compared with control EMOs (Figures 1B and S1D). However, this increase was significantly higher in V-EMOs, increasing from 1.6 ± 0.5% in control EMOs to 3.9 ± 0.2% and 2.9 ± 0.3% of the total EMOs’ area, in V-EMOS and P-EMOs, respectively (Figure 1B). In V-EMOs, we observed a robust vascular network of connected tubules at the interface of epicardial-myocardial areas (Figures 1C and S2A, and Video S1) with branches sprouting toward the myocardium zone (Figures 1D and S2B). The EC-network was also lined by NG2^+^ pericytes and the ECM protein laminin (LAM), essential to form a functional vascular plexus (Figures 1E and S2C). To validate that the vascular network has a perfusable structure, we exposed V-EMOs to fluorescently labeled dextran particles and observed the overlap of the fluorescent signal with the vascular network (Figure S2D), reinforcing the presence of a functional vascular plexus. Supporting also the presence of a functional vasculature in V-EMOs, when exposed to hypoxia injury stimuli, we observed vascular network remodeling with increased vascularization (Figure S2E). In comparison with control, P-EMOs show an increased diameter (696 ± 2 μm EMOs vs. 762 ± 4 μm P-EMOs) (Figures S2F–S2G), with an apparent impact on the epicardial region, where we clearly identified two distinct zones (Figure 1F), similar to what is observed in vivo, in which a mesothelial-cell layer lines a subepicardial region.19

Video S1. Coronary vascularization of V-EMOs, related to Figure 1

Interestingly, although in both V-EMOs and P-EMOs it was clearly identified a compact myocardial-like region adjacent to the epicardial zone, comprising more proliferative CMs (NKX2.5^+^/Ki-67^+^) (Figure S2H), we identified a significant increase in the percentage of these CMs only in V-EMOs (from 8.8 ± 1.0% of NKX2.5^+^/Ki-67^+^within the myocardium region in control EMOs to 19.7 ± 2.0% in V-EMOs) (Figure S2I). In P-EMOs it was observed a stronger staining for the NG2 marker, the mesenchymal/cardiac fibroblast marker vimentin (VIM), and the ECM proteins fibronectin, collagen IV and laminin, in the compact region of the myocardium (Figures 1G and 1H), suggesting a more active epithelial-mesenchymal transition (EMT) in P-EMOs, compared to V-EMOs (Figures S2J and S2K).

Collectively, these observations indicate that EMO exposure to VEGFA stimuli supports the establishment of a functional vascular network and suggest also an impact of PDGFBB signal on promoting 1) the development of an epicardial region with an in vivo-like structural organization and 2) EMT of epicardial-derived cells (EPDCs).

Epicardium-myocardium organoids recapitulate the in vivo developmental fate of epicardium

To decipher the different cell components present in V-EMOs and P-EMOs, we performed a whole-transcriptome analysis of D10CC V-EMOs and P-EMOs, using single-cell RNA sequencing (scRNA-seq), and inferred cellular trajectories based on transcriptional dynamics.

From this analysis, we identified eight distinct clusters present in both V-EMOs and P-EMOs (Figures 2A, S3A, and S3B; Table S3). These clusters were broadly defined as: (1) ventricle CMs, identified by the TNNT2 gene; (2) ECs, expressing the pan-EC markers PECAM1 and CDH5; (3) epicardium/early EPDCs marked by WT1, TBX18, and SEMA3D genes; (4) epicardium-derived mesenchyme/fibroblast-like cells, expressing both EMT and fibroblast markers (POSTN, PTN, DCN, and TCF21); (5) proliferative fibroblast-like cells, expressing the TOP2A gene; (6) mural cells defined by RGS5, MGP, ACTN1 and PDGFRB; and (7–8) two endodermal-derived clusters, both expressing the TTR gene, representing endoderm-derived epithelial cells and hepatic-derived cells.Figure 2. Single-cell RNA sequencing analysis identifying cell types in V-EMOs and P-EMOs(A) t-SNE plot of filtered cells (n = 19,024) from combined V-EMOs and P-EMOs conditions, indicating eight unbiased clusters. Dot plot shows the average expression and percentage of cells expressing canonical markers for each cluster.(B) t-SNE plot (n = 501) highlights the compact myocardial-like subpopulation within the cardiomyocyte cluster.(C) Violin plots show the expression of compact myocardium markers within cardiomyocyte cluster.(D) t-SNE plot (n = 11,229) of combined Epicardium/early EPDC and epicardium-derived mesenchyme/fibroblast-like cells clusters, with pseudotime trajectories inferred from transcriptional kinetics.(E) t-SNE plot highlights the subclustering analysis of Epicardium/early EDPCs cluster.(F) Violin plots show the expression of epicardium mesothelium (cluster 1) genes within the epicardium/early EPDCs cluster.(G) Representative dot and t-SNE plots for epicardial mesothelium, early EPDCs, and epicardial mesenchyme markers, highlighting their transcriptional kinetics along pseudotime trajectory.(H) Violin plots show the expression of myogenic (cluster 2) and adipogenic (cluster 3) genes within the epicardium/early EPDCs cluster.(I) Representative IF staining of V-EMOs (organoid slices), showing the expression of the fatty acid synthase marker (green) and the vascular network (red). Scale bars, 100 μm. See also Figure S3.

Consistent with the IF analysis of V- and P-EMOs, which suggested the presence of a compact myocardial layer adjacent to the epicardial region of the organoids, scRNA-seq data further revealed two distinct subpopulations of TNNT2-positive cells within the CM cluster (Figure 2B). One of these subclusters was enriched in genes previously linked with compact myocardium development in human embryos (MYL2, HEY2, FTH1, and FTL)20^,^21 (Figure 2C).

To assess whether EMOs recapitulate the in vivo developmental fate of epicardium, specifically the transition from epicardial cells to an epicardium-derived mesenchymal cell -population, we performed subclustering and cellular trajectory inference on clusters annotated as “Epicardium/early EPDCs” (cluster 0) and “Epicardium-derived mesenchymal/fibroblast-like cells” (cluster 1) (Figure 2D). For this analysis, we examined the expression of key epicardium and EPDC signature genes curated from recently human fetal epicardium and human heart organoids sequencing datasets.20^,^22^,^23 Within the “Epicardium/early EPDCs” cluster (Figure 2E), a subpopulation of cells expressed canonical epicardial mesothelial markers (“Epicardium mesothelium”, cluster 1), including DSP, CDH1, and CKB (Figure 2F). Notably, genes associated with early or uncommitted EPDCs, such as TBX18, DAPK1, FLRT2, and GPC3, demonstrated transient expression, with higher levels in cluster 0 and gradually decreasing along inferred pseudotime toward cluster 1 (Figures 2G and S3C). In contrast, the expression of DCN, PTN, TCF21, IGFBP7, KRT18, KRT8, and PRDX4 progressively increased along pseudotime, consistent with their established roles in epicardium mesenchyme specification (Figures 2G and S3C). We also identified two smaller cell populations expressing markers linked with epicardial-derived adipose tissue formation (PPARG, GPAM, FASN, and ADIPOR2),24 and with CMs (TNNT2) (Figure 2H), suggesting a myocytic differentiation of EPDCs.23 To confirm the presence of epicardial-derived adipose-like cells, we stained EMOs using an anti-fatty acid synthase (FAS) antibody. The IF results clearly revealed the presence of FAS^+^ cells between the epicardial mesothelium layer and the myocardium region (Figure 2I), recapitulating in vivo observations.25

Following the developmental progression of EPDCs into fibroblast and vascular cells, we incorporated the “Proliferative fibroblasts”, “Mural cells”, and “Endothelial cells” clusters into the pseudotime trajectory analysis (Figures 3A and 3B). We identified two different trajectories: (1) a fibroblast projection characterized by proliferative fibroblast-like cells expressing previously assigned markers for this early differentiation stage, including TOP2A, CCNB2, UBE2C20 (Figures 3C and S3D), and (2) a vascular trajectory progressing toward the mural cells and ECs clusters, enriched for genes associated with smooth muscle cell and cardiovascular system development (MGP, RGS5, LAMA4)20 (Figures 3D and S3E).Figure 3. Single-cell RNA sequencing analysis of epicardial and epicardial-derived cells in V-EMOs and P-EMOs(A and B) t-SNE plot of combined Epicardium/early EPDCs, Epicardium-derived mesenchyme/fibroblast-like cells, mural cells, endothelial cells, and proliferative fibroblast clusters (n = 16,490) (A), pseudotime trajectories inferred from transcriptional kinetics (B). Prediction data distinguished two different trajectories: the fibroblast trajectory (clusters 0, 1, and 2) and the vascular trajectory (clusters 0, 1, 3, and 4).(C) Dot plot shows the transcriptional kinetics along pseudotime of key markers through the fibroblast trajectory, and t-SNE plots for representative markers of proliferative fibroblasts.(D) Dot plot shows the transcriptional kinetics along pseudotime trajectory of key markers through the vascular trajectory, and t-SNE plots for representative markers of mural cells and ECs. See also Figure S3.

Collectively, these results suggest that both V-EMOs and P-EMOs recapitulate the epicardium-derived mesenchyme formation and specification into epicardial-derived fibroblasts, vascular cells, and adipose-like cells, mirroring the developmental potential of the epicardium in vivo.

The vascular network observed in epicardium-myocardium organoids mimics the primitive coronary vascular plexus in the developing human heart

To assess the specific impact of PDGFBB and VEGFA stimuli on EMOs development, we performed subclustering analysis of “Epicardium-derived mesenchymal/fibroblast-like cells”, “Mural cells”, and “Endothelial cells”, in both V-EMOs and P-EMOs.

The subclustering analysis of the “Epicardium-derived mesenchyme/fibroblast-like” cluster revealed a heterogeneous population of cells, with three subclusters (clusters 0, 2, and 3) expressing genes associated with fibroblast/smooth muscle-like cells (SEMA3D, ACTA2, IGFBP7, LUM, TGFBI, BGN, UACA, ITGA1, and TPM2)23 and one subcluster (cluster 1) enriched in genes linked to vascular development (VEGFA, EFNB2, ACKR3, EGR1, and PLXNA4)20^,^22 (Figures 4A and 4B). The latter suggests the presence of perivascular/angiogenic fibroblast-like EPDCs involved in the regulation of vessel formation. Interestingly, most of the V-EMO cells were localized within cluster 1, whereas the P-EMO cells were predominantly localized within cluster 0 and 2 (Figures 4C and S4A). These data suggest a more active vascular network development process in V-EMOs, consistent with immunostaining results, and a stronger commitment of EPDCs into fibroblast/mural-like cells in P-EMOs. In line with this, although mural cells were detected -under both conditions, differential expression analysis between P- and V-EMOs revealed that P-EMOs harbor more mature, extracellular matrix-producing mural cells, with an up-regulated expression of key markers, including COL6A6, MGP, OGN, MFAP4, TLN2, and ELN (Figure S4B and Table S4), and in V-EMOs we observed an up-regulation of markers related with functional pericytes, including CXCL12, RGS5, S1PR3 and NG2.Figure 4. Impact of PDGFBB and VEGF stimuli on EMOs cell composition(A) t-SNE plot shows the subclustering analysis of epicardium-derived mesenchyme/fibroblast-like cells cluster.(B) Violin plots show the expression level of representative genes for each identified subcluster in epicardium-derived mesenchyme/fibroblast-like cells.(C) Relative frequency of each epicardium-derived mesenchyme/fibroblast-like cells subcluster in V-EMO and P-EMO populations.(D) t-SNE plot shows the subclustering analysis of EC cluster.(E and F) Violin plots show the expression level of representative coronary (E) and artery/venous (F) ECs markers.(G) t-SNE plots showing the subclustering analysis of the “Endothelial cells” cluster for each condition (V-EMOs and P-EMOs).(H) Representative IF staining of V-EMOs (organoid slices), highlighting the expression of DACH1, WT1, and NR2F2 in CD31^+^ cells, corroborating a coronary vascular artery-like phenotype for this population.(I) Representative IF staining of P-EMOs (organoid slices), highlighting the generation of two distinct layers in the epicardium region, the mesothelium (WT1^+^/NR2F2^+^/DACH1^+^ cells) and the sub-epicardium space (WT1^−/+low^, NR2F2^+^ cells) (a) and the two distinct sub-populations of ECs, namely CD31^+^/NR2F2^+^ cells within the sub-epicardium space (b and c, ∗1) and in the myocardium surface, and CD31^+^/WT1^+^/NR2F2^-^ cells within the myocardium region (b and d, ∗2). Scale bars, 100 μm. See also Figure S4.

It is known that the primary coronary plexus is mainly derived from both the endocardium and the sinus venosus,19^,^26^,^27 and develops in the subepicardial space. This early coronary EC network then continues to spread and mature at the myocardium surface, forming the coronary veins, and invades the myocardium to form a branching network of coronary arteries.14^,^28 To further confirm the identity of ECs observed in V- and P-EMOs, we performed a detailed characterization of the EC cluster.

When comparing the differentially expressed (DE) markers of our “Endothelial cells” and “Mural cells” clusters with the DE genes assigned to vascular-related clusters from single cell RNAseq data of human embryonic hearts (Figure S4C), we observed a strong correlation, particularly with ECs categorized as coronary vascular cells, supporting the idea that our organoid support the development of a coronary-like early vascular plexus. Additionally, subclustering analysis of EC cluster identified two major EC populations, both expressing markers previously assigned to coronary vasculature development in human hearts (APLNR, FLT4, and ITGA6)20^,^29 (Figures 4D, 4E, and S4D). Interestingly, genes linked to early stages of coronary vascular development and venous-like phenotype30^,^31 were higher expressed in EC subcluster 0 and, on the opposite side, markers linked to coronary arteries specification30 were more strongly expressed in EC subcluster 1 (Figure 4F). When comparing the EC cluster in V- and P-EMOs, we observed that ECs from V-EMOs were mainly located in EC cluster 1, whereas P-EMO ECs were enriched in cluster 0 (Figures 4G and S4E). DE analysis between V- and P-EMOs corroborates this finding, with an up-regulation of genes linked to venous-like ECs (GATA4, FLRT2, and NR2F2)in P-EMOs, while V-EMOs demonstrated up-regulation of coronary artery-like ECs markers (DLL4, CXCR4, EFNB2, and JAG1)30 (Figure S4F and Table S5). It is known that in vivo the immature coronary plexus initially exhibits a venous-like phenotype, which is progressively lost during specification into coronary arteries.28 Consistently, trajectory inference of the EC cluster suggested a temporal transition from cluster 0 to cluster 1 (Figure S4G), indicating that the two clusters may represent different stages of the coronary endothelium development. Together, this data suggests that both V- and P-EMOs develop a coronary-like vascular plexus, with V-EMO enriched in the artery-like ECs and P-EMOs enriched in the venous-like EC phenotype. The latter may either represent primitive coronary veins or reflect a less mature coronary plexus that is still transitioning away from a venous identity before differentiating into coronary arteries.

To confirm the presence of distinct EC subtypes in both conditions, we further characterized the endothelial network observed in V- and P-EMOs by IF analysis. Although the contribution of epicardium to the ECs that comprise the coronary vasculature of the heart remains debatable, studies in mouse and human fetal hearts have shown WT1 expression in both arterial and venous ECs at early development stages.32^,^33 Lineage tracing of pro-epicardium and septum transversum cells also identified WT1^+^/CD31^+^ cells in the coronary vasculature of the developing mouse heart.34 Therefore, we used WT1, together with NR2F2 (COUP-TFII), a coronary venous EC marker,26^,^27 and DACH1, which specifically labels coronary vessels in mouse models,29^,^35^,^36^,^37 for further analysis.

IF analysis confirmed that in V-EMOs the majority of CD31^+^ ECs are positive for both WT1 and DACH1 markers, but not the venous marker NR2F2 (Figures 4H, S4H, and S4I). In contrast, P-EMO contained predominantly CD31^+^/DACH1^+^/NR2F2^+^ ECs within the subepicardial space and on the myocardial surface, with a small subset of CD31^+^/WT1^+^/NR2F2^-^ ECs within the myocardial region (Figures 4I and S4J).

Assembly of vascular plexus epicardium-myocardium organoids and sympathetic neuron spheroids generates functional heart-innervated organoids

To incorporate innervation into vascularized EMOs, we first established a differentiation platform to attain sympathetic neuron spheroids (SNSs) from hiPSCs. To achieve that, we adapted a previously established 2D differentiation protocol38 to a 3D environment. Successful differentiation was confirmed by RT-qPCR analysis, which demonstrated a significant up-regulation of SNs’ markers SOX10, ASCL1 and PHOX2B at the end of the neural induction period (Day 12 of differentiation) (Figure S5A), and by positive immunostaining for the SN enzyme tyrosine hydroxylase (TH) and for the peripheral neuronal marker peripherin (PRPH) after 30–40 days of differentiation (Figure 5A).Figure 5. Assembly of V-EMOs and sympathetic neuron spheroids generates functional heart-innervated organoids(A) Representative BF images of SNSs and IF staining of replated SNSs (7 days after replating) after 40 days of differentiation, highlighting the expression of the sympathetic neuronal markers TH and PRPH. Scale bars, 100 μm.(B) Representative BF images of assembloids 2 days post-assembly, highlighting SNs migration from the SNS region toward the EMOs region of V-EMOs.(C and D) Representative IF z stack 3D reconstruction images of V-EMO-SNS assembloids 12-day post-assembly, highlighting SNs (Tuj1^+^) network in the myocardial region of the assembloids (cTnT^+^) (C), and the interplay with the coronary-like vascular (CD31^+^) network (D). See also Figure S5.

Functional activity of SNs was assessed, using a MEA system, after 50 days of 3D differentiation followed by 8–12 additional days of culture on the electrodes. This analysis confirmed the electrical activity of the neurons, with an average conduction velocity of 0.52 ± 0.03 m/s and a firing rate of 3 ± 0.3Hz (Figures S5B and S5C). SNs also responded to nicotine stimulation, with an increased neuron conduction velocity (11 ± 2%) and an increased number of spikes/burst (58 ± 12%) (Figure S5D). Additionally, we confirmed that SNSs produce the neurotransmitters dopamine (269 ± 78 pg/mL) and norepinephrine (56 ± 4 pg/mL) (Figure S5E), reinforcing their physiological activity. After platform optimization and functional validation of the SNs, we selected V-EMOs, owing to their more robust vascular network, to fuse with fluorescently labeled SNSs (>D50 of differentiation) in 96-well plates, to generate innervated EMOs (iV-EMOs). It was possible to observe neuronal projections toward the V-EMOs region of the assembloids within 48 h (Figure 5B). By days 6 and 12, a dense network of SNs had formed at the surface of the myocardial region (Figure S5F). Additionally, to evaluate the presence of neuro-vascular pattern integration in iV-EMOs and to identify the vascular network, IF analysis of assembloids was performed using beta tubulin III (Tuj1) and CD31, respectively. We observed a well-developed network of interconnected neurons and vascular cells in iV-EMOs (Figures 5C, 5D, S5G, and S5H; Video S2), similarly to what is observed in vivo.

Video S2. Sympathetic innervation of V-EMOs, related to Figure 5

To further confirm the functional interaction between SNs and ventricle CMs, we examined the effect of SN stimulation on the contraction activity of iV-EMOs 2 weeks after initiating the co-culture. These assembloids were treated with 100 μM of nicotine, which activates nicotinic acetylcholine receptors, and the beating rate was assessed after 30 min. We registered a significant increase in beating rate upon nicotine treatment, from 14 ± 4 to 27 ± 5 BPM (Figure S5I). As control, we tested nicotine stimulation on age-matched V-EMOs without external sympathetic innervation, as well as myocardium organoids (MOs). Although, as expected, the contraction of MOs was not affected by nicotine stimulation, in V-EMOs, we observed an impact on organoid contraction. Specifically, we observed that 32 ± 5% of the organoids stopped contraction, which was rescued after changing the culture media, and the remaining suffered an increase in the beat frequency from 13 ± 2 BPM to 18 ± 2 BPM. This observation confirms that the SNs can modulate CM function in iV-EMOs, while also suggesting V-EMOs can be modulated by exogenous stimuli, without external sympathetic innervation, although less robustly. The observed response of V-EMOs to nicotine stimulation can be seen as an in vitro artifact or can be indicative of the presence of an intrinsic cardiac nervous system in these organoids. In fact, the scRNA-seq data identified cells expressing the CALB1 gene within the epicardium mesothelium cluster. Recent work has shown that calbindin 1-positive cells localize to epicardial ganglia, a part of the intrinsic cardiac nervous system, and influence cardiac function independently of the central nervous system.39 Previous studies further suggest that epicardial ganglia control myocardium contraction locally, independently of extrinsic autonomic nerves, and that the nicotinic stimulation of cardiac intrinsic ganglia impacts heart function.40^,^41 This evidence may explain the observed results in our model, although further validation is needed.

Discussion

Normal heart function is ensured by functional vasculature and the autonomic nervous system. Alterations to the structure and function of the coronary vascular network and in the excitability, density, distribution, and neurotransmitter content of heart innervation are well-known events that compromise cardiac function. Therefore, the development of in vitro human heart organoid models that not only reproduce the myocardium but also the complementary cell systems that support the normal function of the heart is crucial for advancing clinical and biomedical applications.

Herein, we describe a physiologically relevant hiPSC-derived heart organoid (EMO) that recreates not only the epicardium-myocardium interaction, but also both the coronary-like embryonic vascular plexus and the sympathetic innervation, featuring an in-vivo-like structural organization and function. Although previously reported hiPSC-derived heart organoids have been shown to partially recapitulate the cardiac vasculature/endocardial vascularization,42^,^43^,^44^,^45 there is no model that reports a self-generated and self-organized coronary-like early vascular network, a limitation reported in recent hiPSC-derived heart organoid models.46 Additionally, these models do not replicate the interplay between vascularization and the establishment of a functional innervated network, whose development remains not completely disclosed. Despite that, it has been shown that the epicardium is an important source of paracrine factors and progenitor cells for the establishment of a functional coronary vasculature in the human heart.38 Moreover, it has been reported that sympathetic ventricle innervation develops in close cross-talk with vascularization, with neurotrophic factors secreted by epicardium-derived vascular cells guiding the extension of axons through the subepicardium and inducing invasion into the myocardium.14^,^47

Our approach leveraged EMOs as a developmentally inspired platform to recreate the in vivo-like generation and interaction between epicardium, coronary vascular cells, and sympathetic neurons. We studied the impact of VEGF and PDGFβ modulation on epicardial-derived cells specification, vascularization development, and spatial organization. scRNA-seq demonstrated that EMOs recapitulate the developmental trajectory of epicardium development and EPDCs specification into fibroblast, mural cells, ECs, and adipose-like cells. Furthermore, our data also revealed that VEGF and PDGFβ signaling pathways activation during EMOs generation supported the generation of a functional and well-developed network of coronary-like endothelium, lined by pericytes and ECM proteins, and potentiated the development of a well-demarked subepicardial space and EPDC specification and differentiation, respectively. These findings highlight the potential of combining VEGF and PDGFβ modulation for fine-tuning epicardial and vascular maturation. Lineage tracing studies in the mouse model provided important insights and established a consensus idea that the primary coronary plexus is mainly derived from both the endocardium and the sinus venosus, and epicardial-derived pericytes are required for coronary vasculature maturation, integrity, growth, and vascular guidance.19 Although the contribution of epicardium to the ECs that comprise the coronary vasculature of the heart is still debatable, it has been demonstrated in mouse and human fetal hearts that WT1 is expressed in ECs of both arteries and veins at early stages of development.32^,^33 In our heart organoid model, we observed that the CD31^+^ ECs, in the sub-epicardial space and myocardial region, express WT1, which is consistent with previously described models.

The positive impact of epicardium on myocardium growth, compaction, and maturation in vivo is well described in the literature.48^,^49 While we previously reported the presence of a compact myocardial region in EMOs,15 the improvement of vascularization significantly increased the robustness of the compact myocardial region in these organoids, showing an increased number of proliferative CMs. It is well documented in mouse models that coronary ECs support myocardial compaction and heart wall expansion, and impairment of these developmental processes leads to congenital heart defects, including ventricular non-compaction cardiomyopathy.50

The control of heart function through autonomic nervous system modulation is a promising therapeutic strategy for several arrhythmogenic cardiomyopathies. As a proof-of-concept, we showed that the assembly of hiPSC-derived SNS with vascularized EMOs promotes the development of a complex vascular interconnected sympathetic neuronal network that exhibits physiological responses to nicotine stimulation, thus proving the functional interaction between the SNs and the CMs present in the iV-EMOs.

In summary, this work presents a multifactorial platform to generate a modular physiologically relevant heart organoid that, although exploring a self-organization environment, is guided by exogenously controlled key steps, which ensure a robust and reproducible platform. Specifically, iV-EMOs are generated in a stepwise and in-vivo-like environments, recreating the sequential stages of human heart embryogenesis, including 1) mesendoderm formation, 2) cardiac progenitor cell specification, 3) ventricular myocardium differentiation, 4) pro-epicardium/septum transversum organ formation, 5) epicardial establishment and EPDCs specification and migration, 6) myocardium maturation and compaction, 7) coronary-like early vasculature development, and 8) sympathetic innervation of the myocardium.

iV-EMOs thus provide a unique platform to address fundamental questions in heart development and pathology. Moreover, the incorporation of the sympathetic autonomic nervous system in EMOs provides an opportunity to explore therapeutic strategies based on SNs activity manipulation/modulation in a human-based setting, which is a promising therapeutic strategy for several arrhythmogenic cardiomyopathies. Additionally, it is well known that altered density and/or function of heart innervation and vascularization are the most prevalent causes of heart failure and sudden cardiac death. In fact, apart from CM death and the subsequent electrophysiological and fibrotic remodeling observed after myocardial infarction events, the level of cardiac sympathetic nerve fiber loss has been considered the most relevant predictor of consequent ventricular arrhythmias.51^,^52 Therefore, iV-EMOs represent a unique model to unravel crucial signals involved in the normal coronary vascular growth and SN innervation of the ventricle myocardium, which can be further explored to develop new regenerative medicine-based therapeutic avenues to restore innervation and promote re-vascularization in the adult heart after ischemic events.

Limitations of the study

In this work, we developed and optimized a human heart organoid model that recapitulates the coordinated development of epicardium, embryonic vascular plexus, and sympathetic innervation. This platform builds upon our previous EMOs, generated through pro-epicardium organoid specification. In vascularized EMOs (V-EMOs), the vascular network is self-generated and self-organizes within the organoid, whereas sympathetic innervation in iV-EMOs is achieved through an assembly-based approach.

Although we demonstrate the formation of a well-developed and interconnected vascular plexus containing ECs and supporting mural cells, further functional and structural characterization of this network is required. Specifically, further assessment of vessel perfusion and quantitative measurement of vessel diameter are needed. Moreover, while our scRNA-seq data indicate the presence of both artery-like and venous-like EC populations, further confirmation using additional subtype-specific markers and immunostaining analyses is necessary to more definitively establish the identity and maturation status of these EC subpopulations. Extended culture periods would also be informative to evaluate the long-term evolution and maturation of the vascular network. The in-vivo-like spatial organization of the vascular network observed in our organoids, together with comparison between our scRNA-seq data and published human embryonic heart datasets, supports a coronary-like phenotype for the vascular network observed in V-EMOs. However, additional validation is required to conclusively demonstrate that this organoid model fully recapitulates embryonic coronary vasculature development. To further dissect the developmental processes occurring in V-EMOs and to clarify the origin of vascular progenitors, scRNA-seq analyses at earlier and sequential stages of differentiation would be highly informative. Finally, this study was conducted using a limited number of hiPSC lines. Validation of the platform across additional genetic backgrounds will be necessary to confirm the robustness and broader applicability of the model.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Mariana A. Branco (marianabranco@accelbio.pt).

Materials availability

All unique/stable reagents generated in this study are available from the lead contact with a complete Materials Transfer Agreement.

Data and code availability

Single-cell RNA-seq data is available through GEO:GSE309829. The raw and processed single-cell RNA-seq data generated in this study have been deposited in the NCBI Sequence Read Archive (SRA:PRJNA1356374). All original code used for data processing and downstream analysis is available on GitHub. Accession numbers and all software resources used in this study are indicated in the key resources table. Microscopy data reported in this article will be shared by the lead contact upon request. Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request.

Acknowledgments

We thank the Gulbenkian Institute for Molecular Medicine (GiMM) for access to the bioimaging facility. We would also like to thank the Flow Cytometry Platform of GiMM for their technical support. We acknowledge funding received from the national foundation 10.13039/501100001871FCT - 10.13039/501100001871Fundação para a Ciência e Tecnologia, Portugal, through 10.13039/501100015620Institute for Bioengineering and Biosciences - 10.13039/501100015620iBB (UIDB/04565/2020 and UIDP/04565/2020), through the Associate Laboratory Institute for Health and Bioeconomy – i4HB (LA/P/0140/2020), and through the project 2022.02251.PTDC granted to M.A.B. We also acknowledge funding received from 10.13039/501100007434Agência Nacional de Inovação (RE-C05-i02 – Missão Interface N.º01/C05-i02/22 and European Union’s Horizon 2022 Research and Innovation programme under the grant agreement NO 101115536.

Author contributions

Conceptualization and methodology, M.A.B and M.M.D.; experimental design, M.A.B; cell culture and analysis, M.A.B, M.M.N., and A.L.R.; sc-RNAseq data processing and analysis, J.M., M.A.B., and M.B.; results discussion, M.A.B; M.M.N., A.L.R., M.F.T., and M.M.D.; writing – original draft, M.A.B; writing – review and editing, M.M.N., A.L.R., M.F.T., M.M.D., J.M., and M.B.; funding acquisition, M.A.B, M.M.D., and J.M.S.C.

Declaration of interests

This article forms the basis of a patent application (PPP 2024200621851) on which M.A.B. and M.M.D are named as inventors.

STAR★Methods

Key resources table

REAGENT or RESOURCESOURCEIDENTIFIERAntibodiesMouse anti-CD31DakoCat# M0823; RRID:AB_2114471Mouse anti-cTnTThermo Fisher ScientificCat# MA5-12960; RRID:AB_11000742Rabbit anti-NKX2.5AbcamCat# ab97355; RRID:AB_10680260Rabbit anti-ISL1AbcamCat# ab178400; RRID:AB_2927537Rabbit anti-WT1AbcamCat# ab89901; RRID:AB_2043201Mouse anti-Tuj1BiolegendCat# 801213; RRID:AB_2313773Rabbit anti- Fatty Acid SyntaseAbcamCat# ab22759; RRID:AB_732316Rabbit anti-Collagen IAbcamCat# ab34710; RRID:AB_731684Mouse anti-FibronectinAbcamCat# ab253288Rabbit anti-LamininAbcamCat# ab11575; RRID:AB_298179Mouse anti-VimentinSigma AldrichCat# V6630; RRID:AB_477627Rabbit anti-Ki-67AbcamCat# ab833; RRID:AB_306483Rabbit anti-DACH1ProteintechCat# 10914-1-AP; RRID:AB_2230330Mouse anti-COUP-TFII/NR2F2R&D SystemsCat# PP-H7147-00; RRID:AB_2155627Mouse anti-THBiolegendCat# 818002; RRID:AB_2734573Rabbit anti-PRPHMilliporeCat# AB1530; RRID:AB_90725Alexa Fluor 488 Donkey Anti-mouse IgGThermo Fisher ScientificCat# A21202; RRID:AB_141607Alexa Fluor 488, Donkey Anti-rabbit IgGThermo Fisher ScientificCat# A21206; RRID:AB_2535792Alexa Fluor 546, Donkey Anti-mouse IgGThermo Fisher ScientificCat# A11030; RRID:AB_2737024Alexa Fluor 546, Donkey Anti-rabbit IgGThermo Fisher ScientificCat# A10040; RRID:AB_2534016Alexa Fluor 647, Donkey Anti-rabbit IgGThermo Fisher ScientificCat# A31573; RRID:AB_2536183Alexa Fluor 647, Donkey Anti-mouse IgGThermo Fisher ScientificCat# A31571; RRID:AB_162542Chemicals, peptides, and recombinant proteinsmTeSR^TM^1StemCell TechnologiesCat# 85850B-27 supplement minus insulinThermo Fisher ScientificCat# A1895601B-27 supplementThermo Fisher ScientificCat# 17504044CHIR99021ReproCELLCat# 04-0004-02IWP-4ReproCELLCat# 25704-0036BMP4Thermo Fisher ScientificCat# 120-05ET-10UGRetinoic AcidSigma AldrichCat# R2625-100MGAdvanced DMEM/F12Thermo Fisher ScientificCat# 12634-010GlutaMAX™ SupplementThermo Fisher ScientificCat# 35050061L-Ascorbic acid 2-phosphate sesquimagnesium salt hydrateSigma AldrichCat# A8960VEGF-165Thermo Fisher ScientificCat# 100-20PDGF-BBR&D SystemsCat# 220BB050LDN-193189ReproCELLCat# 04-0074SB431542Sigma AldrichCat# S4317-5MGNeurobasal™ Plus MediumThermo Fisher ScientificCat# A3582901Advanced DMEM/F-12Thermo Fisher ScientificCat# 12634010RPMI 1640 MediumThermo Fisher ScientificCat# 11875093Essential 6™ MediumThermo Fisher ScientificCat# A1516401N2 supplementThermo Fisher ScientificCat# 1750248B27 supplement plusThermo Fisher ScientificCat# A3582801NGFThermo Fisher ScientificCat# 450-01BDNFThermo Fisher ScientificCat# 450-02-100UGGDNFThermo Fisher ScientificCat# 450-10-50UGFluo-4-AM dyeThermo Fisher ScientificCat# F23981Bacterial and virus strainspAAV-CAG-tdTomatoAddgeneCat# 59462-AAV1Chemicals, peptides, and recombinant proteinsFluorescein isothiocyanate–dextranSigmaCat# FD40Critical commercial assaysNorepinephrine ELISA kitAbnovaCat# KA1891Dopamine ELISA kitsAbnovaCat# KA1887Deposited dataRaw single-cell RNA-seq dataNCBI Sequence Read Archive (SRA)PRJNA1356374Processed single-cell RNA-seq dataGene Expression Omnibus (GEO)[GSE309829](GSE309829)Original code used for data processing and downstream analysisGitHubhttps://github.com/CoLAB-AccelBio/heart_organoids_single_cell_analysis.Experimental models: Cell linesHuman induced pluripotent stem cell line DF6WiCelliPS-DF6-9-9T.BSoftware and algorithmsFiji/ImageJ V2.0Schindelin et al., 2012https://imagej.net/Fiji.htmlPrism 8Graphpad Software Inc.https://www.graphpad.comZEN Blue 2.5Zeisshttps://www.zeiss.com/microscopy/en/products/software/zeiss-zen.htmlCell Ranger v9.0.110x Genomicshttps://www.10xgenomics.com/support/software/cell-ranger/latestscStudioBica et al. 2025 (bioRxiv)https://github.com/DiseaseTranscriptomicsLab/scStudioR (v4.1.2)R Projecthttps://www.r-project.orgR (v4.3.3)R Projecthttps://www.r-project.orgSeurat (v4.1.1)Satija Labhttps://satijalab.org/seurat/scDblFinder (v1.8.0)Bioconductorhttps://bioconductor.org/packages/scDblFinderscran (v1.22.1)Bioconductorhttps://bioconductor.org/packages/scranscater (v1.22.0)Bioconductorhttps://bioconductor.org/packages/scaterTSCAN (v1.44.0)Bioconductorhttps://bioconductor.org/packages/TSCANclustree (v0.5.1)CRANhttps://cran.r-project.org/web/packages/clustree/index.htmldevtools (v2.4.5)CRANhttps://cran.r-project.org/web/packages/devtools/index.htmlggridges (v0.5.6)CRANhttps://cran.r-project.org/web/packages/ggridges/index.htmlggplot2 (v3.5.2)CRANhttps://cloud.r-project.org/web/packages/ggplot2/index.htmlRColorBrewer (v1.1-3)CRANhttps://cran.r-project.org/web/packages/RColorBrewer/index.htmldplyr (v1.1.4)CRANhttps://cran.r-project.org/web/packages/dplyr/index.htmltidyr (v1.3.1)CRANhttps://cran.r-project.org/web/packages/tidyr/index.htmlDT (v0.33)CRANhttps://cran.r-project.org/web/packages/DT/index.html

Experimental model and study participant details

This research uses the hiPSC line iPS-DF6-9-9T.B. provided by WiCell Bank. This cell line is vector free and was reprogramed from foreskin fibroblasts with a karyotype 46, XY that were collected from healthy donors using defined factors in the Laboratory of Dr. James Thomson, at University of Wisconsin.

Method details

Cell maintenance

hiPSCs were maintained in mTeSR^TM^1 (StemCell Technologies) in six-well plates coated with Matrigel™ (Corning). The medium was changed daily. Cells were routinely passaged every three to four days using 0.5 mM EDTA solution (Thermo Fisher Scientific).

Generation of EMOs

The generation of EMOs starts with the combination of hiPSC-myocardium organoids (MOs) and hiPSC-pro-epicardium organoids (PEO), previously described by our group.15^,^53 Briefly, the differentiation of hiPSC-spheroids (4000 cells/spheroid) was initiated on day 0 (D0) by replacing the mTeSR^TM^1 by RPMI 1640 medium (Thermo Fisher Scientific) with 2%(v/v) B-27 minus insulin (Thermo Fisher Scientific) (RB27-) supplemented with 11 μM CHIR99021 (CHIR). After 24 h, the medium was changed and on D3, cells were cultured in RB27- supplemented with IWP-4 at a final concentration of 5 μM, for two days. At D5, in the case of MOs differentiation, the medium was changed to RB27-, and in the case of PEOs differentiation, 3 μM CHIR, 25 ng/mL BMP4, and 4 μM RA were added from D5 to D7, using DMEM/F12 + 2.5 mM Glutamax + 100 μg/mL of ascorbic acid (DMEM/F12) as basal medium. At D7, aggregates were flushed from the AggreWell™800 plate (StemCell) and transferred to 6-Well ultra-Low Attachment (ULA) plates. Thereafter, the medium was changed every two days until cell harvest. In the case of MOs, the medium used was RPMI 1640 medium supplemented with 2%(v/v) B-27 and in the case of PEOs it was used DMEM/F12 medium.

To establish the co-culture system, D11 MOs and D11 PEOs were dissociated to single cells using 0.25% (v/v) Trypsin-EDTA for 7 min at 37 °C. After cell counting, both cell types were combined at a proportion of 90%MOs:10%PEOs cells and reaggregated in ultra-low 96-well microwell plates. Organoids were cultured in DMEM/F12 medium for six days. After that time, the culture medium was supplemented with VEGFA (50 ng/mL) and PDGFBB (5 ng/mL) factors, alone or combined, from D6 to D10 of co-culture, to evaluate the sole and combined effect of these factors to promote epithelial to mesenchymal transition (EMT) of epicardial progenitor cells and coronary vascularization.

Differentiation of sympathetic neuron spheroids (SNSs) from hiPSCs

To initiate differentiation of SNSs, hiPSCs were harvested with accutase (Sigma) for 7 min at 37 °C. After dissociation, spheroids of hiPSCs were generated using microwell plates (AggreWell™800, StemCell Technologies) according to the manufacturer’s instructions. Cells were plated at a cell density of 0.3 × 10^6^ cells/well in mTeSR™1 supplemented with 10 μM ROCKi. 24 h later, the differentiation process (based on11) was initiated. Briefly, from D0 to D11 we used 1) E6 basal media from D0 to D3, 2) mixed E6 and N2B27 (Neurobasal™ Plus Medium + B27 Plus + N2 supplement) at (75:25)%, respectively, from D3 to D6 3), and mixed E6 and N2B27 at (50:50)%, from D6 to D8, and 4) only N2B27 from D8-D11. At D0, basal medium was supplemented with 500 nM LDN193189 (LDN) and 10 μM SB431542 (SB) and at D2, the medium was supplemented with 500 nM LDN + 10 μM SB + 3 μM CHIR + 10 μM DAPT. At D3, the medium was supplemented with 10 μM SB + 3 μM CHIR + 10 μM DAPT + 1 μM SAG and aggregates were transferred to ULA 6-well plates. From D4 to D7, the medium was supplemented every day with 3 μM CHIR + 10 μM DAPT + 1 μM SAG. From D7 to D11, the medium was supplemented every day with only 1 μM SAG. At D11, SN progenitor spheroids were singularized using accutase for 10 min at 37 °C. After cell counting, cells were resuspended in N2 maturation media (Neurobasal™ Plus medium + 2% (v/v) B27Plus + 1% (v/v) Glutamax + 1% (v/v) N2 supplement + 0.2 mM Ascorbic Acid + 10 ng/mL NGF + 10 ng/mL BDNF + 10 ng/mL GDNF) supplemented with 10 μM of ROCKi (100 μL/well), and then seeded on a U-shaped ULA 96-well plate. From this point on, the medium was changed every three-four days and maturation was prolonged for 30-40 days.

Generation of innervated heart assembloids

Heart-innervated assembloids were generated by combining EMOs, obtained after 10 days of co-culture, with SNSs obtained after 30-40 days of maturation in U-shaped ULA 96-well plate (1:1 SNSs/EMOs). These assembloids were cultured in 50:50 N2 maturation medium and DMEM/F12 for 10 days and the medium was changed every 3 days.

Viral labeling of neural spheroids

To follow SN projections towards the EMO region of the assembloids, day 30 SNSs were labeled with pAAV-CAG-tdTomato virus (Addgene, cat. no. 59462-AAV1). For labeling SNSs, each spheroid was cultured in 40 μL of N2 maturation medium supplemented with 0.5 μL of virus solution in a 96-well plate, overnight. On the next day, 100 μL of N2 maturation medium was added to each well. After 24 h, the SNSs were washed twice with PBS and cultured in fresh N2 maturation medium. Fluorescently labeled SNSs were usually observed seven days after virus infection.

Immunostaining analysis of EMOs and SNSs

EMOs and SNSs were fixed in 4% PFA in PBS at 4 °C overnight. After PFA removal, cells were stored in PBS at 4 °C for further analysis.

Slices

EMOs were incubated in 15%(w/v) sucrose in PBS, at 4 °C overnight. Afterwards, the 3D models were embedded in 7.5%(w/v)/15%(w/v) gelatin/sucrose and frozen in isopentane at −80 °C. Organoid/spheroid sections of 10-12 μm of thickness were cut on a cryostat-microtome (Leica CM3050S, Leica Microsystems), collected on Superfrost™ Microscope Slides and stored at −20 °C. Sections were then de-gelatinized for 45 min in PBS at 37°C. Organoid/spheroid sections were incubated in 0.1 M Glycine (Millipore) for 10 min at RT, permeabilized with 0.1%(v/v) Triton X-100 (Sigma) in PBS at RT for 10 min and blocked with 10%(v/v) fetal bovine serum (FBS, Gibco) in TBST (20 mM Tris-HCl pH 8.0, 150 mM NaCl and 0.05%(v/v) Tween-20, Sigma), at RT for 30 min. Organoids/spheroids were then incubated with the primary antibody diluted in 10%(v/v) FBS in TBST solution (Table S1) at 4 °C overnight. Secondary antibodies were added for 30 min in dark, and nuclear counterstaining was performed using 4,6-diamidino-2-phenylindole (DAPI, 1.5 μg/mL, Sigma), at RT for 5 min.

Whole mount

EMOs were incubated in 0.1 M Glycine for 30 min at RT, permeabilized with 0.1%(v/v) Triton X-100 (Sigma) in PBS at RT for 30 min and blocked with 10%(v/v) FBS in TBST, at RT for 1 hour. Organoids were then incubated with the primary antibody diluted in 10%(v/v) FBS in TBST solution (Table S1) at 4 °C overnight, followed by secondary antibody incubation at 4 °C overnight, in the dark.

Imaging quantitative analysis

For quantification of the CD31-positive area in EMOs, organoid sections stained for CD31 and DAPI were analyzed. Using ImageJ software, the DAPI signal was used to define the total area of each organoid section, while the area delineated by CD31 staining was measured to calculate the percentage of the total area occupied by CD31.

Light sheet microscopy of vascularized and innervated EMOs

Fluorescence images of vascularized and innervated EMOs were acquired using a Light Sheet Fluorescence Microscope (LSFM), the Zeiss Lightsheet Z.1, with solid states laser lines 488 nm and 561 nm selected for excitation of Alexa-488 and tdTomato or Alexa-546 and associated filter sets 505-545 nm and 575-615 nm, respectively, equipped with a W Plan-Apochromat 20x/1.0 Corr DIC M27 75 mm objective and a pco.edge 5.5 sCMOS camera. The organoids were embedded in 1 %(m/v) low melting-point agarose, and a glass capillary was used as a sample embedding container in the water chamber. Z-stacks from 4 viewing angles were acquired and fused with ZEN (Zeiss) software to generate a 3D volume rendering image of the entire sample.

EMOs vascular network functional assessment

To assess the permeability of the vasculature in EMOs, after 10 days of co-culture, EMOs were incubated overnight with fluorescein-conjugated dextran (40 kDa; 500 μg/ml; Sigma, FD40) dissolved in culture medium. The following day, images were acquired using a Leica DMI 3000B microscope equipped with a Nikon DXM 1200F digital camera. Subsequently, the samples were fixed for immunostaining analysis. To evaluate vascular remodelling in response to hypoxic stimulus, EMOs after 10 days of co-culture were cultured under hypoxic conditions (5% O_2_) for 48 hours. After this period, EMOs were fixed, sectioned, and subjected to immunostaining analysis.

Multielectrode array (MEA) recording of SNSs

SNSs were plated in a Matrigel-coated high-density MEA system comprising 26400 electrodes per well (MaxWell Biosystems AG; MX2-S-6W). SNSs were cultured for 12-20 days before data recording. Network and axon tracking analysis was performed to attain information regarding axon length, neuron conduction velocity and neuron firing rate, using a MaxWell Biosystems software. During recordings, the temperature was kept at 37°C. To assess SNSs response to nicotine, the media of SNSs cultured in MEA chips was supplemented with 10 μM of nicotine and SNs activity was analyzed after 30 min of incubation.

Norepinephrine and dopamine ELISA

SNSs media after three days of culture was collected and centrifuged at 300 g for 5 mins to eliminate cells or debris. Samples were stored at -80°C until analyzed. Total norepinephrine and dopamine levels in the samples were quantified with a norepinephrine and dopamine ELISA kits (Abnova) according to manufacturer’s instruction.

scRNA-seq analysis

Cell preparation for scRNA-seq

V-EMOs and P-EMOs were dissociated to single cells using 0,25% Trypsin-EDTA and incubation for 10 min at 37°C. After centrifugation for 3 min at 200g, suspended cells were stained with LIVE/DEAD™ Fixable Near-IR Dead Cell Stain (ThermoFischer Scientific, USA) for 15 min. Cells were washed with PBS, centrifuged at 200g for 3 min and resuspended in 300 μl of 0.5% BSA in PBS. Stained cells were sorted using an Aurora CS spectral cell sorter (Cytek Biosciences, USA) at the Flow Cytometry Platform of Gulbenkian Institute for Molecular Medicine. For sorting, a 100 μm nozzle was used with a Purity sort mode (100% Purity Mask and 25% Yield Mask). The sort chamber was previously cooled to 4°C. After exclusion of debris and doublets, live cells (negative for the Dead Cell Stain) were sorted into 1.5 mL tubes previously filled with 30 μL of PBS with 2% BSA. Sorted populations were carefully homogenized with a micropipette and centrifuged at 4°C and 200g for 3 min. Cell pellet was resuspended in PBS and cells were quantified using the Trypan Blue Exclusion Method.

scRNA-sequecing

To generate single-cell gene expression libraries, the Chromium Next GEM Single Cell 3' Reagent Kits v3.1 (10x Genomics) was used, according to the manufacturer’s protocol. Final libraries were quantified using Qubit (Thermo Fisher Scientific) and assessed for size distribution on an Agilent TapeStation (Agilent Technologies). Libraries were pooled and sequenced on a NextSeq 2000 using a paired-end run with the following read configuration: 28 bp (Read 1, cell barcode and UMI), 8 bp (i7 index), and 91 bp (Read 2, transcript). Libraries were sequenced to a depth of 30,000–50,000 reads per cell.

scRNA-seq data pre-processing

The Cell Ranger pipeline (v9.0.1) was used for sample demultiplexing, barcode processing, and generation of the single-cell gene–barcode UMI count matrix. Briefly, samples were demultiplexed to produce paired FASTQ files for each library. Reads were aligned to the GRCh38 human reference transcriptome (Ensembl gene annotation) provided with Cell Ranger. PCR duplicates were removed by matching the same unique molecular identifier (UMI), 10x barcode, and gene, and collapsing them into a single UMI count.

scRNA-seq data post-processing

Downstream analyses, including quality control, normalization, clustering, and differential expression, were performed using scStudio, a web-based application for comprehensive and modular scRNA-seq data analysis (https://www.biorxiv.org/content/10.1101/2025.04.17.649161v1), in R version 4.1.2. A total of 22,400 cells (11,561 from V-EMOs and 10,939 from P-EMOs) were imported, merged, and processed collectively. Genes detected in fewer than five cells and cells expressing fewer than 200 genes were removed. Cells with >40% of transcripts mapping to mitochondrial genes were excluded to eliminate putative dead or lysed cells. This permissive cutoff was chosen due to the naturally high mitochondrial content of cardiomyocytes. Doublets were identified and removed using the scDblFinder package (version 1.8.0, https://bioconductor.org/packages/release/bioc/html/scDblFinder.html). In total, 19,024 cells remained for downstream analysis. Library size normalization was performed using the scran package (version 1.22.1, https://bioconductor.org/packages/release/bioc/html/scran.html) to account for variations in sequencing depth per cell.

Dimensionality reduction and clustering analysis

We performed clustering and dimensionality reduction analyses using the Seurat package (version 4.1.1, https://satijalab.org/seurat/). Highly variable genes (n = 1,000) were selected with the FindVariableFeatures function using the “vst” method. Principal component analysis (PCA) was performed on these genes, and the top 25 principal components (PCs) were used to generate t-SNE and UMAP embeddings. Clustering was performed using shared nearest neighbor (SNN) modularity optimization based on the same 25 PCs. To explore cluster stability, we tested resolutions from 0.1 to 2.0 (step size = 0.1) and visualized the results with the clustree package (version 0.5.1, https://cran.r-project.org/web/packages/clustree/index.html), which highlights relationships between clusters across resolutions. Final clusters were selected based on both clustree output and the expression patterns of canonical markers recently identified markers from available sequencing datasets of human fetal epicardium and human heart organoids.20^,^22^,^23

Differential expression analysis

Differential gene expression analysis was performed on the normalized dataset using the FindAllMarkers or FindMarkers functions in Seurat, applying the MAST statistical framework. Analyses were conducted at the cluster level for each subset to identify genes significantly enriched or depleted in specific cell populations.

Trajectory inference analysis

Cellular differentiation trajectories were reconstructed using the TSCAN package (version 1.44.0, https://bioconductor.org/packages/3.14/bioc/html/TSCAN.html) in R version 4.3.3. Following cells clustering based on PCA components, an unsupervised minimum spanning tree (MST) was generated to infer pseudotime ordering. Trajectories were visualized with cluster annotations, and pseudotime-associated genes were identified using TSCAN’s built-in analysis of variance (ANOVA) function.

RT-qPCR analysis of SN progenitors

Total RNA was isolated from SN progenitors both at days 6 and 12 of differentiation, using a High Pure RNA Isolation Kit (Roche) according to the manufacturer's instructions. Total RNA was converted into cDNA with High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific). PCR reactions were performed with Syber Green Master Mix (nzytech). The primer sequence is available in the Supplementary Data (Table S2). Reactions were run in triplicate in ViiA7 Real-Time PCR Systems (Applied BioSystems). For each analyzed time point, gene expression was normalized against the expression of the housekeeping gene glyceraldehyde-3-phosohate dehydrogenase (GAPDH) and results were analysed with the QuantStudioTM RT-PCR Software.

Calcium transient imaging of EMOs

Fluo-4-AM dye (ThermoFisher) was used for calcium imaging of EMOs. For that, Fluo-4-AM dye-loading solution was prepared according to the manufacturer’s instructions and EMOs were incubated with that solution for 30 min at 37°C. Afterwards, the Fluo-4-AM dye-loading solution was exchanged by DMEM/F12 medium. Before starting image acquisition, EMOs were incubated at 37°C for another 30 min to stabilize. Videos of beating EMOs were taken for a period of 30–60 s with a ZEISS Cell Observer SD confocal microscope. Calcium transient profiles were analyzed using an in-house Python-developed software.

Quantification and statistical analysis

All presented data regarding IF staining, MEA and PCR analysis were representative of at least three biologically independent experiments per condition/marker. All raw data was collected in Microsoft Excel and statistical analysis was performed with GraphPad software. Statistical significance was evaluated with a two-tailed, unpaired Student’s t test (p < 0.05) when appropriate. The data are presented as mean ± SEM and represent a minimum of three independent experiments. Statistical significance was assigned as not significant (ns), p > 0.05; ∗p ≤ 0.05; ∗∗p ≤ 0.005; ∗∗∗p ≤ 0.0005. All micrograph images are representative of at least four independent experiments per condition/marker.