An integrative single-nucleus multiomic atlas of the human left ventricle identifies gene regulatory network dynamics across cardiac development, aging, and disease.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Genome Biol . 2026 Apr 6;27:126. doi: 10.1186/s13059-026-04061-7 Search in PMC Search in PubMed View in NLM Catalog Add to search An integrative single-nucleus multiomic atlas of the human left ventricle identifies gene regulatory network dynamics across cardiac development, aging, and disease William Gao William Gao 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 7 Medical Scientist Training Program, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by William Gao 1, 2, 3, 4, 7, ✉, # , Peng Hu Peng Hu 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 9 Present Address: Key Laboratory of Exploration and Utilization of Aquatic Genetic Resources, Ministry of Education, Shanghai Ocean University, Shanghai, 201306 China Find articles by Peng Hu 1, 9, # , Brittney Wick Brittney Wick 8 Genomics Institute, University of California, Santa Cruz, Santa Cruz, CA 95060 USA Find articles by Brittney Wick 8 , Qi Qiu Qi Qiu 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Qi Qiu 1, 2, 3, 4 , Hongjie Zhang Hongjie Zhang 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Hongjie Zhang 1, 2, 3, 4 , Ying Li Ying Li 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Ying Li 1, 2, 3, 4 , Xiangjin Kang Xiangjin Kang 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Xiangjin Kang 1 , Kenneth Bedi Kenneth Bedi 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 5 Department of Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Kenneth Bedi 4, 5 , Maximilian Haeussler Maximilian Haeussler 8 Genomics Institute, University of California, Santa Cruz, Santa Cruz, CA 95060 USA Find articles by Maximilian Haeussler 8 , Kotaro Sasaki Kotaro Sasaki 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 6 Department of Pathology and Laboratory Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 10 Department of Biomedical Sciences, University of Pennsylvania, School of Veterinary Medicine, Philadelphia, PA 19104 USA Find articles by Kotaro Sasaki 3, 6, 10 , Kenneth Margulies Kenneth Margulies 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 5 Department of Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Kenneth Margulies 4, 5 , Hao Wu Hao Wu 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA Find articles by Hao Wu 1, 2, 3, 4, ✉ Author information Article notes Copyright and License information 1 Department of Genetics, University of Pennsylvania, Philadelphia, PA 19104 USA 2 Penn Epigenetics Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 3 Penn Institute of Regenerative Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 4 Penn Cardiovascular Institute, University of Pennsylvania, Philadelphia, PA 19104 USA 5 Department of Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 6 Department of Pathology and Laboratory Medicine, University of Pennsylvania, Philadelphia, PA 19104 USA 7 Medical Scientist Training Program, University of Pennsylvania, Philadelphia, PA 19104 USA 8 Genomics Institute, University of California, Santa Cruz, Santa Cruz, CA 95060 USA 9 Present Address: Key Laboratory of Exploration and Utilization of Aquatic Genetic Resources, Ministry of Education, Shanghai Ocean University, Shanghai, 201306 China 10 Department of Biomedical Sciences, University of Pennsylvania, School of Veterinary Medicine, Philadelphia, PA 19104 USA ✉ Corresponding author. # Contributed equally. Received 2025 Oct 20; Accepted 2026 Mar 24; Collection date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/ . PMC Copyright notice PMCID: PMC13067603  PMID: 41937210 Abstract Background As the first organ to develop in utero, the human heart undergoes extensive molecular, structural and metabolic remodeling during development and must sustain its function throughout life. Results We generate an integrated multiomic atlas of human cardiac cells, combining newly generated and publicly available single-nucleus RNA sequencing datasets from 299 donors and single nucleus ATAC-seq datasets from 106 donors. Developmental and disease-associated processes drive far more extensive molecular remodeling than sex-associated or aging-dependent effects. Across nearly all cardiac cell types, developmental and disease-driven changes exhibit strong overlap at both the transcriptomic and epigenomic levels, revealing widespread reactivation of fetal-associated gene programs beyond cardiomyocytes. Both cardiac development and disease show convergent shifts in intercellular communication, including increased TGFβ signaling. Integration of gene expression and chromatin accessibility data reveals putative cell-type–specific transcriptional factors driving fetal reactivation in major cardiac diseases. Spatial transcriptomics data orthogonally identifies localization of this fetal reactivation signature within spatially distinct niches in ischemic and fibrotic zones of acute myocardial infarction. Finally, we construct a cell-type–resolved enhancer-to-gene linkage map that refines the association of dilated and hypertrophic cardiomyopathy genetic risk loci to downstream target genes. Conclusions This study presents a comprehensive multimodal, cell-type–resolved atlas of the human heart, providing a foundation for understanding human cardiac gene regulation across the human lifespan and in cardiac diseases. Supplementary Information The online version contains supplementary material available at 10.1186/s13059-026-04061-7. Keywords: Cardiovascular disease, Single cell biology, Multiomics, Epigenomics, Transcriptomics, Single-cell atlas integration Background Cardiovascular disease (CVD) is the leading cause of death worldwide. Age is among its greatest risk factors, with CVD mortality rising exponentially after age 65 and with earlier and higher incidence observed in males [ 1 ]. Given its global burden, numerous studies have investigated the mammalian heart across development, aging, and disease states, including analyses of human donor hearts. Early transcriptomic studies using microarrays and bulk RNA-sequencing [ 2 – 4 ] provided valuable insights but lacked resolution. The advent of single cell RNA-sequencing (scRNA-seq) has since enabled cell-type–specific analysis of cardiac tissues [ 5 ]. In particular, the development and application of high-throughput single-nucleus RNA-seq (snRNA-seq) has made it possible to profile all major cardiac cell types in postnatal mammalian hearts in an unbiased manner, including large-sized cardiomyocytes [ 6 ]. Since the publication of the first large-scale human heart snRNA-seq studies in 2020 [ 7 , 8 ], there has been rapid expansion of single-cell/nucleus transcriptomic studies using human donor hearts. In parallel, single-cell/nucleus epigenomic technologies such as the assay for transposase-accessible chromatin (scATAC-seq) have enabled genome-wide mapping of chromatin accessibility at single-cell/nucleus resolution [ 9 ]. These methods can reveal transcription factor (TF) binding activity and gene regulatory dynamics and link cis -regulatory elements to downstream genes within complex tissues such as human hearts [ 10 , 11 ]. Despite these advances, no study has comprehensively compared transcriptional and epigenomic landscapes of the human heart across its full phenotypic spectrum, encompassing early development, sex differences, aging, and disease states, largely because individual studies remain limited by donor numbers and thus lack sufficient statistical power. To address this gap, we performed an integrated single-cell multiomic analysis of the human heart, combining newly generated datasets with publicly available studies. Performing such large-scale integration presents unique challenges, including pronounced inter-study batch effects [ 12 ] and computational bottlenecks, as many analytical tools do not scale efficiently to atlas-scale datasets typically comprising millions of cells [ 13 ]. Nonetheless, integrated analysis provides several key advantages. Most prior studies included fewer than 50 donors, limiting statistical power for robust differential gene analyses [ 14 ]. By aggregating multiple datasets, we increase sensitivity to detect subtle biological signals and distinguish reproducible, “ground truth” molecular signals from study-specific artifacts. Furthermore, recent advances in computational methods now allow more efficient, scalable and accurate integration and interpretation of large single-cell datasets. Here, by jointly analyzing 299 snRNA-seq, 106 snATAC-seq, and 27 spatial transcriptomic samples from human hearts, we identify robust molecular changes across cardiac development, aging, and disease. These analyses systematically identify key genes, pathways, and cell-type–specific regulatory networks that may underlie developmental remodeling, age-dependent changes, and disease states in the human heart. Results Integration of ~ 2.3 million snRNA-seq profiles of human cardiac nuclei across 299 donors To investigate gene expression dynamics underlying human cardiac maturation, aging, and disease states in a cell-type–specific manner, we constructed an integrated snRNA-seq atlas combining newly generated and publicly available datasets (Additional file 1 : Fig. S1a). Our new snRNA-seq dataset has 19 donors (Additional file 1 : Fig. S2), including two fetal, 15 non-diseased (ND) postnatal, and two dilated cardiomyopathy (DCM) donors, profiled using the droplet microfluidics-based sNucDrop-seq method [ 6 ]. Because cardiomyocytes are too large for standard droplet-based encapsulation protocols [ 5 , 6 ], most studies used snRNA-seq rather than scRNA-seq and focused on the left ventricle (LV) of the human heart (Fig. 1 a and Additional file 1 : Fig. S1a). Accordingly, we focused our analysis on LV snRNA-seq datasets, except for fetal samples in which the whole heart was typically sampled. Fig. 1. Open in a new tab Integrated human cardiac single-nucleus RNA-seq atlas. a Transcriptomic atlas includes single-nucleus RNA-sequencing (snRNA-seq) from 299 human donors, with phenotypic diversity across age, sex and disease states. b Uniform manifold approximation and projection (UMAP) plot of 2,275,105 snRNA-seq nuclei across 13 cell types. c Marker gene expression for each annotated snRNA-seq major cardiac cell type. d Subclustering of cardiomyocytes based on snRNA-seq. e Subclustering of fibroblasts based on snRNA-seq. Abbreviations: LEC: lymphatic epithelial cells; vSMC: vascular smooth muscle cells; DCM: dilated cardiomyopathy; HCM: hypertrophic cardiomyopathy To minimize batch effects, we prioritized datasets that met one of these criteria: (1) ≥ 10 donors, or (2) ≥ 3 donors across at least two disease/developmental categories (fetal ND, postnatal ND, and postnatal diseased). Ten studies encompassing 280 donors satisfied these criteria [ 8 , 11 , 15 – 24 ] (Additional file 2 : Table S1), which together with our new data, yielded a total of 299 snRNA-seq LV datasets. Collectively, these datasets are diverse by sex (111 female, 188 male), age group (13 fetal, 57 young (≤ 40 yr), 128 middle-aged (40–59 yr), 101 old (≥ 60 yr)), and disease status (167 ND, 132 diseased) (Fig. 1 a; Additional file 3 : Table S2). We selected age cutoffs of “young” as age ≤ 40 and “old” as age ≥ 60 to potentially capture non-linear changes that may occur between 40 and 60 [ 25 ] and since CVD risk rises markedly after age 60 (1). These age group cutoffs were similar to cutoffs determined using two other approaches (k-means with k = 3, or tertiles) (Additional file 1 : Fig. S1b), confirming robustness to alternative definitions of aging cohorts. Diseased donors included many types of heart disease: acute myocardial infarction (AMI), arrhythmogenic cardiomyopathy (ARVC), dilated cardiomyopathy (DCM), hypertrophic cardiomyopathy (HCM), ischemic cardiomyopathy (ICM), and non-compaction cardiomyopathy (NCCM), with the most representation for DCM, HCM, and ICM (Additional file 1 : Fig. S1c). After combining processed count matrices provided by each individual study and performing standard preprocessing for our newly generated dataset ( Methods ), we retained 2,275,105 nuclei after quality control filtering (referred from hereafter as the “aggregated dataset”; Fig. 1 b; Additional file 1 : Fig. S3a). Using the scIB integration benchmarking framework [ 13 , 26 ], we confirmed that scVI [ 27 ] achieved the best performance across metrics measuring batch effect correction while preserving true biological heterogeneity (Additional file 1 : Fig. S3b-c). Therefore, we used scVI, which scales efficiently to millions of cells/nuclei, for data integration and batch correction in both the Uniform Approximation and Manifold Projection (UMAP) embedding and the count matrix [ 28 ]. At the embedding level, scVI integration demonstrated batch removal of technical artifacts such as study and technology (Additional file 1 : Fig. S3d-e), while retaining biological heterogeneity, such as within cell type differences across developmental age and disease status (Additional file 1 : Fig. S3f). Additionally, though we focused on nuclear transcriptomes for this atlas given the much larger number of snRNA-seq datasets compared to scRNA-seq (limited by the depletion of cardiomyocytes due to their large size), scVI partially integrated cells and nuclei of the same cell type together. In a subsampled dataset containing 50,000 single cells and 50,000 single nuclei, cells and nuclei for the same cell type did not fully mix at the embedding level, but this may reflect the expected biological preservation of differences between nuclear and whole cell transcriptomes (Additional file 1 : Fig. S3g). Using this carefully integrated and benchmarked human cardiac snRNA-seq atlas (Additional file 1 : Fig. S3), we re-annotated cell types independently of the original studies, using literature-curated canonical marker genes. We identified 13 major cardiac cell types: adipocyte, cardiomyocyte, endocardial, endothelial, epicardial, fibroblast, lymphatic epithelial (LEC), lymphoid, mast, myeloid, neuronal, pericyte, and vascular smooth muscle (vSMC) cells (Fig. 1 c). Atrial cardiomyocytes, which were identified in some fetal whole heart samples, separated from ventricular cardiomyocytes in the embedding and were not included in further analysis. Consistent with previous studies, we find that the most prevalent snRNA-seq cell types in the human heart are cardiomyocytes, fibroblasts, endothelial cells, pericytes, and myeloid cells. However, since some cardiomyocytes contain more than one nucleus [ 29 ], the number of cardiomyocyte nuclei overestimates the number of cardiomyocyte cells. The rarest cell types are LECs, adipocytes, mast cells, and epicardial cells (Additional file 1 : Fig. S3h). Our revised annotations showed > 95% concordance between our revised annotation and the original annotations for most cell types (Additional file 1 : Fig. S3i). Discrepancies primarily reflected differences between transcriptionally similar cell types (e.g., endothelial vs. endocardial vs LEC; vSMC vs. pericyte). In rarer cell types such as Mast cells where the original and new annotation have higher discordance (Additional file 1 : Fig. S3j), the nuclei uniquely annotated in our study were more transcriptionally correlated with consensus-labeled cells than those annotated only in the original datasets (Additional file 1 : Fig. S3k), suggesting that our revised annotations for rarer cell types may be more accurate through large-scale integration across independent datasets. We next performed subclustering within each major cell type to identify biologically distinct subpopulations. To ensure biological relevance rather than study-specific technical artifacts, we focused on subclusters that were consistently represented across multiple studies and supported by literature-reported marker genes in the relevant cell type ( Methods ; Additional file 1 : Fig. S4). Using this approach, we identified seven subclusters of ventricular cardiomyocytes (vCMs), including fetal, HCM, and DCM-enriched states. Among these subclusters was a rare (0.05% of cells) population of “senescent” vCMs expressing canonical senescence markers such as CDKN1A and CDKN2A which is enriched in DCM and HCM (Fig. 1 d; Additional file 1 : Fig. S4b). In fibroblasts, we also identified seven subpopulations, including an “activated” subcluster enriched in AMI, DCM, and HCM ( POSTN/TNC/FAP ) and an “inflammatory” subcluster ( IL6/THBS1/NR4A1 ) enriched in pediatric heart failure (Fig. 1 e; Additional file 1 : Fig. S4e). Endothelial cells included six subclusters largely based on vessel identity (arterial, capillary, venous), but included an “activated” subpopulation expressing inducible nitric oxide synthase ( NOS2 ) [ 30 ] and inflammatory markers such as TNFRSF4 , that is enriched in fetal, DCM, HCM, and ICM hearts (Additional file 1 : Fig. S4d). In lymphoid cells, we identified eight subpopulations spanning natural killer (NK) cells, T cells, and B cells (Additional file 1 : Fig. S4g). Notably, we observed an age-associated decline in non-plasma B and plasma B cells, consistent with the features of immunosenescence [ 31 ]. In myeloid cells, we identified six subclusters, including anti-inflammatory ARG1 + macrophages, CCR2 + monocyte-derived macrophages, cycling ( CDK1/EZH2/TOP2A ) resident macrophages ( FOLR2 + /LYVE1 + ) (Additional file 1 : Fig. S4h ) . Monocyte-derived macrophages are enriched in AMI [ 32 ], cycling cells are most enriched in fetal hearts, while the ARG1 + anti-inflammatory macrophages were moderately enriched in DCM. Many of these subclusters comprise < 5% of each cell type, demonstrating how large-scale integration across studies enables detection of rare yet reproducible cell states in human hearts. Integration of ~ 700,000 snATAC-seq profiles of human cardiac nuclei across 106 donors Gene expression is orchestrated by cell-type–specific transcription factors (TFs) that bind to cis -regulatory elements. To better understand epigenetic regulation of cell-type–specific gene expression in the human heart, we constructed an integrated epigenomics atlas consisting of LV snATAC-seq datasets across developmental and disease contexts, resulting in snATAC-seq datasets from 95 human donors [ 10 , 11 , 16 ] and newly generated snATAC-seq data from 11 additional donors (5 fetal and 6 postnatal ND donors) (Additional file 1 : Fig. S5). Altogether, the aggregation of 106 total donors (Fig. 2 a; Additional file 1 : Fig. S6a; Additional file 4 : Table S3) enables a large, integrated analysis of human cardiac snATAC-seq datasets to date. This integrated dataset is also diverse in terms of sex (38 female, 68 male), age groups (17 fetal, 14 young, 44 middle-aged, 31 old) and disease status (93 ND, 13 diseased). All diseased snATAC-seq donors were from a single study that profiled ICM and AMI [ 11 ]. Fig. 2. Open in a new tab Integrated human cardiac single-nucleus ATAC-seq atlas. a Epigenomic atlas includes single-nucleus ATAC-sequencing (snATAC-seq) from 106 human donors, also with phenotypic diversity. b UMAP plot of 690,044 snATAC-seq nuclei across 11 cell types. c Pseudobulked chromatin accessibility of marker genes for each annotated snATAC-seq cell type. d Genomic feature annotation of snATAC-seq peaks using ChIPSeeker. e Observed overlap between non-promoter distal peaks (> 1 kb away from TSS) and Spurrell et al. 2022 enhancers, compared to a null distribution of expected overlap of shuffled peaks of same sizes. Abbreviations: RPKM: reads per kilobase million, UTR: untranslated region; TSS: transcription start site We preprocessed raw fragment files from each individual study using SnapATAC2 [ 33 ] and applied Harmony [ 34 ] for batch correction, which has recently been shown to outperform other ATAC integration approaches [ 35 ] (Additional file 1 : Fig. S6b). We then performed Leiden clustering and cell type annotation using two complementary approaches: (1) MAGIC, which provides an imputed gene expression based on chromatin accessibility, and (2) label transferring from matched RNA-based annotations in multimodal datasets. While these approaches were largely concordant, the label transfer method better recovered rarer cell types such as adipocytes and neuronal cells (Additional file 1 : Fig. S6c). After stringent quality control filtering and batch integration, we retained 690,044 high quality nuclei across 11 major cardiac types (Fig. 2 b; Additional file 1 : Fig. S6d-e, Methods ). We then identified cell-type–specific peaks of accessible chromatin using MACS3 [ 36 ]. Promoters of known marker genes displayed markedly higher accessibility within their corresponding cell types (Fig. 2 c), confirming that snATAC-seq robustly captures cell-type–specific cis -regulatory elements. About a quarter of peaks were distal, located > 1 kilobase from transcriptional start sites (TSSs) (Fig. 2 d). These distal intergenic peaks displayed 97.6% overlap with enhancers identified in bulk H3K27ac ChIP-seq analysis of human heart tissues (2) (Fig. 2 e), demonstrating snATAC-seq’s ability to capture active distal cis -regulatory elements. Together, our integration of ~ 2.3 million snRNA-seq nuclei and ~ 700 K snATAC-seq nuclei establishes the most comprehensive single-nucleus multiomic atlas of the human heart to date. This resource provides an unprecedented opportunity to dissect transcriptional and epigenetic regulation across development, aging, and disease. The integrated atlas is available as an online portal hosted by the UCSC Cell Browser [ 37 ] ( https://multiomic-human-heart.cells.ucsc.edu ; Methods ). Systematic benchmarking of statistical approaches for robust differential gene expression and chromatin accessibility analysis Leveraging the phenotypic diversity in our multiomic atlas, we next sought to identify cell type-specific differentially expressed genes (DEGs) and differentially accessible regions (DARs) across biological states, starting with sex-associated differences in non-diseased donors. To perform this analysis at the cell-type level, we pseudobulked all nuclei from the same donor and cell type (Fig. 3 a). This computationally efficient strategy substantially decreased the number of false positive DEGs compared to treating single cells/nuclei as individual observation units, which can lead to inflated significance due to within-donor correlations [ 12 , 14 , 38 ]. Additionally, prior work has shown that performing DEG and DAR analysis on raw counts while treating batch effects as covariates, yields more accurate results than using batch-corrected “denoised” counts [ 12 , 39 ]. Fig. 3. Open in a new tab Cell-type–resolved molecular changes across sex and aging. a Pseudobulked profiles for each cell type were obtained by combining raw counts for all cells belonging to the same cell type for each donor. This pseudobulked count matrix was used for differential gene expression analysis with DESeq2. b Number of sex DEGs by cell type. c Number of sex DARs by cell type. d FRMD5 cardiomyocyte pseudobulked promoter accessibility is higher in males. e FRMD5 cardiomyocyte pseudobulked gene expression is higher in males. f FRMD5 bulk gene expression in left ventricle (from GTEx left ventricle) is not statistically significant between sexes. g FRMD5 is highly expressed in cardiomyocytes and neurons, with sex-associated differences only present in cardiomyocytes. h Number of age-associated DEGs by cell type. i Number of age-associated DARs by cell type. j Promoter accessibility of sterol response element binding factor 1 ( SREBF1 ) significantly decreases across age groups. k Gene expression of SREBF1 significantly decreases across age groups. l Low-density lipoprotein receptor ( LDLR ), a target gene of SREBF1 , also significantly decreases across age groups. Abbreviations: ND: non-diseased; DEGs = differentially expressed genes; DARs = differentially accessible regions; * ( p < 0.05), ** ( p < 0.01), *** ( p -value < 0.001); n.s. = non-significant As an alternative statistical approach, we evaluated the mixed effects model NEBULA on the single-nucleus raw count matrix, treating donor, sex, age, technology, and study as fixed effects [ 40 ]. However, we found that NEBULA was sensitive to ambient RNA contamination, often detecting spurious upregulation of cardiomyocyte marker genes ( TTN , RYR2 ) with age in many other cell types. Conversely, NEBULA showed reduced sensitivity for detecting sex chromosome genes as sex-associated in rarer cell types (Additional file 1 : Fig. S7). Since the mixed effect model leverages single cell-level variability, NEBULA may be more prone to false negatives when the expression is low (as for sex chromosome genes), and to false positives in the presence of ambient RNA, an artifact common in nuclear preparations and incompletely corrected by current denoising methods like SoupX, which was used in our analysis [ 41 , 42 ]. Based on systematic benchmarking between these approaches, we therefore identified DESeq2 [ 43 , 44 ] for pseudobulk-based analysis as the more robust and interpretrable approach for downstream DEG and DAR analysis. Notably, while the DESeq2 model may perform better here due to its lower complexity, which leads to more stable parameter estimation, the increased complexity of the NEBULA may be advantageous in other contexts. Since the number of pseudobulked transcripts varies with cell number and sequencing depth, we performed downsampling of pseudobulked cell-type–specific profiles using donors with the highest number of transcripts. This established that 50,000 transcripts were sufficient for obtaining a Spearman correlation > 0.8 to the full pseudobulked expression regardless of cell type (Additional file 1 : Fig. S8a-b). Similarly, 300,000 fragments were sufficient for pseudobulked analysis of snATAC-seq datasets (Additional file 1 : Fig. S8d-e). We therefore restricted DEG and DAR analyses to donors satisfying these coverage thresholds (≥ 50,000 pseudobulked transcripts and 300,000 fragments for each cell type). To ensure statistically well-powered analysis, we included only cell types with at least 30 donors, thereby performing DEG analysis for all 13 cell types and DAR analysis for five major cardiac types (cardiomyocyte, endothelial, fibroblast, myeloid, and pericyte) (Additional file 1 : Fig. S8c, f). For sex- and age-specific analyses, we modeled expression as a function of sex (male or female) and age group (young, middle, old) using data from non-diseased donors across all studies. We regressed out technology and study (“tech_plus_study”) as covariates, which were identified as the major batch effect (Additional file 1 : Fig. S9a-b). We refer to this as the “batch as covariate” approach. Because batch regression may inadvertently remove part of the biological signal, we also tested a “meta-analysis” strategy: running DESeq2 for sex-associated analysis individually for studies with at least three male and female donors using a weighted Fisher’s meta-analysis for aggregating p -values. To benchmark performance against a silver standard, we compared DEG results from fully pseudobulked “batch as covariate” and “meta-analysis” against the LV bulk RNA-seq from GTEx [ 4 ] (Additional file 1 : Fig. S9c). Both methods demonstrated high specificity (> 0.99) (Additional file 1 : Fig. S9d-g), but the “batch as covariate” method identified more “true” DEGs overall, with higher sensitivity, consistent with a recent study [ 45 ]. The “meta-analysis” approach yielded fewer DEGs, largely due to the limited number of studies with balanced donor sex (Additional file 1 : Fig. S9l). For both approaches, overall sensitivity was modest, reflecting technical differences in assay type as well as biological differences between whole cells in bulk tissue and single nuclei [ 46 ]. For sex chromosome genes (“gold standard”), both approaches achieved 100% sensitivity. Comparable trends were observed in aging analyses (“old” vs “young”) (Additional file 1 : Fig. S9h-k, m), where the reduced number of age-balanced studies further limited meta-analysis sensitivity. Given the similar specificity and superior sensitivity for “batch as covariate” approach, we selected this method for DEG and DAR analysis in the rest of this study. Atlas-wide analysis reveals sex- and aging-associated differences in gene expression and chromatin states Using the “batch as covariate” approach, we first examined sex-associated transcriptomic and epigenomic differences between male and female donors across major cardiac cell types. Overall, we detected fewer than 30 sex-associated DEGs in most cell types (Fig. 3 b; Additional file 5 : Table S4). The number of identified sex- and age-associated DEGs was lower than in a recent study [ 47 ]. This discrepancy likely reflects our larger donor cohort ( n = 154 for this study vs. n = 73 for the previous study), which favors reproducible signals across studies, as well as methodological differences in the previous study, which used NEBULA (Additional file 1 : Fig. S10). As expected, DEGs were predominantly sex-chromosome genes, including X-linked (e.g. XIST, TSIX, KDM6A ) and Y-linked (e.g., DDX3Y, KDM5D, LINC00278 ) genes (Additional file 1 : Fig. S11a-d). Similarly, DAR analysis also identified few sex-associated chromatin changes outside of the expected sex chromosomes (Fig. 3 c; Additional file 1 : Fig. S11e-g; Additional file 5 : Table S4). Within cardiomyocytes, we consistently identified FRMD5 , an autosomal gene associated with ataxia [ 48 ] and lipid metabolism [ 49 ], as a reproducible, male-upregulated DEG in every individual study with at least three male and female donors (Additional file 1 : Fig. S11h). Both promoter accessibility and gene expression of FRMD5 were significantly higher in males (Fig. 3 d-e). Interestingly, FRMD5 showed no significant sex-associated difference at the bulk tissue level in GTEx data (Fig. 3 f), likely because of higher but non-differential sex expression in neurons (Fig. 3 g). This underscores the power of our integrated snRNA-seq atlas to uncover cell-type–specific expression regulatory changes masked in bulk tissue analyses. Finally, while the number of sex-differentially expressed genes was low, we also performed gene set enrichment analysis to identify pathway-level differences between males and females. This identified that immune pathways and oxidative phosphorylation were upregulated in males for several cell types perhaps due to different metabolic demands across sexes (Additional file 1 : Fig. S11i), which have known differences in LV mass and ejection fraction [ 50 ]. We next examined cell-type–specific transcriptomic and epigenomic changes associated with aging by comparing old versus young non-diseased donors. DEG analysis revealed that cardiomyocytes, adipocytes, and fibroblasts exhibited the most transcriptional changes, followed by myeloid cells and neuronal cells (Fig. 3 h; Additional file 6 : Table S5). Notably, in contrast to sex DEGs for all cell types besides cardiomyocytes, the number of aging DEGs across cell types increased as a function of number of donors and did not saturate even when all 299 donors were included, suggesting that additional donors would likely reveal further age-associated signals (Additional file 1 : Fig. S12a). This trend also helps explain why rarer cell types with fewer donors tended to have fewer significant DEGs. DAR analysis in cell types with enough snATAC-seq donors corroborated that cardiomyocytes and fibroblasts showed the most aging-associated chromatin changes (Fig. 3 i; Additional file 6 : Table S5). Most age-associated DEGs were cell-type–specific (Additional file 1 : Fig. S12b). Specifically in cardiomyocytes, the sterol response element binding factor 1 ( SREBF1 ), a key regulator of lipid metabolism, was significantly downregulated with age (Fig. 3 j-k; Additional file 1 : Fig. S12c). Several of its known downstream targets (based on SREBF1 ChIP-seq), such as LDLR (low-density lipoprotein receptor) also significantly decreased with age (Fig. 3 l; Additional file 1 : Fig. S12c-d), indicating age-associated reduction in cholesterol uptake that may contribute to metabolic shifts in aging cardiomyocytes. Additionally, we constructed a sex * age interaction model to identify genes that change differentially across age by sex. Given the overall low number of sex and aging genes, there were fewer than 5 such genes across any cell type. However, CSMD1, a gene recently linked to DCM [ 51 ], is a notable gene that increases with age in males while decreasing in females (Additional file 1 : Fig. S12f-g). Though most aging DEGs were cell-type–specific, a subset of age-associated DEGs were shared across multiple cardiac cell types. For example, PTCHD4 , which we identified as a top marker gene for the “senescent” cardiomyocyte subcluster (Fig. 1 d), increased with age across seven human cardiac cell types, including cardiomyocytes and fibroblasts (Additional file 1 : Fig. S12c, e). While PTCHD4 has been implicated in cellular senescence [ 4 , 52 ], we did not observe a significant age-associated increase in canonical senescence markers CDKN1A (p16) and CDKN2A (p21) or in a composite senescence score derived from the SenMayo gene set [ 53 ] in any cell type (Additional file 1 : Fig. S13a), likely due to the low percentage of senescent cardiomyocytes (Fig. 1 d). Given the data sparsity in snRNA-seq, it is also possible that PTCHD4 may be a more reliable senescence marker for single-cell analyses. Aging has also been associated with increased transcriptional noise [ 54 ] and loss of the Y chromosome in males in other tissues [ 55 ]. However, we observed no consistent changes in transcriptional noise, consistent with prior findings that noise estimates are highly metric-dependent [ 56 ] (Additional file 1 : Fig. S13b). We also did not observe age-associated loss of the Y chromosome across cardiac cell types (Additional file 1 : Fig. S13c). At a pathway level, aging was associated with decrease of proliferation and epithelial-mesenchymal transition but surprisingly showed an increase in oxidative phosphorylation across many cell types (Additional file 1 : Fig. S13d). While oxidative phosphorylation is typically associated with decreases across aging, this finding may suggest that non-diseased cardiac aging is not intrinsically associated with pathological mitochondrial dysfunction or could reflect compensatory activation of nuclear-encoded mitochondrial genes downstream of mitochondrial dysfunction, as reported elsewhere [ 57 , 58 ]. Together, these results suggest that healthy human cardiac aging is characterized by relatively limited transcriptional and chromatin remodeling, with the most pronounced effects in long-lived, non-dividing cardiomyocytes. Cardiac cell types exhibit shared regulatory programs during fetal-postnatal transition To better understand human cardiac maturation, we further explored chromatin and gene expression changes during the fetal-postnatal transition. For this analysis, we performed DEG and DAR analysis between fetal and postnatal, non-diseased young donors, identifying thousands of DEGs and DARs across multiple cardiac cell types during the fetal-postnatal transition (Fig. 4 a-b; Additional file 7 : Table S6). The number of shared DEGs across cell types was higher than in sex- or aging-based contrasts, suggesting that many cell types undergo coordinated remodeling of transcriptional landscapes during this critical developmental window (Additional file 1 : Fig. S14a). Several cell types displayed strong fetal-biased expression of insulin-like growth factor binding proteins such as IGF2BP1 and IGF2BP3 , which are critical for cardiac proliferation during development [ 59 ]. In contrast, multiple cell types in the young postnatal heart exhibited increased expression of THRB , thyroid receptor B, which inhibits proliferation and drives perinatal cardiac maturation [ 60 ] (Additional file 1 : Fig. S14b-d). Gene set enrichment analysis further identified that several pathways related to growth and cell division (MYC/E2F targets, G2/M checkpoint) are enriched in the fetal state (Additional file 1 : Fig. S14e), whereas immune and inflammatory pathways such as TNF-alpha and interferon signaling are low during development and increase postnatally, consistent with the emergence of immune competence following a tolerogenic fetal environment [ 61 ]. Trajectory analysis across major cell types corroborated these results, identifying several pathways at different stages along the fetal to postnatal transition in a cell-type–specific manner. For example, the fetal cardiomyocyte stage was enriched for proliferation and cell cycle terms, while the postnatal stage was enriched for fatty acid oxidation and inflammatory responses in cardiomyocytes and endothelial cells (Additional file 1 : Fig. S15; Methods ). Fig. 4. Open in a new tab Molecular changes across human cardiac development and disease. a Number of DEGs across development. b Number of DARs across development. c Transcription factor-driven regulons with increased and decreased activity across development. d Number of disease-associated DEGs, comparing binarized diseased hearts to non-diseased hearts. e Number of disease-associated DARs, comparing binarized diseased hearts to non-diseased hearts. f Gene set enrichment analysis identifies cell-type–specific pathways that are positively or negatively enriched in three disease subtypes: DCM, HCM, and ICM. Abbreviations: DCM: dilated cardiomyopathy; HCM: hypertrophic cardiomyopathy, ICM: ischemic cardiomyopathy We next inferred TF-driven GRN activity using SCENIC +, identifying numerous TFs associated with significant changes during developmental transitions from fetal to postnatal hearts, several of which are shared across multiple cell types. At the fetal stage, we observed significantly higher fetal GRN activity for canonical cardiac developmental TFs such as MEF2A , SOX5/6 [ 62 ], and other related TFs such as TEAD1 , a known regulator of cardiac progenitor proliferation and organ size via the Hippo signaling pathway [ 63 ] (Fig. 4 c). In contrast, postnatal hearts showed increased activity of progesterone steroid receptor ( PGR) and androgen steroid receptor ( AR ). We also observed postnatal up-regulation of FOXO1/FOXO3 , which repress fetal proliferative programs and promote metabolic maturation [ 64 ]. Together, these results reveal that the fetal-postnatal transition involves coordinated, cell-type–specific and shared changes in both gene expression and chromatin states, reflecting a global shift from proliferative to metabolically mature and immunocompetent states in postnatal human hearts. Heart failure subtypes display shared and distinct transcriptional and pathway remodeling To investigate the molecular alterations underlying disease conditions, we next explored disease-associated transcriptional and chromatin state changes. First, we performed DEG analysis between diseased and non-diseased donors in a binarized manner. Like the fetal-postnatal transition, we also identified hundreds to thousands of disease-associated DEGs and DARs across most cardiac cell types (Fig. 4 d-e; Additional file 8 : Table S7). In cardiomyocytes, we observed marked upregulation of known heart failure markers including atrial natriuretic peptide ( NPPA ) [ 65 ] and decrease of MYH6 [ 66 ] (Additional file 1 : Fig. S16a). ACE2 , which reduces blood pressure, inflammation, and adverse cardiac remodeling through the degradation of angiotensin II [ 67 ] and serves as the human receptor for the SAR-CoV2 virus [ 68 ], is specifically downregulated in disease-associated cardiac fibroblasts. Additionally, fibroblasts showed decreased expression of endothelin receptor B, EDNRB , and monoamine oxidase, MAOA , both of which regulate blood pressure (Additional file 1 : Fig. S16b). Surprisingly, the pro-inflammatory alarmins S100A8 and S100A9 are the most significant downregulated myeloid transcripts in diseased hearts (Additional file 1 : Fig. S16c), which may suggest immune dysregulation in end-stage heart failure rather than strong pro-inflammatory state observed in early stages of disease [ 69 ]. Though our snATAC-seq dataset only included AMI/ICM donors, the inclusion of many heart disease subtypes in our snRNA-seq dataset allowed us to interrogate transcriptional differences across three heart failure subtypes: ICM, DCM, and HCM. To this end, we performed cell-type–specific DEG analysis by individually comparing DCM, HCM, and ICM against non-diseased donors, again identifying hundreds of DEGs for each disease subtype (Additional file 1 : Fig. S16d; Additional files 9 – 11 : Tables S8-10). This analysis revealed that most DEGs are distinct to each subtype. Using the Jaccard similarity to measure the degree of intersection between DEG sets, we find that the Jaccard similarity is low (< 0.5; values range from 0 for no intersection to 1 for total intersection) across disease subtypes (Additional file 1 : Fig. S16e) with greater similarity between the non-ischemic subtypes DCM and HCM. To better understand disease subtype convergence and divergence at the pathway level, we performed gene set enrichment analysis on genes ranked by decreasing expression fold change relative to non-diseased donors (Fig. 4 f). This analysis revealed that all heart failure subtypes are associated with well-characterized, shared molecular derangements such as decreased oxidative phosphorylation [ 70 ], which contrasts with the aging analysis that revealed many cell types with increased oxidative phosphorylation during aging (Additional file 1 : Fig. S13d). At the same time, we identified divergent pathways, such as the increase of interferon signaling specifically in HCM, increase of epithelial-mesenchymal transition in ICM, and decreased angiogenesis, interferon signaling, and androgen response in DCM (Fig. 4 f). Trajectory analyses of specific cardiac cell types further supported these observed enrichments during non-diseased to disease transitions (Additional file 1 : Fig. S17). For example, in cardiomyocytes, genes that increase along the trajectory from ND to HCM were associated with increased inflammation, while the reverse occurred in DCM (Additional file 1 : Fig. S17a). Although our diseased samples represent end-stage heart failure, where convergent signatures may obscure early disease-specific triggers, these results demonstrate that integrating multiple subtypes within a unified analytical framework reveals both shared and divergent molecular remodeling programs that would not be discernible in smaller, single-subtype datasets. Cardiac disease is associated with fetal reactivation across diverse cell types To explore how transcriptional and epigenetic changes across development, aging, and disease are interrelated, we performed an intersection analysis to quantify overlap between DEGs that change concordantly across pairs of contrasts (i.e., up- or down-regulated in the same direction) across pairs of biological contrasts (aging vs. development, development vs. disease, and aging vs. disease) for each cell type. For each comparison, we computed the fraction of co-directional DEGs relative to the total number of shared DEGs (Fig. 5 a) and evaluated statistical significance by generating an empirical null distribution given the observed DEG set sizes ( Methods ; Additional file 1 : Fig. S18a-b). Fig. 5. Open in a new tab Diverse cardiac cell types in disease subtypes display reactivation of fetal genes. a The Z-score for overlap of DEGs across all cell types for all pairs of contrasts (aging & development, aging & disease, disease & development). A threshold of significance (dotted line) indicates |Z-score|> 3. b The Z-score for overlap of DARs across all cell types for all pairs of contrasts. c Overrepresented pathways per cell type for DEGs that are shared between fetal and disease subtypes. d TFs per cell type that have significantly higher or lower GRN regulon activity in fetal and diseased states, relative to non-diseased hearts. e Cell–cell communications were inferred using liana x tensorcell2cell, which divides the overall cell–cell communication network in several programs and compares their relative contributions across biological states (fetal, non-diseased (ND), diseased). Program 5 has significant higher activity in fetal and diseased cells. f Cell type communication strength in Program 5. g The enriched pathways for each program, with positive enrichment indicating increased pathway activity. Program 5 is enriched for TGFβ signaling. h Major ligand-receptor interactions in program 5 that contribute to TGFβ signaling As expected, no significant overlap was observed between developmental and aging DEGs (Fig. 5 a). The overlap between aging and disease was marginally significant in cardiomyocytes, endocardial cells, fibroblasts, and myeloid cells (Fig. 5 a), suggesting that a subset of age-dysregulated genes may predispose to disease. In contrast, we observed strong and significant overlap between developmental and disease DEGs across nearly all cardiac cell types, both at a disease-binarized level (Fig. 5 a) and within each major heart failure subtype (DCM, HCM, and ICM) (Additional file 1 : Fig. S18c). Using our integrated snATAC-seq datasets, which included AMI/ICM donors, we also observed significant overlap between fetal and disease-enriched DARs (Fig. 5 b). These results strongly support the “fetal reactivation” hypothesis in heart disease, in which the stressed adult heart re-engages developmental gene programs [ 71 , 72 ]. Notably, while prior studies have largely focused on cardiomyocytes [ 73 ], our integrated cell-type–resolved multiomic analysis demonstrates that fetal reactivation is a pervasive, multi-lineage phenomenon, with around 20–40% of all disease-associated DEGs per heart failure subtype overlapping with fetal-enriched genes (Additional file 1 : Fig. S18d). Consistently, multiomics atlas-based results extend findings from a smaller study of pediatric heart failure [ 74 ] and prior bulk-tissue analyses [ 75 ]. The functional role of fetal reactivation in diseased adult hearts, whether adaptive or maladaptive, remains largely unresolved [ 71 , 72 ]. To further characterize fetal reactivation and its differences across disease subtypes, we identified DEGs shared between fetal and disease states and performed overrepresentation analysis (ORA) separately for co-upregulated and co-downregulated DEGs. We found that fetal reactivation gene sets were often shared among disease subtypes (Fig. 5 c). For instance, fetal-diseased co-upregulated genes were enriched for cell growth and hypertrophy-related pathways (e.g., E2F targets, epithelial-mesenchymal transition, apical junction), most prominently in fibroblasts and cardiomyocytes. In contrast, co-downregulated genes were enriched for inflammatory pathways (TNF-alpha, interferon gamma) across multiple cell types. Together, these findings suggest that advanced heart failure is characterized by partial reactivation of fetal-like transcriptional programs involving proliferative and hypertrophic activation in structural cell types, coupled with an anti-inflammatory, fibrotic polarization across non-myocyte populations [ 76 ]. Nomination of transcription factors driving fetal reactivation To identify the transcriptional regulators underlying the disease-associated fetal reactivation programs, we leveraged our integrated multiomic dataset to reconstruct cell type–specific gene regulatory networks (GRNs). We first identified TF-driven GRNs showing significant activity changes in diseased and fetal states relative to non-diseased adult hearts. We then focused on GRNs that were concordantly altered (either increased or decreased) in both fetal and diseased conditions. This analysis revealed distinct sets of TFs potentially mediating fetal reactivation in a cell type–specific manner (Fig. 5 d). Across non-cardiomyocyte cell types, thyroid hormone ( THRB ) was decreased in fetal and diseased states, reflecting the role of thyroid hormone in postnatal maturation and subclinical hypothyroidism observed in about one-third of heart failure patients [ 60 , 77 ]. In cardiomyocytes, we observed decreased activity of BACH2 , which has been previously associated with protection against cardiac hypertrophy [ 78 ]. Conversely, in fibroblasts, we identified upregulation of KLF12 , which has been shown to worsen fibrotic load after pressure overload induced heart failure in mouse models, and TEAD1 , which is implicated in fibroblast activation [ 79 ]. In endothelial cells, we identified upregulation of ZEB1 , a master TF critically involved in initiation and progression of atherosclerotic cardiovascular disease by regulation of endothelial cell and vascular smooth muscle cell functions [ 80 ]. Therefore, the fetal reactivation programs appear to involve a mixture of adaptive and maladaptive stress responses. Overall, these results nominate key TFs that may coordinate the reactivation of developmental gene programs across diverse cardiac cell types. Intercellular TGFβ signaling increases in both fetal and disease states Beyond cell-intrinsic regulation such as TF-driven GRNs, intercellular signaling plays a critical role in shaping cardiac cell states. To identify signaling pathways altered in fetal and diseased hearts, we also performed cell–cell communication (CCC) analysis, which infers ligand-receptor interactions based on statistically significant co-expression patterns. We applied liana x tensorcell2cell [ 81 ], a computational framework that integrates and harmonizes communication scores across multiple CCC inference methods [ 82 – 86 ] to compare cell–cell communication across biological states, focusing on which CCC programs differ between fetal and diseased hearts relative to ND hearts. Among the inferred CCC programs, program 5 had higher activity in both fetal and diseased hearts relative to ND hearts (Fig. 5 e). This program predominantly involved interactions among fibroblasts, pericytes, and vSMCs (Fig. 5 f; Additional file 1 : Fig. S19) and was highly enriched for TGFβ signaling (Fig. 5 g). Program 5 also featured increased communication through collagen-integrin interactions (Fig. 5 h), consistent with extracellular matrix remodeling and fibrotic activation. TGFβ plays dual roles in cardiac biology, promoting tissue maturation during development but driving fibrosis and adverse remodeling in disease [ 87 ]. Our cell–cell communication analysis thus reveals a shared upregulation of TGFβ-mediated intercellular signaling across both fetal and pathological hearts, suggesting that reactivation of developmental communication networks may underlie maladaptive remodeling in heart failure. Fetal reactivation is focalized in spatial niches We next explored whether fetal reactivation genes, previously examined only for a few candidate genes [ 88 ], are spatially localized within specific regions of diseased hearts or broadly distributed. To address this, we analyzed spatial transcriptomic datasets from recently published studies of non-diseased and AMI donors generated using the Visium platform from 10X Genomics (Fig. 6 a) [ 11 , 16 ]. Using the set of shared fetal + disease DEGs for each cell type, we computed a “fetal reactivation score” for each spatial spot that was weighted by the per-cell type fetal reactivation score and its spot-level cell type proportion computed by the cell2location method (Fig. 6 b). Specifically, this score represents the average expression of fetal reactivation genes minus the average expression of a reference set of other genes in the same corresponding expression quantile, with a positive score indicating enrichment of fetal reactivation genes in each spot. Across multiple diseased donors, we observed regions with positive fetal reactivation scores, which was often spatially distinct from the localization of gene activity scores obtained using either a random set of genes or non-fetal reactivated diseased DEGs (Additional file 1 : Figs. S20-28). Notably, fibrotic and ischemic regions exhibited a significantly higher proportion of spots with positive fetal reactivation scores compared to both ND donors and the remote, non-infarct regions in diseased AMI samples [ 16 ] (Fig. 6 c). Moreover, multiple donors displayed significant spatial focalization as measured by the spatial autocorrelation metric Moran’s I, suggesting localized spatial niches of strong fetal reactivation (Fig. 6 d). In AMI, these fetal reactivation spots were most enriched for cardiomyocytes and myeloid cells (Fig. 6 e), consistent with the increased proliferation of myeloid cells during acute cardiac injury. Fibroblasts were depleted from fetal reactivation scores, which may reflect differences between their transcriptional profiles and cell–cell communication during the acute infarct period and chronic heart failure. At the gene level, one of the most correlated genes with high fetal reactivation scores was NPPA , which has been well-established as a cardiac stress hormone (Fig. 6 f). Together, these analyses provide independent validation of fetal gene reactivation in human heart disease and reveal unexpected spatial focalization of disease-associated fetal reactivation within distinct pathological niches, particularly those associated with fibrosis and ischemia. Fig. 6. Open in a new tab Spatial patterns of fetal reactivation signatures in diseased human hearts. a Overview of spatial transcriptomics datasets analyzed. b Cell-type proportion weighted fetal reactivation scores (left) and abundance scores of five major cell types (right: Cardiomyocyte – yellow, Fibroblast – red, Myeloid – blue, Endothelial – green, and Pericyte – purple) in representative non-diseased (ND) and diseased regions. Regions of fetal reactivation in diseased hearts are highlighted by arrow heads. c Proportion of fetal reactivation spots in ND and diseased regions. d Spatial localization, as calculated by Moran’s I, in ND and diseased regions. e Cell type enrichment in fetal reactivation spots, calculated using a model that assigns positive coefficients to cell types that have higher proportion in spots with positive fetal reactivation scores. f Genes with strongest spot-level correlation to fetal reactivation score Cell-type–specific multiomic atlas refines genetic risk loci to target gene linkages Finally, we leveraged our integrated single-nucleus multiomics atlas to improve the functional interpretation of heart failure-associated loci identified from recently reported genome-wide association studies (GWAS). Across traits, the majority of GWAS risk loci (> 90%) are in non-coding regions of the genome, where causal variants often modulate cell-type–specific expression of downstream genes through altered activity of cis -regulatory elements. Recently, three large-scale GWAS have uncovered ~ 50 new genetic loci associated with DCM [ 89 , 90 ] and HCM [ 91 ]. These studies primarily prioritized downstream genes through expression quantitative trait loci (eQTL) or rare protein-coding rare variant co-localization, without incorporating cell-type–specific chromatin accessibility data. We hypothesized that intersecting GWAS loci with cell-type–specific enhancer-to-gene linkages predicted from combined snRNA-seq and snATAC-seq datasets could refine target gene identification and determine the relevant cardiac cell types. Leveraging the multiomic non-diseased datasets, we identified peak-to-gene linkages using scE2G [ 92 ], a computational method that has recently been shown to outperform alternative methods in recovering experimental validated enhancers (Fig. 7 a). Performing scE2G identified ~ 70,000 to 100,000 peak-to-gene linkages per each cell type, including ubiquitous and cell-type–specific linkages (Fig. 7 b-c; Methods ). We next intersected these linked peaks with an expanded SNP set containing both the lead variants (i.e., DCM and HCM GWAS SNPs) and other variants in high linkage disequilibrium (r 2 > 0.8) to these lead SNPs (Fig. 7 d). Through integration with this expanded SNP set, we identified 78/117 DCM (67%) and 58/68 HCM loci (85%) overlapping peak-to-gene linkages, consistent with previously reported overlap rates of cardiovascular GWAS loci in accessible regions [ 93 ]. For each locus, we then defined the epigenomic prioritized gene(s) as the top two genes in terms of enhancer-to-gene confidence score ( Methods ). When comparing our results with the candidate gene(s) reported in the original studies using orthogonal approaches, we observed strong concordance, with 81% (63/78) full or partial concordance for DCM loci and 72% (41/58) for HCM loci (Fig. 7 e; Additional file 1 : Figs. S29-30; Additional files 12–13: Tables S11-12; Methods ). In several cases, our analysis provided cell-type–specific enhancers providing support for cell type(s) that may be relevant for disease association. For example, the lead SNP chr6:54,092,790, which was associated with both DCM and HCM, is located within a cardiomyocyte-specific peak linked to MLIP , which is important for cardiac stress adaptation [ 94 ] (Fig. S29). Fig. 7. Open in a new tab Cell-type–specific enhancers overlap with GWAS loci. a We used scE2G to identify cell-type–specific peak-to-gene linkages from multiomic nuclei in our atlas. b Number of peak-to-gene linkages identified for each cell type. c Clustering of peak-to-gene linkages, which reveals some linkages that are shared by many cell types and other linkages that are cell type-specific. d Approach for intersecting GWAS variants with peak-to-gene linkages. e Concordance rate between the linked genes using our multiomic epigenomic approach versus the prioritized gene linked by the original GWAS, which did not use single-cell epigenomic data. f rs4616, located in the 3’ UTR of TRPC4AP , is a lead DCM GWAS SNP that lies within a cardiomyocyte-specific enhancer linked to MYH7B . GSS was the prioritized gene by the original GWASs, but was not linked in our analysis. g Cell-type–specific expression of genes near rs4616 ( GSS , MYH7B , and TRPC4AP ). h SNP chr1:211,933,964 is linked to DTL in original GWAS and using our epigenomic method, based on highest scE2G (linkage prediction) score. i Donor-pseudobulked expression of DTL in cardiomyocytes (top) and accessibility of the peak containing chr1:211,933,964 (bottom). j SNP chr15:78,760,672 is linked to ADAMTS7 in original GWAS and using our epigenomic method. k Donor-pseudobulked expression of ADAMTS7 in endothelial cells (top) and accessibility of the peak containing chr15:78,760,672 (bottom). Abbreviations: DCM: dilated cardiomyopathy; GWAS: genome-wide association study; HCM: hypertrophic cardiomyopathy; SNP: single nucleotide polymorphism; UTR: untranslated region For cases where our epigenomic approach is discordant with the prioritized gene by other methods, one reason is that the causal SNPs may not act through non-coding regulatory mechanisms. In other cases of discordance, our multiomic framework identifies more biologically plausible downstream target genes. For example, both DCM GWAS linked SNP rs4616, located in the 3’ untranslated region of TRPCA4P , to glutathione synthetase ( GSS ), a broadly expressed gene involved in redox signaling (Additional file 1 : Fig. S31a). In contrast, we identify that rs4616 lies within a cardiomyocyte-specific peak that is most strongly linked to MYH7B (Fig. 7 f). Unlike GSS , MYH7B encodes a muscle-restricted, non-translated myosin isoform (Fig. 7 g; Additional file 1 : Fig. S31b) previously implicated in dilated cardiomyopathy through non-coding mechanisms [ 95 , 96 ]. Interestingly, some of the genes that are predicted as regulated by disease-associated SNPs are also fetal reactivation genes. For example, DTL , which has been implicated in DNA damage repair as a regulator of TP53 [ 97 ], is prioritized by the original GWAS and our epigenomic method as linked to the DCM-associated SNP chr1: 211,933,964 (Fig. 7 h). Integrating the scE2G predictions with differential expression and chromatin accessibility analysis, we demonstrate that this gene is specifically upregulated in both fetal and DCM hearts (Fig. 7 i). Moreover, the chromatin accessibility of fetal hearts for the putative enhancer containing this SNP is more accessible in fetal hearts. Another example is the HCM-associated SNP chr15:78,760,672, which is linked to ADAMTS7 and has higher fetal chromatin accessibility and higher expression in fetal and HCM compared to non-diseased endothelial cells. Further work will be required to probe chromatin accessibility of these loci in DCM and HCM hearts, to causally perturb these SNPs or SNP-overlapping peaks to confirm the effect on the nominated genes, and to investigate how overexpression of these genes in specific cell types affect cardiac development and disease progression. Together, these results demonstrate that our multiomic enhancer-to-gene linkage framework can more robustly identify causal target gene and enhance the biological interpretation of cardiac GWAS loci by integrating cell type–specific regulatory contexts. Discussion Here, we present a large, integrated single-nucleus transcriptomic and epigenomic analysis of the human cardiac cell-types to date. In addition to aggregating publicly available datasets, we generated 19 new snRNA-seq and 11 snATAC-seq datasets from human donors spanning the full life span, establishing a comprehensive multimodal resource for the cardiac research community. This large-scale integration enabled us to address key challenges associated with dataset size, batch heterogeneity, and donor variability, thereby substantially improving statistical power. Nonetheless, certain confounding factors inherent to human donor studies, such as postmortem interval, tissue sampling strategy, and inclusion criteria, cannot be fully controlled. Additionally, while we performed ambient RNA removal for datasets in which we had both raw and filtered counts, for some studies, we do not have access to the raw count matrix or raw sequencing data due to patient identifiability constraints. While ambient RNA was already performed by the individual studies, the lack of the same ambient RNA removal steps may have resulted in inadequate ambient RNA removal, which is a major limitation of the current snRNA-seq method. Furthermore, genetic variation among donors remains an important potential confounder. While most datasets currently lack genotyping information, future efforts combining single-cell multiomics with genotype-resolved eQTL analyses will enable direct investigation of genetic variants on gene regulation at cell-type resolution [ 98 ], complementing prior bulk eQTL studies. By systematically comparing cell-type-resolved transcriptional and chromatin accessibility profiles across sex, aging, development, and disease, we uncovered distinct patterns of transcriptional and epigenetic remodeling. Sex- and age-associated transcriptional differences were modest, whereas developmental and disease-associated changes were extensive, involving hundreds to thousands of genes per cell type. Surprisingly, we identified that some inflammatory pathways such as TNF-alpha signaling are decreased in heart failure subtypes, most notably DCM. Whether this reflects a true shift toward an anti-inflammatory state in end-stage heart failure or instead represents immune dysregulation rather than active suppression will require further investigation. The limited overlap between age- and disease-associated genes suggests that aging and disease may represent molecularly distinct trajectories, with aging reflecting gradual, compensatory remodeling and disease reflecting acute, end-stage perturbations. Despite including 299 donors, the number of detected aging-associated DEGs did not plateau, implying that even larger cohorts are required to capture subtle aging signatures across diverse cardiac cell types. In contrast, we found strong overlap between developmental and disease-associated transcriptional programs, supporting the concept of fetal gene reactivation in heart failure. Through integrative analysis of cell-type–resolved chromatin accessibility and gene expression profiles, we identified TF-driven GRNs and cell–cell communication pathways, notably TGFβ signaling, that are concordantly regulated in both fetal and diseased hearts. However, given the relatively modest accuracy of computational methods for predicting gene regulatory and cell–cell communication networks, future work involving experimental perturbations of these intercellular interactions and TFs in a cardiac context will be critical. Furthermore, re-analysis of orthogonal spatial transcriptomics datasets revealed that fetal reactivation is not spatially uniform but rather focalized within distinct pathological niches such as fibrotic and ischemic regions. Whether these fetal reactivation foci represent adaptive, pro-regenerative microenvironments or maladaptive remodeling hubs remains an open question warranting experimental validation. Additionally, we explored whether we could leverage our single-nucleus multiomics atlas to provide further support that links GWAS DCM and HCM loci to downstream target genes. We demonstrate that integration of snRNA-seq and snATAC-seq provides more robust evidence often corroborating the genes prioritized using gene expression only approaches and nominating involved cell types. Importantly, our analysis identifies more biologically plausible downstream genes. Still, further predictive modeling and experimental work is required for identifying which SNPs are causal and validating their downstream targets. Overall, our study provides new insights into the molecular logic of human cardiac remodeling and establishes a foundation for future mechanistic studies. While the direct human relevance of these data is a major strength, we only analyzed nuclear transcriptomes, which does not account for post-transcriptional, translational, and metabolic regulation outside the nucleus. Finally, our findings are correlative and would require experimental validation to establish causality. Model organisms may enable some of these follow-up studies; however, differences in cardiac biology across species may limit their applicability. As protocols for hiPSC-derived cardiac organoids continue to advance, they offer a promising avenue to experimentally test hypotheses generated by our atlas. Thus, coupling this integrated multiomic resource with targeted mechanistic and perturbation studies will deepen our understanding of how the human heart develops, remodels, and fails across its full phenotypic spectrum. Conclusions In summary, this work delivers an integrated multiomic framework for studying human cardiac biology across development, aging, and disease. Leveraging the large-scale multiomic resource, we systematically uncovered differentially expressed genes and gene regulatory networks that contribute to this phenotypic diversity. Among many other findings, our analysis revealed strong reactivation of fetal programs across nearly all cell types in the diseased human heart, which we further characterized in epigenomic, cell–cell communication, and spatial contexts. Finally, we generated a comprehensive, cell-type–specific enhancer-to-gene regulatory map that further prioritizes downstream genes as mediators of dilated and hypertrophic cardiomyopathy. Overall, this work provides the most extensive multiomic resource to date for studying human cardiac biology and provides a groundwork for further mechanistic and translational studies aimed at understanding and treating cardiac disease. Methods Human cardiac tissue procurement Newly generated snRNA-seq and snATAC-seq were prepared from the adult hearts deposited in the Penn Heart Tissue Biobank. Hearts were harvested at time of cardiac surgery from subjects with heart failure (2 DCM patients) or from non-diseased donor hearts that were unable to be transplanted within the surgical window. Hearts were transported to the research lab in cardioplegic solution to prevent ischemic damage, and single nucleus suspensions of tissue were biopsied from the left ventricular wall. Written consent for research use of these tissues was obtained from the subject providing tissue or next of kin and approved by the University of Pennsylvania institutional review board. Fetal human hearts (8–18 weeks post fertilization) were obtained from donors undergoing elective abortion at the University of Pennsylvania, with all donors providing written consent for the use of the conceptus in research. Fetal age was determined via ultrasonographic measurement of crown-rump length or head circumference. Post-conception week estimates were obtained from obstetric records and tissue pathology. As an orthogonal metric, we performed ultrasonographic measurement of crown-rump length or head circumference. The heart tissues were submerged in RPMI-1640 medium and meticulously dissected under a stereomicroscope to remove surrounding connective and adipose tissues. Human cardiac nuclei isolation and fluorescence activated nuclear sorting (FANS) Human cardiac tissues were flash-frozen in liquid nitrogen and stored at −80 °C before nuclear isolation. Nuclei were isolated and purified as previously described with some modifications [ 6 ]. Briefly, 14 ml sucrose cushion 1.8 M sucrose (Sigma-Aldrich, RNase and DNase free, ultra-pure grade), 10 mM Tris–HCl pH 8.0, 3 mM MgAc 2 (Sigma-Aldrich), and protease inhibitor cocktail (Roche) were added to the bottom of centrifuge tubes (Beckman Coulter). Using a glass homogenizer (Wheaton), the tissue was dounce homogenized (19 times with loose pestle followed by 4 times with tight pestle) in 12 mL homogenization buffer (0.32 M sucrose, 5 mM CaCl 2 (Sigma-Aldrich), 3 mM MgAc 2 , 10 mM Tris–HCl pH 8.0, 0.1% Triton X-100 (Sigma-Aldrich), 0.1 mM EDTA (Invitrogen), protease inhibitor cocktail (Roche). Homogenates (~ 12 ml) were carefully layered onto the sucrose cushion in the centrifuge tubes, and 10 mL homogenization buffer was added atop the homogenates. The tubes were then balanced and centrifuged in a Beckman Coulter L7-65 Ultracentrifuge at 82,705 g at 4 °C for 120 min using a Beckman Coulter SW28 swing bucket rotor (Beckman Coulter). The supernatant was carefully removed via aspiration. 1 mL chilled DPBS was added to resuspend the nuclear pellet, and nuclei were subsequently transferred to a 1.5-ml tube. Nuclei were pelleted at 500 g for 10 min at 4 °C and then resuspended in 0.01% BSA (Sigma-Aldrich) in DPBS. After resuspension, nuclei were filtered through a 40-μm cell strainer (Fisher Scientific), visually inspected for morphology and quality assurance, and counted using a Fuchs-Rosenthal counting chamber before FANS. Isolated nuclei were resuspended in 1.5 mL blocking buffer (1 × PBS, 0.5% BSA (Sigma A4503), Roche Complete Protease Inhibitor without EDTA), and blocked for 20 min at 4 °C with rotation. After blocking, nuclei were incubated with DAPI (1:1000) for 5 min. Nuclei were then washed, pelleted, and resuspended in 1 mL FACS buffer (1 × PBS, 1% BSA (Sigma A4503), Roche Complete Protease Inhibitor without EDTA), and passed through a 35 µm strainer (Corning #352,235) in preparation for flow cytometry. FANS was performed using a BD Biosciences Influx cell sorter at the University of Pennsylvania Flow Cytometry and Cell Sorting Facility. Singlet nuclei were identified using DAPI fluorescence. Single-nucleus RNA sequencing (sNucDrop-seq) Single nucleus RNA sequencing analysis of FANS human cardiac nuclei was performed using the sNucDrop-seq protocol [ 6 , 99 ]. Briefly, FANS nuclei were individually diluted to a concentration of 100 nuclei/mL in DPBS containing 0.01% BSA. Approximately 1.25 ml of this single-nucleus suspension was loaded for each sNucDrop-seq run. The single-nucleus suspension was then co-encapsulated with barcoded beads (ChemGenes) using an Aquapel-coated PDMS microfluidic device (μFluidix) connected to syringe pumps (KD Scientific) via polyethylene tubing with an inner diameter of 0.38 mm (Scientific Commodities). Barcoded beads were resuspended in lysis buffer (200 mM Tris–HCl pH8.0, 20 mM EDTA, 6% Ficoll PM-400 (GE Healthcare/Fisher Scientific), 0.2% Sarkosyl (Sigma-Aldrich), and 50 mM DTT (Fermentas; freshly made on the day of run) at a concentration of 120 beads/ml. The flow rates for nuclei and beads were set to 4,000 ml/hour, while QX200 droplet generation oil (Bio-rad) was run at 15,000 mL/hr. A typical run lasts 20 min. Droplet breakage with Perfluoro-1-octanol (Sigma-Aldrich), reverse transcription and exonuclease I treatment were performed, as previously described, with minor modifications. For up to 120,000 beads, 200 μl of reverse transcription (RT) mix (1 × Maxima RT buffer (ThermoFisher), 4% Ficoll PM-400, 1 mM dNTPs (Clontech), 1 U/mL RNase inhibitor, 2.5 mM Template Switch Oligo (TSO: AAGCAGTGGTATCAACGCAGAGTGAATrGrGrG), and 10 U/mL Maxima H Minus Reverse Transcriptase (ThermoFisher) were added. The RT reaction was incubated at room temperature for 30 min, followed by incubation at 42 °C for 120 min. After determining an optimal number of PCR cycles for amplification of cDNA using real-time PCR analysis (Applied Biosystems QuantStudio 7 Flex), aliquots of 6,000 beads were amplified by 50-μl PCR reactions (25 μl of 2 × KAPA HiFi hotstart readymix (KAPA biosystems), 0.4 μl of 100 mM TSO-PCR primer (AAGCAGTGGTATCAACGCAGAGT), 24.6 μl of nuclease-free water) with the following thermal cycling parameter (95 °C for 3 min; 4 cycles of 98 °C for 20 s, 65 °C for 45 s, 72 °C for 3 min; 10–12 cycles of 98 °C for 20 s, 67 °C for 45 s, 72 °C for 3 min; 72 °C for 5 min, hold at 4 °C). After two rounds of purification with 0.6 × SPRISelect beads (Beckman Coulter), amplified cDNA was eluted with 10 μL of water. We then tagmented cDNA using the Nextera XT DNA sample preparation kit (Illumina, FC-131–1096), starting with 550 pg of cDNA pooled in equal amounts, from all PCR reactions for a given run. Following cDNA tagmentation, we further amplified the tagmented cDNA libraries with 12 enrichment PCR cycles using the Illumina Nextera XT i7 primers along with the P5-TSO hybrid primer. After quality control analysis by Qubit 3.0 (Invitrogen) and a Bioanalyzer (Agilent), libraries were sequenced on an Illumina NextSeq 500 instrument using the 75-cycle High Output v2 Kit (Illumina). We loaded the library at 2.0 pM and provided Custom Read1 Primer (GCCTGTCCGCGGAAGCAGTGGTATCAACGCAGAGTAC) at 0.3 mM in position 7 of the reagent cartridge. The sequencing configuration was 20 bp (Read1), 8 bp (Index1), and 60 bp (Read2). Single-nucleus ATAC-seq (10 × Genomics) Using the nuclei isolation protocol detailed for sNucDrop-seq protocol, nuclei were aliquoted for single nucleus ATAC sequencing analysis of human cardiac nuclei. snATAC-seq was performed using the 10X Genomics ATAC v1 kit. A target of 10,000 nuclei was incubated with Tn5 transposase per reaction. After tagmentation, nuclei were encapsulated using Chromium Chip E in the manufacturer’s Chromium Controller. The library was sequenced on an Illumina NextSeq 500 using these read parameters: 50 bp for Read1, 8 bp i7 index, 16 bp i5 index, and 50 bp for Read2. Penn snRNA-seq dataset preprocessing Raw FASTQ files were aligned against the hg38 reference genome (refdata-cellranger-arc-GRCh38-2020-A-2.0.0) using STARSolo (STAR v2.7.11a), using the parameters –soloUMIlen 8 –soloCBlen 2 –soloFeatures Gene GeneFull –soloCellFilter EmptyDrops_CR (Additional file 1 : Fig. S2a). The raw and filtered (based on the UMI threshold set by EmptyDrop_CR) directories produced were passed as inputs to SoupX (v1.6.2) for ambient RNA removal. Doublets were removed using scrublet (v1) by removing cells with predicted_doublet = = True. The filtered, ambient RNA removed count matrices were concatenated together and analyzed together using scanpy (v1.10.2). We performed scVI-based integration using donor_id as the batch key. Using the KNN graph in the scVI embedding, we called Leiden clusters with a resolution of 0.5 and annotated clusters based on marker genes. The primary marker genes that assisted in this manual annotation included known markers reported in the literature [ 7 , 8 , 16 ] including MYH7/RYR2/TTN (Cardiomyocyte); COL3A1/COL1A1/DCN/LUM (Fibroblast); ADIPOQ/CIDEX/DGAT2/FABP4/PPARG (Adipocyte); ECMN/PECAM1/VWF (Endothelial); POSTN + other endothelial markers (Endocardial); ITLN1 / MSLN/WT1 (Epicardial); MMRN1/NRP2/PROX1 (LEC); MS4A1/SKAP1/THEMIS (Lymphoid); CPA3/HDC/KIT (Mast); CD163/F13A1/MS4A6A (Myeloid); NRXN1/NRXN3/XKR4 (Neuronal); DLC1/EGLAM/PDGFRB/RSG5 (Pericyte); ACTA2/CARMN/LMOD1/MYH11 (vSMC). More information about the Penn snRNA-seq dataset is provided in Fig. S2. Penn snATAC-seq dataset preprocessing Raw FASTQ files were aligned against the hg38 reference genome (refdata-cellranger-arc-GRCh38-2020-A-2.0.0) using cellranger-atac 2.1.0 under default parameters to produce fragment files (Fig. S5a). These fragment files were inputted to SnapATAC2 (v2.8.0) for preprocessing, which included calculating the number of fragments and TSS enrichment score per barcode. We filtered cells to those with at least 1000 unique fragments and a TSS enrichment score > 3 and combined all samples together. We then performed harmony (python v0.0.10) batch integration (using donor_id as the batch key), followed by SnapATAC2’s spectral embedding. Leiden clusters were identified using a resolution of 0.5. Cluster annotation was performed using MAGIC (v3.0.0) imputed gene expressions of marker genes noted in the snRNA-seq section above. Peak calling was performed using MACS3 (v3.0.2). More information about the Penn snATAC-seq dataset is provided in Additional file 1 : Fig. S5. External snRNA-seq dataset processing External datasets [ 8 , 11 , 15 – 21 , 24 ] were downloaded according to the data availability sections of each paper (Additional file 1 : Fig. S3a). Most of these studies only sampled the LV, which was the focus of this study. For studies that included other datasets such as [ 8 ] and [ 16 ], we only included datasets sampled from the LV. The link for downloading the processed data is provided in Additional file 2 : Table S1 under the “RNA” section. The accessions of the 62 ENCODE v4 snRNA-seq datasets we downloaded are provided under the “ENCODE RNA” section [ 100 ]. snRNA-seq integration and broad cell type annotation The raw count matrices (filtered to remove droplets containing non-cells) across all the studies were concatenated together using Scanpy [ 101 ]. For ambient RNA removal for our new dataset and published datasets in which unfiltered count matrices were readily available, we performed ambient RNA removal using SoupX [ 102 ]. For datasets with only processed filtered count matrices, these count matrices were used, as ambient RNA removal was already performed by each study using either SoupX or cellbender (v0.2.2) [ 103 ]. Doublets were removed using scrublet (v1) by removing cells with predicted_doublet = = True. For sNucDrop-seq data, we kept nuclei with greater than 100 unique molecular identifiers (UMIs) due to the lower library complexity and sequencing depth of these samples. For 10X Genomics datasets, we kept all nuclei with greater than 500 UMIs. For postnatal nuclei, we removed those with > 5% mitochondrial transcripts, > 5% ribosomal protein transcripts, or > 5% hemoglobin transcripts. For fetal nuclei, we removed those with > 15% mitochondrial transcripts, > 15% ribosomal protein transcripts, or > 5% hemoglobin transcripts (Additional file 1 : Fig. S3a). After benchmarking the performance of 11 integration methods using scIB (v1.1.6) [ 13 ], we used the top-performing method scVI (v1.2.0) [ 27 ] to remove batch effects at the UMAP embedding level in two rounds: whole atlas to define broad cell types, and then at a cell-type– specific level for subclustering. For the first scVI round, we trained for a maximum of 10 epochs. The layer of the anndata object used was the raw counts. The other parameters were all the defaults (n_layers = 1, dropout_rate = 0.1, n_latent = 10). The batch size for training used the default scVI settings, which scales the batch sizes for training depending on the size of the dataset. “tech_plus_study” was specified as the batch and “donor_id” was included as a categorical covariate. “pct_counts_ribo”, “pct_counts_mt”, and “pct_counts_hb”, and log10_UMI were used as continuous covariates in our model. Donor age was not included as a covariate in the latent space, as we considered this a biological signal that we did not want to regress out. scVI batch corrected counts for all genes were stored, using the ENCODE v4 study as the conditional batch as it had the largest number of non-diseased donors. These batch-corrected counts were used for batch-corrected analysis with liana x tensorcell2cell (v1.2.1) [ 81 ], as recommended. Raw counts were used for differential gene expression, as batch corrected counts have been shown to distort differential gene expression analysis [ 12 ]. In the scVI embedding, leiden clustering at a resolution of 0.5 was performed using Scanpy. To avoid annotating subclusters that may arise from technical rather than biological differences, Leiden clusters in which > 80% of cells came from a single study were not annotated and removed. After this filtering, the top 50 marker genes per Leiden cluster, including marker genes noted above in “Penn snRNA-seq dataset preprocessing”, were used to independently annotate each Leiden cluster, agnostic to the cell annotations from the original studies. For each of these cell types, we re-examined the UMAP embedding and confirmed marker gene expression for all cells in those leiden clusters. In this cell-type–filtered count matrix, cells that segregated away from the main cluster and/or expressed marker genes for other cell types may represent doublets/high ambient RNA cells, so we removed these cells. At this stage, atrial cardiomyocytes, which separated from ventricular cardiomyocytes, were removed, as they arose mostly from fetal hearts and only a few studies. Cell type subclustering for snRNA-seq datasets After this first round of integration, we subsetted the count matrix to each individual cell type and reran scVI using 500 epochs (with early stopping allowed using elbo_validation) per cell type. All cell types converged within 500 epochs. For this cell-type–specific scVI integration, we added three additional covariates (“cm_score”, “fib_score”, and “endo_score”) which reflect 10 highly-expressed genes specific to these three cell types. The “cm_score” used TNNT2, TTN, MYH6, MYH7, ACTC1, TNNI3, MYL2, MYL7, RYR2, and PLN. The “endo_score” used PECAM1, VWF, KDR, FLT1, ESAM", "EMCN, RAMP2, CD34, ENG, and PLVAP. The “fib_score” used COL1A1, COL1A2, COL3A1, DCN, LUM, COL5A1, COL6A1, COL6A2, COL6A3, and FBLN1. We found that added these additional covariates, which were calculated using scanpy’s score_genes(), successfully removed artifactual subcluster levels of ambient RNA (most from cardiomyocytes, rather than fibroblasts or endothelial cells) in less common cell types. We performed subclustering using the Leiden resolution of 0.25, as this was associated with highest Silhouette score, for the eight less common cell types. For the five major cell types (Cardiomyocyte, Endothelial, Fibroblast, Myeloid, Pericyte), we used a higher resolution of 0.5. Using the top 50 marker genes per cell type, we performed a literature search to identify genes with clear functions ascribed within the specific cell types/cell states. After this subclustering step, we again removed additional subclusters that appeared to be doublets. Finally, we combined all filtered cell type count matrices to produce a final dataset of 2,275,105 nuclei across 13 major cardiac cell types. To assess whether cell-type–specific subclusters were enriched or depleted for specific metadata categories (age group, developmental and disease status), we performed a chi-square-based enrichment test. Briefly, the global frequencies of each metadata category were computed across all cells. This was compared against the observed counts of cells belonging to each metadata category. Expected counts were derived assuming random distribution of cells based on global frequencies and subcluster size. A chi-squared goodness-of-fit test was applied to each subcluster to evaluate whether the observed distribution of metadata categories deviated significantly from expectation. Enrichment was quantified using the log2 fold change between observed and expected. The adjusted p -values after FDR correction were almost always significant due to the large number of cells, so our enrichment analysis focused on the effect size (log change). External snATAC-seq dataset processing Fragment files from external datasets [ 10 , 11 , 16 , 24 ] were downloaded according to the data availability sections of each paper (Additional file 1 : Fig. S6b). The link for downloading the processed data is provided in Additional file 2 : Table S1 under the “ATAC” section. The accessions of the 71 ENCODE v4 snATAC-seq datasets that we downloaded are provided under the “ENCODE ATAC” section [ 100 ]. snATAC-seq integrative preprocessing We combined Penn fragment files with previously published fragment files from the other datasets as input to SnapATAC2 [ 33 ]. Using SnapATAC2, we performed quality control filtering by TSS enrichment score and number of fragments per cell were used to define per-sample cutoffs. After combining high quality cells across samples, we performed spectral embedding on the tile count matrix (5 kb tiles across the genome) and Harmony integration. Leiden clustering with resolution of 1 was performed using the spectral embedding nearest neighbor graph. Cell type annotation of these clusters was performed using two approaches. First, gene expression using imputed using MAGIC [ 104 ] and marker genes for snRNA-seq were used to assist with cell type annotation. While this approach worked for more common cell types, we found the second approach, which involved label transferring of multiome cells, to be more effective at annotating smaller clusters with rarer cell types (e.g., adipocyte, epicardial) (Additional file 1 : Fig. S6c). In this second approach, we filtered the snATAC-seq nuclei to those produced with 10X multiome and examined the proportion of snRNA-seq annotations for each leiden cluster. Leiden clusters with > 90% of cells sharing a single snRNA-seq annotation were annotated as that cell type. For clusters with > 20% of a rare cell type, these leiden clusters were kept and then marker gene imputation with MAGIC was performed to identify likely true cell types. All other leiden clusters were removed, as they may represent doublets. The inability to identify distinct clusters for two cell types present in snRNA-seq (LEC, Endocardial) may be due to less pronounced epigenetic vs. transcriptomic differences for LEC and Endocardial cells, which were located within the Endothelial cluster. After cell type annotation, reads were aggregated at the cell type level for peak calling with MACS3 [ 36 ], producing a cell by peak count matrix. Harmony batch integration and UMAP dimensionality reduction was reperformed on this matrix with peaks as features to generate the final UMAP plot in Fig. 2 b. Intersection of snATAC-seq peaks with other feature sets The cell-type–specific peaks called by MACS3 were merged into 500 nucleotide regions using SnapATAC2. The peak coordinates were converted to BED format [ 105 ]. The Spurrell et al. 2022 list with 33,317 enhancers [ 2 ], called from healthy donors was also converted to BED format. For intersection of Spurrell enhancers and snATAC-seq peaks, snATAC-seq peaks were filtered to only those that are located > 1 kb from a transcriptional start site. This filtered set of distal peaks was used for intersection with Spurrell enhancers, as this was used for their definition of enhancers. The distance from the nearest TSS was computed using ChIPSeeker [ 106 ]. Due to the significant length discrepancy between snATAC-seq peaks (fixed length of 500 nt) and Spurrell enhancers (median length 4244 nt), we calculated the number of Spurrell enhancers with at least one snATAC-seq overlapping it, using: bedtools (v2.31.1) intersect -a < Spurrell_bed > -b < snATAC_enhancer_peaks_bed > -c. Using bedtools shuffle, we generated shuffled snATAC-seq peak sets and then used bedtools intersect to obtain a null distribution. Differential expressed gene (DEG) analysis We tested three differential methods for differential expressed gene (DEG) and differential accessibility region (DAR) analysis: (1) pseudobulked fixed effect model using all donors with “batch as covariate” with DESeq2 (v1.44.0), (2) pseudobulked fixed effect model on individual studies using DESeq2 followed by effect size Weighted Fisher’s meta-analysis, (3) mixed effect model using NEBULA (v1.5.2). After demonstrating that approach 1 captures the most “ground truth” DEGs when performing sex and aging DEG analysis, we performed all analyses using approach 1. For approach 1 (“batch as covariate”), we generated donor-based pseudobulked count matrices, by aggregating the raw counts of all nuclei of the same cell type annotation for each donor, using our revised cell type annotation. This produces a donor x gene count matrix. We included only donors that have at least 50,000 total transcripts, since this reliably had a Spearman correlation > 0.8 for all cell types (Additional file 1 : Fig. S8a-b). We then performed differential gene expression with pydeseq2 (v0.5.2) [ 44 ], a faster and more scalable python-based implementation of DESeq2 [ 43 ], using this design matrix: expression ~ sex + age_group + disease_binary + tech_plus_study. We opted to include age as a categorical variable rather than a continuous variable as all other variables were categorical and to account for non-linear changes with age. The age groups therefore included fetal, young age (≤ 40), middle age [ 40 – 59 ], and old age (≥ 60). Age cutoffs reflect transition periods of age-related cardiovascular disease risk. As several of the different types of cardiac disease are represented only by one or two studies, the disease status was initially binarized to either diseased or non-diseased. For sex and aging analysis, we did not include diseased donors, so this covariate was not included in the design matrix (expression ~ sex + age_group + tech_plus_study). Many different droplet-based technologies were used for snRNA-seq data generation, including Dropseq, 10X Genomics 3’ v1-v3, and 10X Multiome v1. As some studies used multiple different technologies, we created a variable called “tech_plus_study” that combines both the study and technology information together. Without regressing out this technical batch effect, the pseudobulked PC1 and PC2 show significant correlation with tech_plus_study. Using limma (v3.60.6) [ 107 ] to remove batch effects at the pseudobulked level reduced the clear segregation of samples by technology and study, leading to segregation by development stage along PC2 (Additional file 1 : Fig. S9a-b). After performing binarized disease DEG analysis, we then compared three disease subtypes with the most representation (DCM, HCM, and ICM) against non-diseased (ND) donors using this design matrix (expression ~ sex + age_group + non_binarized_disease + tech_plus_study). DESeq2 returns a log 2 fold change (FC) and false discovery rate (FDR) [ 108 ] adjusted p -value. We considered significant upregulated genes to be those with log 2 FC > 0.5 (corresponding to 1.41-fold change) and adjusted p < 0.05. We considered significant downregulated genes to be those with log 2 FC < −0.5 (0.707-fold change) and adjusted p < 0.05. The interpretation of the “up” and “down” depends on the contrast. The contrasts we focused on for this study include sex: male vs. female, development: fetal vs. young (age group), aging: young vs. old (age group), and disease: Y vs. N (disease_binary), as indicated in Fig. 2 a. As an example, an “upregulated” gene in the developmental contrast (fetal vs. young) would represent a gene significantly higher in fetal donors (group 1 of contrast) compared to the postnatal young donors (group 2 of contrast). A “downregulated” gene would represent a gene significantly higher in the postnatal young donors (or equivalently a gene down in fetal donors). For approach 2 (“meta-analysis”), we performed DESeq2 for each study individually, rather than including “tech_plus_study” as a covariate (Additional file 1 : Fig. S9c). Therefore, our DESeq2 design matrix for sex and aging analysis is expression ~ sex + age_group. To obtain reliable effect sizes and p -values, we only included studies that had at least three donors for each contrast category (e.g., three males and three females). Importantly, this reduced the overall number of donors that we could include when compared to approach 1. We obtained an aggregate fold change and p -value for each gene as follows. First, the fold changes were computed per gene for each study individually. Then, using the log fold change standard error (lfcSE), we weighted each study inversely (1/lfcSE 2 ). The aggregate effect size was calculated the sum of each individual effect size (log 2 FC) divided by the sum of the weights. The standard error of this estimate was calculated as the square root of the 1/(sum of weights). The z-statistic was calculated as the effect size divided by the SE, and a p -value was obtained as p meta = 2 x [1- Φ (|z meta |)] where Φ is the inverse cumulative density function for the normal distribution. Finally, this meta-analysis p -value was adjusted using the Benjamini–Hochberg method to correct for multiple hypothesis testing. To call DEGs, we used the same cutoffs of abs(log 2 FC) > 0.5 and p meta-adj. < 0.05 as approach 1. This approach was more conservative than approach 1 and resulted in slightly lower sensitivity at similar specificity, so we did not continue with this approach (Additional file 1 : Fig. S9d-k). For approach 3 (“mixed-effect”), the input was the single-cell count matrix with corresponding metadata, including donor metadata (sex and age group, tech_plus_study) and cell metadata (number of UMIs). The NEBULA-LN method was used with a mixed effect model treating donor, sex, age group, and tech_plus_study as fixed effects and log 10 (UMIs) as an offset term. Since the effect size is at the single-cell level, we considered DEGs to be those with abs(log 2 FC) > 0.1 and p adj. < 0.05. We did not continue with this approach, which seemed to be sensitive to ambient RNA and did not recover known sex chromosome genes as DEGs for the sex analysis (Additional file 1 : Fig. S7). Differential accessible region (DAR) analysis Similar to DEG analysis, we performed donor-based pseudobulked accessibility analysis by aggregating the raw counts of all nuclei of the same cell type annotation for each donor. We used the approach 1 outlined in the “Differential gene expression analysis” section. We performed differential accessibility analysis with pydeseq2 [ 44 ] using this design matrix: expression ~ sex + age_group + disease_binary + technology. Study was not included in the design matrix to reduce collinearity, as all the diseased cells come from the same study [ 128 ]. snRNA-seq trajectory analysis Trajectories from fetal to non-diseased were calculated by aggregating gene expression along the density of fetal cells and non-diseased in the embedding space for cardiomyocytes, endothelial cells, fibroblasts, and myeloid. We found that this density gradient approach performed better in terms of constructing a trajectory from non-diseased states to diseased states than other pseudotime-based methods such as Palantir and DPT, which work better for developmental trajectories. Additionally, trajectories from non-diseased to diseased were calculating by aggregating gene expression along the density from non-diseased cells to diseased cells in the embedding space, after filtering to include cells belonging to only two categories (e.g., non-diseased and DCM). Genes were hierarchically clustered based on expression dynamics. Clusters were assigned using the fcluster function from scipy using the “distance” threshold set to 25% of the maximum linkage distance in the linkage matrix. Over-representation analysis against the Hallmark pathways was performed for each cluster. Gene set enrichment analysis (GSEA) and over-representation analysis (ORA) GSEA was performed using gseapy (v1.0.3) [ 109 ], a fast python implementation of GSEA [ 110 ] against the Hallmark gene sets [ 111 ]. After performing the DESeq2 approach 1 for each contrast and cell type, we obtained a results table for each gene, including log 2 fold change. Genes were sorted in terms of decreasing log 2 fold change and then passed into gseapy prerank() function. Significant GSEA results were defined as gene sets with an FDR-adjusted p value < 0.05. For over-representation analysis, the list of fetal reactivation genes was inputted against the Hallmark gene sets. Fetal reactivation genes were defined as directionally concordant DEGs that are shared between the developmental and disease contrasts. For example, fetal reactivation “up” genes are those that are higher in the fetal hearts relative to postnatal young heart and higher in diseased hearts relative to postnatal non-diseased hearts. DEG/DAR concordant overlap significance across contrasts To determine whether there is a significant concordant overlap between DEGs in two contrasts within a given cell type, we simulated expected overlap under a null distribution that accounts for the size of the gene sets (Additional file 1 : Fig. S18a-b). For example, for two contrasts 1 (e.g., disease vs non-diseased) and 2 (e.g., fetal vs young), there are four total gene sets of varying size: up in 1 (W), down in 1 (X), up in 2 (Y), down in 2 (Z). W and X are disjoint (no overlapping elements) as no gene can both be up- or down-regulated in the same contrast. Y and Z are also disjoint. The amount of concordant overlap is the sum of the total number of genes that intersect between W and Y (referred to as “a”) plus the total number of genes that intersect between X and Z (referred to as “d”). Discordant overlap is the sum of the total number of genes that intersect between W and Z (referred to as “b”) and those that intersect between X and Y (referred to as “c”). To obtain a proportion of concordant overlap, we divide the total number of directionally concordant DEGs (“a” + “d”) by the total number of shared DEGs (“a” + “b” + “c” + “d”). To simulate a null distribution, we randomly sample W genes (“up in 1”) from the gene universe (“G”; all genes represented in the count matrix for that cell type) and then remove those W genes from “G” to sample X (“down in 1”) genes. We then sample Y genes (“up in 2”) from the gene universe and then remove those Y genes from “G” to sample Z genes (“down in 2”) genes. We calculate the degree of overlap across 10,000 simulations to obtain an approximately normal null distribution with estimated μ and σ. The Z-score of the observed degree of overlap [O = (a + d)/a + b + c + d)] is calculated as (O—μ)/σ. A positive Z-score would indicate concordant DEG overlap between two contrasts greater than expected by chance, while a negative Z-score would indicate overlap lower than expected by chance (mutual exclusivity). A Z-score greater than 3 was considered significant. Cell–cell communication analysis Liana x tensorcell2cell (v1.2.1) was performed on scVI-batched integrated counts, as recommended by the software developers. The counts were all transformed as if generated from the ENCODE v4 Multiome batch, using the scVI transform_batch() function. Without performing this batch correction, we obtained spurious significant differences in cell–cell communication that seemed to be driven by sequencing depth differences and other study-specific technical artifacts. We followed the tutorial titled “Intercellular Context Factorization with Tensor-Cell2cell” to build a 4D tensor, representing contexts, interactions, sender cell types, and receiver cell types. To reduce jargon, “program” was used in this manuscript instead of “factor”. The number of programs identified using optimal rank estimation was 7. Next, we identified programs that showed significant differences across three age + disease status conditions: fetal non-diseased, postnatal non-diseased, and postnatal diseased using the c2c.plotting . context_boxplot() function. Enriched pathways were identified using decoupleR [ 112 ] with PROGENy [ 113 ] gene sets. SCENIC + gene regulatory network analysis SCENIC + [ 114 ] (v1.0a1) integrates snATAC-seq and snRNA-seq to infer gene regulatory networks. As SCENIC + is not highly scalable beyond > 30 K cells, we subsampled nuclei for this analysis. For sex and aging analysis, we subsampled 3000 nuclei per cell type divided evenly across six age group + sex combinations (male-young, female-young, male-middle, female-middle, male-old, female-old) from non-diseased donors. We only performed analysis on five cell types that had enough cells for each combination: cardiomyocyte, endothelial, fibroblast, myeloid, and pericyte, resulting in 15,000 total cells. For fetal and disease analysis, we subsampled 1000 cells for each age_status + disease category (fetal, non-diseased postnatal, diseased postnatal) per cell type, using 15,000 cells as well. The rest of the analysis was performed using the SCENIC + tutorial ( https://scenicplus.readthedocs.io/en/latest/tutorials.html ). Using the area under the curve (AUC) metric that measures the activity of a TF-driven GRN per cell, we compared cells for each contrast (e.g., old vs. young) in each cell type. Significantly differential GRNs were defined as those with log 2 FC > 0.25 and adjusted p -value < 0.05 based on a Wilcoxon ranked sum test. For this analysis, we only examined GRNs for activators where both TF expression and chromatin accessibility of TF binding sites is positively correlated with gene expression (+/+ in SCENIC + notation). Repressor regulons were removed from analysis, as recommended by the SCENIC + authors, due to difficulty in interpretation. Additionally, we added two filtering steps to remove possible false positive regulons. First, as described here [ 114 ], we computed the correlation between meta-cell TF expression and meta-cell regulon activity in a cell-type–specific manner. Regulons with a correlation lower than 0.7 were removed. The most common regulons filtered out in this step were AP-1 TFs, which may reflect the redundancy of these TFs or their rapid expression due to dissociation-induced stress artifacts [ 115 ]. Secondly, we removed regulons that have a lower activity than 0.05 in the non-diseased meta-cells for a given cell type. This removes regulons that have low activity from further analysis. Spatial transcriptomics analysis Processed spatial transcriptomics datasets from Kuppe et al. 2022 [ 11 ] and Kanemaru et al. 2023 [ 16 ] were downloaded. Combined spot level expression of fetal reactivation genes were calculated using scanpy’s score_gene() function [ 101 ] using the gene set of shared upregulated DEGs in the developmental and disease contrasts. A weighted fetal reactivation score was computed by multiplying the per cell-type fetal reactivation score by the cell-type proportion, computed using cell2location (v0.1) [ 116 ]. The proportion of fetal reactivation spots for each sample was determined as the (# spots with fetal reactivation score > 0)/(# total spots). For spatial autocorrelation, Moran’s I was calculated for the fetal reactivation scores for each cell type. To determine which cell types are enriched for fetal reactivation spots, we performed multivariable linear regression analysis in which we model the fetal reactivation score s f ~ x 1 β 1 + x 2 β 2 + … x n β n , where s f corresponds to the fetal reactivation score within a given spot and where x i corresponds to the centered log-ratio transformed cell-type abundance within a given spot. This analysis was only performed for diseased sections with greater than 5% of spots with positive fetal reactivation scores. The cell-type proportions were transformed to remove the compositional constraint. Prioritization of target genes for GWAS loci We used scE2G (v1.2) [ 92 ] to perform peak-to-gene linkages in a cell type specific manner. Using 10X Multiome cells in our atlas (from Kanemaru 2023 and ENCODE v4, which were all non-diseased donors), we downsampled the number of cells to a maximum of 10,000 cells per cell type randomly selected from all donors. This threshold was chosen based on method runtime scalability, and since peak-to-gene linkage does not substantially improve with more cells. We ran the scE2G Multiome method for all cell types that had at least 1 million pseudobulked UMIs and 2 million ATAC fragments based on method’s recommendations. 9 cell types (excluding Adipocyte and Mast cells) met these criteria. Per scE2G recommendations, we defined high confidence enhancer-to-promoter linkages as having a score > 0.17, as this threshold was previously shown to yield 70% recall when evaluating K562 CRISPR inhibition-validated enhancer-gene pairs. GWAS lead SNPs were obtained from recent GWAS studies for dilated cardiomyopathy (DCM) [ 89 , 90 ] and hypertrophic cardiomyopathy (HCM) [ 91 ]. For the two DCM GWAS, the lead SNP sometimes differed between the studies. Therefore, we merged nearby lead SNPs within 50 kb that associated with the same prioritized gene into 117 distinct loci. Then, using PLINK (v1.90b6.21) (plink –r2 –ld-window-kb 1000 –ld-window 99,999), we identified all SNPs within 1 megabase that have an r 2 > 0.8, using the 1 K Genomes high coverage [ 117 ] EUR variant call files to construct linkage disequilibrium blocks. With this expanded set of SNPs, we used bedtools [ 105 ] to identify SNPs located in linked peaks called using scE2G. For these SNPs, we designated the gene linked to each enhancer as the prioritized gene. In cases where multiple enhancers per SNP map to different genes, we considered the prioritized genes as the top two genes by E2G score. This gene/set of genes was used to compute concordance between our approach and the prioritized genes from the original GWAS studies. “Full concordance” occurs if the set of scE2G and original study prioritized genes is an exact match. “Partial concordance” occurs when the prioritized gene list shared at least one gene, but do not fully match (e.g., we prioritized RRAS2/SPON1 for DCM locus chr11:14,353,533 but the original study prioritizes RRAS2/COPB1 ). “Discordance” occurs when no genes are shared between our scE2G-prioritized gene set and the original study gene set. We observed that most SNPs within strong LD (r 2 > 0.8) had similar overlap with scE2G peaks, so their overlap results were merged for visualization in Additional file 1 : Fig. S29 and Additional file 1 : Fig. S30. In cases, where there are SNPs corresponding to the same locus (lead SNP and expanded set of SNPs) that have distinct scE2G peak overlap profiles, these results are shown separately but included together in a dotted box. Visualization of genome browser tracks in Fig. 7 f was performed using the Integrated Genomics Viewer browser [ 118 ]. Integrative multiomic atlas web browser Our multiomic atlas is hosted by the developers at UCSC Cell Browser and available here: https://multiomic-human-heart.cells.ucsc.edu [ 119 ]. Besides the significant increase of cells, our resource includes the following features that are specific to the UCSC Cell Browser and not present in other existing human heart atlases. These features include (1) ability to color cells by gene expression or metadata, (2) default display of cell types by dot plots, violin plots, custom color palette, marker genes, dataset genes, (3) display of metadata histograms when cells are selected, (4) a split screen mode, which improves user experience. Additionally, for the ATAC-seq dataset, bigwig tracks for each cell type (divided into fetal, non-diseased, and diseased) are easily visualized on the UCSC Genome Browser. The peaks are annotated by distance between the peak and the nearest gene, increasing the user interpretability. Supplementary Information 13059_2026_4061_MOESM1_ESM.docx (79.6MB, docx) Additional file 1: Supplemental Figures (Figs. S1-S31) [ 133 – 174 ]. 13059_2026_4061_MOESM2_ESM.xlsx (15.3KB, xlsx) Additional file 2: Table S1. List of data sources for integration. 13059_2026_4061_MOESM3_ESM.xlsx (32.4KB, xlsx) Additional file 3: Table S2. Donor metadata for snRNA-seq datasets. 13059_2026_4061_MOESM4_ESM.xlsx (20.3KB, xlsx) Additional file 4: Table S3. Donor metadata for snATAC-seq datasets. 13059_2026_4061_MOESM5_ESM.xlsx (56.2MB, xlsx) Additional file 5: Table S4. Sex differentially expressed genesand differentially accessible regionsper cell type. 13059_2026_4061_MOESM6_ESM.xlsx (56.2MB, xlsx) Additional file 6: Table S5. Aging DEGs and DARs per cell type. 13059_2026_4061_MOESM7_ESM.xlsx (57.4MB, xlsx) Additional file 7: Table S6. Developmental DEGs and DARs per cell type. 13059_2026_4061_MOESM8_ESM.xlsx (54.1MB, xlsx) Additional file 8: Table S7. Binarized disease DEGs and DARs per cell type. 13059_2026_4061_MOESM9_ESM.xlsx (13MB, xlsx) Additional file 9: Table S8. DCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM10_ESM.xlsx (12.6MB, xlsx) Additional file 10: Table S9. HCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM11_ESM.xlsx (13.3MB, xlsx) Additional file 11: Table S10. HCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM12_ESM.xlsx (16.8KB, xlsx) Additional file 12: Table S11. Prioritized genes for DCM GWAS loci. 13059_2026_4061_MOESM13_ESM.xlsx (14.3KB, xlsx) Additional file 13: Table S12. Prioritized genes for DCM GWAS loci. Acknowledgements We thank members of the Wu Lab and the Leducq Network on “The Inflammatory-Fibrosis Network in Ischemic Heart Failure” for helpful discussions. Peer review information Clint Miller and Claudia Feng were the primary editors of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team. The peer-review history is available in the online version of this article. Authors’ contributions W.G and H.W. conceived of this project, in collaboration with K.M. and K.S. K.B. and K.M. procured all postnatal human heart tissue samples. K.S. procured all fetal human heart tissue samples. W.G. combined Penn datasets with external datasets and performed all integrative sn/scRNA-seq, snATAC-seq, and spatial transcriptomics analyses, with suggestions from H.Z. and Y.L. P.H. generated all new snRNA-seq and snATAC-seq datasets and performed initial Penn dataset analyses, with help from Q.Q and K.X. B.W. and M.H. created the UCSC Cell Browser. W.G. and H.W. wrote the manuscript, with contributions from all authors. All authors read and approved the final manuscript. Funding This work was supported by a National Heart Lung and Blood Institute (NHLBI) grant DP2- HL142044 (to H.W.). W.G. and H.W. were supported by a NHGRI grant U01-HG012047 and a Leducq Foundation grant 20CVD02. W.G. was supported by NIH Medical Scientist Training Program grants T32GM007170 & T32GM148377 and a T32 Training Grant in Computational Genomics T32HG000046-26. B.W. and M.H. were supported by NIMH RF1MH132662, CIRM DISC0-14514, and NHGRI U24HG002371. Data availability The processed aggregated snRNA-seq dataset with 2,275,105 nuclei, the aggregated snATAC-seq dataset with 690,044 nuclei, and the newly generated raw expression matrices and raw sequence files generated for this study (labeled as “This study” in Additional file 1: Fig. S1a and Additional file 1: Fig. S6a) will be available on the Gene Expression Omnibus (GEO) under the accession GSE290367 [ 120 ]. All code used for analysis in this study are provided in Github ( https://github.com/wgao688/Human-cardiac-multiome-analysis ) [ 121 ] and are deposited under Zenodo DOI ( https://zenodo.org/uploads/18061421 ) [ 122 ] under the MIT license. The previously published datasets included in this integrated analysis are included in Additional file 2: Table S1. For snRNA-seq, they include Chaffin 2022 [ 123 ], ENCODE v4 [ 124 ], Hill 2022 [ 125 ], Kanemaru 2023 [ 126 ], Koenig 2022 [ 127 ], Kuppe 2022 [ 128 ], Litvinukova 2020 [ 129 ], Reichart 2022 [ 130 ], Sim 2021 [ 131 ], and Simonson 2023 [ 132 ]. For snATAC-seq, they include ENCODE v4 [ 124 ], Kanemaru 2023 [ 126 ], and Kuppe 2022 [ 128 ]. Declarations Ethics approval and consent to participate For all new snRNA-seq and snATAC-seq datasets generated in this study, we have received approval through the Penn Heart Tissue Biobank, which has been approved by the Penn IRB (#848421) since 2005. For fetal donors, the IRB number is 832470. Written informed consent for research use of donated tissue was obtained from next of kin in all cases. The consent forms explicitly provide permission to generate and disseminate genomic data from de-identified subjects. The experimental methods used in this work comply with the Helsinki Declaration. Competing interests The authors declare no competing interests. Footnotes Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. William Gao and Peng Hu contributed equally to this work. Contributor Information William Gao, Email: [email protected]. Hao Wu, Email: [email protected]. References 1. Berry JD, Dyer A, Cai X, Garside DB, Ning H, Thomas A, et al. Lifetime risks of cardiovascular disease. N Engl J Med. 2012;366(26):321–9. 10.1056/NEJMoa1012848. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 2. Spurrell CH, Barozzi I, Kosicki M, Mannion BJ, Blow MJ, Fukuda-Yuzawa Y, et al. Genome-wide fetalization of enhancer architecture in heart disease. Cell Rep. 2022;40(20):111400. 10.1016/j.celrep.2022.111400. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 3. Cordero P, Parikh VN, Chin ET, Erbilgin A, Gloudemans MJ, Shang C, et al. Pathologic gene network rewiring implicates PPP1R3A as a central regulator in pressure overload heart failure. Nat Commun. 2019;10(24):2760. 10.1038/s41467-019-10591-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 4. GTEX CONSORTIUM. The GTEx consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369(11):1318–30. 10.1126/science.aaz1776. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 5. Miranda AMA, Janbandhu V, Maatz H, Kanemaru K, Cranley J, Teichmann SA, et al. Single-cell transcriptomics for the assessment of cardiac disease. Nat Rev Cardiol. 2023;20(5):289–308. 10.1038/s41569-022-00805-7. [ DOI ] [ PubMed ] [ Google Scholar ] 6. Hu P, Liu J, Zhao J, Wilkins BJ, Lupino K, Wu H, et al. Single-nucleus transcriptomic survey of cell diversity and functional maturation in postnatal mammalian hearts. Genes Dev. 2018;32(19–20):1344–57. Available from: https://genesdev.cshlp.org/content/32/19-20/1344.short . [ DOI ] [ PMC free article ] [ PubMed ] 7. Tucker NR, Chaffin M, Fleming SJ, Hall AW, Parsons VA, Bedi KC, et al. Transcriptional and cellular diversity of the human heart. Circulation. 2020;142(4):466–82. 10.1161/CIRCULATIONAHA.119.045401. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 8. Litviňuková M, Talavera-López C, Maatz H, Reichart D, Worth CL, Lindberg EL, et al. Cells of the adult human heart. Nature. 2020;588(7838):466–72. 10.1038/s41586-020-2797-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 9. Buenrostro JD, Wu B, Litzenburger UM, Ruff D, Gonzales ML, Snyder MP, et al. Single-cell chromatin accessibility reveals principles of regulatory variation. Nature. 2015;523(23):486–90. 10.1038/nature14590. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 10. Ameen M, Sundaram L, Shen M, Banerjee A, Kundu S, Nair S, et al. Integrative single-cell analysis of cardiogenesis identifies developmental trajectories and non-coding mutations in congenital heart disease. Cell. 2022;185(22):4937-4953.e23. 10.1016/j.cell.2022.11.028. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 11. Kuppe C, Ramirez Flores RO, Li Z, Hayat S, Levinson RT, Liao X, et al. Spatial multi-omic map of human myocardial infarction. Nature. 2022;608(7924):7924. 10.1038/s41586-022-05060-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 12. Nguyen HCT, Baik B, Yoon S, Park T, Nam D. Benchmarking integration of single-cell differential expression. Nat Commun. 2023;14(1):1. 10.1038/s41467-023-37126-3. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 13. Luecken MD, Büttner M, Chaichoompu K, Danese A, Interlandi M, Mueller MF, et al. Benchmarking atlas-level data integration in single-cell genomics. Nat Methods. 2022;19(1):41–50. 10.1038/s41592-021-01336-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 14. Zimmerman KD, Espeland MA, Langefeld CD. A practical solution to pseudoreplication bias in single-cell studies. Nat Commun. 2021;12(1):738. 10.1038/s41467-021-21038-1. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 15. Chaffin M, Papangeli I, Simonson B, Akkad AD, Hill MC, Arduini A, et al. Single-nucleus profiling of human dilated and hypertrophic cardiomyopathy. Nature. 2022;608(7921):174–80. 10.1038/s41586-022-04817-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Kanemaru K, Cranley J, Muraro D, Miranda AMA, Ho SY, Wilbrey-Clark A, et al. Spatially resolved multiomics of human cardiac niches. Nature. 2023;619(7971):801–10. 10.1038/s41586-023-06311-1. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 17. Sim CB, Phipson B, Ziemann M, Rafehi H, Mills RJ, Watt KI, et al. Sex-specific control of human heart maturation by the progesterone receptor. Circulation. 2021;143(16):1614–28. 10.1161/CIRCULATIONAHA.120.051921. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 18. Hill MC, Kadow ZA, Long H, Morikawa Y, Martin TJ, Birks EJ, et al. Integrated multi-omic characterization of congenital heart disease. Nature. 2022;608(7921):181–91. 10.1038/s41586-022-04989-3. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 19. Koenig AL, Shchukina I, Amrute J, Andhey PS, Zaitsev K, Lai L, et al. Single-cell transcriptomics reveals cell-type-specific diversification in human heart failure. Nat Cardiovasc Res. 2022;1(3):263–80. 10.1038/s44161-022-00028-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 20. Reichart D, Lindberg EL, Maatz H, Miranda AMA, Viveiros A, Shvetsov N, et al. Pathogenic variants damage cell composition and single cell transcription in cardiomyopathies. Science. 2022;377(6606):eabo1984. 10.1126/science.abo1984. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 21. Simonson B, Chaffin M, Hill MC, Atwa O, Guedira Y, Bhasin H, et al. Single-nucleus RNA sequencing in ischemic cardiomyopathy reveals common transcriptional profile underlying end-stage heart failure. Cell Rep. 2023;42(2):112086. 10.1016/j.celrep.2023.112086. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 22. Hitz BC, Jin-Wook L, Jolanki O, Kagda MS, Graham K, Sud P, et al. The ENCODE Uniform Analysis Pipelines. bioRxiv. 2023 Apr 6;2023.04.04.535623. 10.1101/2023.04.04.535623 PubMed PMID: 37066421; PubMed Central PMCID: PMC10104020. 23. Luo Y, Hitz BC, Gabdank I, Hilton JA, Kagda MS, Lam B, et al. New developments on the encyclopedia of DNA elements (ENCODE) data portal. Nucleic Acids Res. 2020;48(D1):D882-9. 10.1093/nar/gkz1062. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 24. ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;6(7414):57–74. 10.1038/nature11247. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 25. Shen X, Wang C, Zhou X, Zhou W, Hornburg D, Wu S, et al. Nonlinear dynamics of multi-omics profiles during human aging. Nat Aging. 2024;14:1–16. 10.1038/s43587-024-00692-2. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 26. Sikkema L, Ramírez-Suástegui C, Strobl DC, Gillett TE, Zappia L, Madissoon E, et al. An integrated cell atlas of the lung in health and disease. Nat Med. 2023;29(6):1563–77. 10.1038/s41591-023-02327-2. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 27. Lopez R, Regier J, Cole MB, Jordan MI, Yosef N. Deep generative modeling for single-cell transcriptomics. Nat Methods. 2018;15(12):1053–8. 10.1038/s41592-018-0229-2. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 28. McInnes L, Healy J, Melville J. UMAP: uniform manifold approximation and projection for dimension reduction. arXiv. 2020. 10.48550/arXiv.1802.03426. [ Google Scholar ] 29. Zafeiridis A, Jeevanandam V, Houser SR, Margulies KB. Regression of cellular hypertrophy after left ventricular assist device support. Circulation. 1998;18(7):656–62. 10.1161/01.cir.98.7.656. [ DOI ] [ PubMed ] [ Google Scholar ] 30. Lowry JL, Brovkovych V, Zhang Y, Skidgel RA. Endothelial nitric-oxide synthase activation generates an inducible nitric-oxide synthase-like output of nitric oxide in inflamed endothelium. J Biol Chem. 2013;8(6):4174–93. 10.1074/jbc.M112.436022. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 31. Frasca D, Diaz A, Romero M, Garcia D, Blomberg BB. B cell immunosenescence. Annu Rev Cell Dev Biol. 2020;6:551–74. 10.1146/annurev-cellbio-011620-034148. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 32. Bajpai G, Bredemeyer A, Li W, Zaitsev K, Koenig AL, Lokshina I, et al. Tissue resident CCR2− and CCR2+ cardiac macrophages differentially orchestrate monocyte recruitment and fate specification following myocardial injury. Circ Res. 2019;18(2):263–78. 10.1161/CIRCRESAHA.118.314028. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 33. Zhang K, Zemke NR, Armand EJ, Ren B. SnapATAC2: a fast, scalable and versatile tool for analysis of single-cell omics data. bioRxiv. 2023. 10.1101/2023.09.11.557221. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 34. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–96. 10.1038/s41592-019-0619-0. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 35. Weinand K, Langan EM, Curtis M, Raychaudhuri S. Defining effective strategies to integrate multi-sample single-nucleus ATAC-seq datasets via a multimodal-guided approach. bioRxiv. 2025. 10.1101/2025.04.02.646871. [ Google Scholar ] 36. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;17(9):R137. 10.1186/gb-2008-9-9-r137. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 37. Speir ML, Bhaduri A, Markov NS, Moreno P, Nowakowski TJ, Papatheodorou I, et al. UCSC cell browser: visualize your single-cell data. Bioinformatics. 2021;37(7):4578–80. 10.1093/bioinformatics/btab503. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 38. Squair JW, Gautier M, Kathe C, Anderson MA, James ND, Hutson TH, et al. Confronting false discoveries in single-cell differential expression. Nat Commun. 2021;12(28):5692. 10.1038/s41467-021-25960-2. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 39. Teo AYY, Squair JW, Courtine G, Skinnider MA. Best practices for differential accessibility analysis in single-cell epigenomics. Nat Commun. 2024;15:8805. 10.1038/s41467-024-53089-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 40. He L, Davila-Velderrain J, Sumida TS, Hafler DA, Kellis M, Kulminski AM. NEBULA is a fast negative binomial mixed model for differential or co-expression analysis of large-scale multi-subject single-cell data. Commun Biol. 2021;4(26):1–17. 10.1038/s42003-021-02146-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 41. Caglayan E, Ayhan F, Liu Y, Vollmer RM, Oh E, Sherwood CC, et al. Molecular features driving cellular complexity of human brain evolution. Nature. 2023;620(7972):7972. 10.1038/s41586-023-06338-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. Janssen P, Kliesmete Z, Vieth B, Adiconis X, Simmons S, Marshall J, et al. The effect of background noise and its removal on the analysis of single-cell expression data. Genome Biol. 2023;24:140. 10.1186/s13059-023-02978-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 43. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(5):550. 10.1186/s13059-014-0550-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 44. Muzellec B, Teleńczuk M, Cabeli V, Andreux M. PyDESeq2: a python package for bulk RNA-seq differential expression analysis. bioRxiv. 2022:2022.12.14.520412. 10.1101/2022.12.14.520412. [ DOI ] [ PMC free article ] [ PubMed ] 45. Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, et al. Improving reproducibility of differentially expressed genes in single-cell transcriptomic studies of neurodegenerative diseases through meta-analysis. Nat Commun. 2025;16:7436. 10.1038/s41467-025-62579-z. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 46. Ding J, Adiconis X, Simmons SK, Kowalczyk MS, Hession CC, Marjanovic ND, et al. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods. Nat Biotechnol. 2020;38(6):6. 10.1038/s41587-020-0465-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 47. Read DF, Booth GT, Daza RM, Jackson DL, Gladden RG, Srivatsan SR, et al. Single-cell analysis of chromatin and expression reveals age- and sex-associated alterations in the human heart. Commun Biol. 2024;7(26):1–14. 10.1038/s42003-024-06582-y. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 48. Lu S, Ma M, Mao X, Bacino CA, Jankovic J, Sutton VR, et al. De novo variants in FRMD5 are associated with developmental delay, intellectual disability, ataxia, and abnormalities of eye movement. Am J Hum Genet. 2022;109(10):1932–43. 10.1016/j.ajhg.2022.09.005. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 49. Guo T, Yin RX, Pan L, Yang S, Miao L, Huang F. Integrative variants, haplotypes and diplotypes of the CAPN3 and FRMD5 genes and several environmental exposures associate with serum lipid variables. Sci Rep. 2017;7(1):45119. 10.1038/srep45119. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 50. Phua AIH, Le TT, Tara SW, De Marvao A, Duan J, Toh DF, et al. Paradoxical higher myocardial wall stress and increased cardiac remodeling despite lower mass in females. J Am Heart Assoc. 2020;9(4):e014781. 10.1161/JAHA.119.014781. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 51. Enzan N, Miyazawa K, Koyama S, Kurosawa R, Ieki H, Yoshida H, et al. Genome-wide analysis of heart failure yields insights into disease heterogeneity and enables prognostic prediction in the Japanese population. Nat Commun. 2025;16(1):9680. 10.1038/s41467-025-64659-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 52. Rossi M, Banskota N, Shin CH, Anerillas C, Tsitsipatis D, Yang JH, et al. Increased PTCHD4 expression via m6A modification of PTCHD4 mRNA promotes senescent cell survival. Nucleic Acids Res. 2024;52(12):7261–78. 10.1093/nar/gkae322. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 53. Saul D, Kosinsky RL, Atkinson EJ, Doolittle ML, Zhang X, LeBrasseur NK, et al. A new gene set identifies senescent cells and predicts senescence-associated pathways across tissues. Nat Commun. 2022;13(1):4827. 10.1038/s41467-022-32552-1. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 54. Bartz J, Jung H, Wasiluk K, Zhang L, Dong X. Progress in discovering transcriptional noise in aging. Int J Mol Sci. 2023;24(4):4. 10.3390/ijms24043701. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 55. Vermeulen MC, Pearse R, Young-Pearse T, Mostafavi S. Mosaic loss of chromosome Y in aged human microglia. Genome Res. 2022;32(10):1795–807. 10.1101/gr.276409.121. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 56. Ibañez-Solé O, Ascensión AM, Araúzo-Bravo MJ, Izeta A. Lack of evidence for increased transcriptional noise in aged tissues. Elife. 2022;11:e80380. 10.7554/eLife.80380. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 57. Mastroeni D, Khdour OM, Delvaux E, Nolz J, Olsen G, Berchtold N, et al. Nuclear but not mitochondrial-encoded oxidative phosphorylation genes are altered in aging, mild cognitive impairment, and Alzheimer’s disease. Alzheimers Dement. 2017;13(5):510–9. 10.1016/j.jalz.2016.09.003. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 58. Rahman FA, Quadrilatero J. Altered senescence and mitochondrial transcriptome defines age-related changes in satellite cells. Front Cell Dev Biol. 2026. 10.3389/fcell.2025.1699206. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 59. Nielsen J, Christiansen J, Lykke-Andersen J, Johnsen AH, Wewer UM, Nielsen FC. A family of insulin-like growth factor II mRNA-binding proteins represses translation in late development. Mol Cell Biol. 1999;19(2):1262–70. 10.1128/MCB.19.2.1262. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 60. Chattergoon NN. Thyroid hormone signaling and consequences for cardiac development. J Endocrinol. 2019;242(1):T145-60. 10.1530/JOE-18-0704. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 61. Li X, Zhou J, Fang M, Yu B. Pregnancy immune tolerance at the maternal-fetal interface. Int Rev Immunol. 2020;39(6):247–63. 10.1080/08830185.2020.1777292. [ DOI ] [ PubMed ] [ Google Scholar ] 62. Cohen-Barak O, Yi Z, Hagiwara N, Monzen K, Komuro I, Brilliant MH. Sox6 regulation of cardiac myocyte development. Nucleic Acids Res. 2003;31(20):5941–8. 10.1093/nar/gkg807. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 63. Benhaddou A, Keime C, Ye T, Morlon A, Michel I, Jost B, et al. Transcription factor TEAD4 regulates expression of Myogenin and the unfolded protein response genes during C2C12 cell differentiation. Cell Death Differ. 2012;19(2):220–31. 10.1038/cdd.2011.87. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 64. Ali M, Liccardo D, Cao T, Tian Y. Natriuretic peptides and forkhead O transcription factors act in a cooperative manner to promote cardiomyocyte cell cycle re-entry in the postnatal mouse heart. BMC Dev Biol. 2021;21(1):6. 10.1186/s12861-020-00236-y. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 65. Song W, Wang H, Wu Q. Atrial natriuretic peptide in cardiovascular biology and disease (NPPA). Gene. 2015;569(1):1–6. 10.1016/j.gene.2015.06.029. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 66. Nakao K, Minobe W, Roden R, Bristow MR, Leinwand LA. Myosin heavy chain gene expression in human heart failure. J Clin Invest. 1997;100(9):2362–70. 10.1172/JCI119776. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 67. Vukusic K, Thorsell A, Muslimovic A, Jonsson M, Dellgren G, Lindahl A, et al. Overexpression of the SARS-CoV-2 receptor angiotensin converting enzyme 2 in cardiomyocytes of failing hearts. Sci Rep. 2022;12(1):965. 10.1038/s41598-022-04956-y. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 68. Beyerstedt S, Casaro EB, Rangel ÉB. COVID-19: angiotensin-converting enzyme 2 (ACE2) expression and tissue susceptibility to SARS-CoV-2 infection. Eur J Clin Microbiol Infect Dis. 2021;40(5):905–19. 10.1007/s10096-020-04138-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 69. Murphy SP, Kakkar R, McCarthy CP, Januzzi JL. Inflammation in heart failure. JACC. 2020;75(11):1324–40. 10.1016/j.jacc.2020.01.014. [ DOI ] [ PubMed ] [ Google Scholar ] 70. Bhullar SK, Dhalla NS. Status of mitochondrial oxidative phosphorylation during the development of heart failure. Antioxidants (Basel). 2023;12(11):1941. 10.3390/antiox12111941. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 71. Dirkx E, da Costa Martins PA, De Windt LJ. Regulation of fetal gene expression in heart failure. Biochimica et Biophysica Acta (BBA) - Molecular Basis of Disease. 2013;1832(12):2414–24. 10.1016/j.bbadis.2013.07.023. [ DOI ] [ PubMed ] [ Google Scholar ] 72. Rajabi M, Kassiotis C, Razeghi P, Taegtmeyer H. Return to the fetal gene program protects the stressed heart: a strong hypothesis. Heart Fail Rev. 2007;12(3):331–43. 10.1007/s10741-007-9034-1. [ DOI ] [ PubMed ] [ Google Scholar ] 73. van der Pol A, Hoes MF, de Boer RA, van der Meer P. Cardiac foetal reprogramming: a tool to exploit novel treatment targets for the failing heart. J Intern Med. 2020;288(5):491–506. 10.1111/joim.13094. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 74. Mehdiabadi NR, Boon Sim C, Phipson B, Kalathur RKR, Sun Y, Vivien CJ, et al. Defining the fetal gene program at single-cell resolution in pediatric dilated cardiomyopathy. Circulation. 2022;146(14):1105–8. 10.1161/CIRCULATIONAHA.121.057763. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 75. D’Antonio M, Nguyen JP, Arthur TD, Matsui H, Donovan MKR, D’Antonio-Chronowska A, et al. In heart failure reactivation of RNA-binding proteins is associated with the expression of 1,523 fetal-specific isoforms. PLoS Comput Biol. 2022;18(2):e1009918. 10.1371/journal.pcbi.1009918. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 76. Shi M, Yuan H, Li Y, Guo Z, Wei J. Targeting macrophage phenotype for treating heart failure: a new approach. Drug Des Devel Ther. 2024;18:4927–42. 10.2147/DDDT.S486816. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 77. Mohamud MA, İbrahim İG, Ahmed SA, Karataş M, Jeele MOO. Prevalence of thyroid dysfunction among patients with heart failure at a tertiary hospital in Mogadishu, Somalia. Int J Gen Med. 2022;15:6335–9. 10.2147/IJGM.S371697. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 78. Jiang X, Cao M, Wu J, Wang X, Zhang G, Yang C, et al. Protections of transcription factor BACH2 and natural product myricetin against pathological cardiac hypertrophy and dysfunction. Front Physiol. 2022;13:971424. 10.3389/fphys.2022.971424. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 79. Peisker F, Halder M, Nagai J, Ziegler S, Kaesler N, Hoeft K, et al. Mapping the cardiac vascular niche in heart failure. Nat Commun. 2022;13(1):3027. 10.1038/s41467-022-30682-0. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 80. Li H, Zou J, Yu XH, Ou X, Tang CK. Zinc finger E-box binding homeobox 1 and atherosclerosis: new insights and therapeutic potential. J Cell Physiol. 2021;236(6):4216–30. 10.1002/jcp.30177. [ DOI ] [ PubMed ] [ Google Scholar ] 81. Baghdassarian HM, Dimitrov D, Armingol E, Saez-Rodriguez J, Lewis NE. Combining lIANA and tensor-cell2cell to decipher cell-cell communication across multiple samples. Cell Rep Methods. 2024. 10.1016/j.crmeth.2024.100758. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 82. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes. Nat Protoc. 2020;15(4):1484–506. 10.1038/s41596-020-0292-x. [ DOI ] [ PubMed ] [ Google Scholar ] 83. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088. 10.1038/s41467-021-21246-9. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 84. Baruzzo G, Cesaro G, Di Camillo B. Identify, quantify and characterize cellular communication from single-cell RNA sequencing data with scSeqComm. Bioinformatics. 2022;38(7):1920–9. 10.1093/bioinformatics/btac036. [ DOI ] [ PubMed ] [ Google Scholar ] 85. Hou R, Denisenko E, Ong HT, Ramilowski JA, Forrest ARR. Predicting cell-to-cell communication networks using NATMI. Nat Commun. 2020;11(1):5011. 10.1038/s41467-020-18873-z. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 86. Raredon MSB, Yang J, Garritano J, Wang M, Kushnir D, Schupp JC, et al. Computation and visualization of cell–cell signaling topologies in single-cell systems data using Connectome. Sci Rep. 2022;12(1):4187. 10.1038/s41598-022-07959-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 87. Derrick CJ, Noël ES. The ECM as a driver of heart development and repair. Development. 2021;148(5):dev191320. 10.1242/dev.191320. [ DOI ] [ PubMed ] [ Google Scholar ] 88. Asp M, Salmén F, Ståhl PL, Vickovic S, Felldin U, Löfling M, et al. Spatial detection of fetal marker genes expressed at low level in adult human heart tissue. Sci Rep. 2017;7(1):12941. 10.1038/s41598-017-13462-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 89. Jurgens SJ, Rämö JT, Kramarenko DR, Wijdeveld LFJM, Haas J, Chaffin MD, et al. Genome-wide association study reveals mechanisms underlying dilated cardiomyopathy and myocardial resilience. Nat Genet. 2024. 10.1038/s41588-024-01975-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 90. Zheng SL, Jurgens SJ, McGurk KA, Xu X, Grace C, Theotokis PI, et al. Evaluation of polygenic scores for hypertrophic cardiomyopathy in the general population and across clinical settings. Nat Genet. 2025. 10.1038/s41588-025-02094-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 91. Tadros R, Zheng SL, Grace C, Jordà P, Francis C, West DM, et al. Large-scale genome-wide association analyses identify novel genetic loci and mechanisms in hypertrophic cardiomyopathy. Nat Genet. 2025. 10.1038/s41588-025-02087-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 92. Sheth MU, Qiu WL, Rosa Ma X, Gschwind AR, Jagoda E, Tan AS, et al. Mapping enhancer-gene regulatory interactions from single-cell data. bioRxiv. 2024:2024.11.23.624931. 10.1101/2024.11.23.624931. 93. Benton ML, Talipineni SC, Kostka D, Capra JA. Genome-wide enhancer annotations differ significantly in genomic distribution, evolution, and function. BMC Genomics. 2019;20(1):511. 10.1186/s12864-019-5779-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 94. Cattin ME, Wang J, Weldrick JJ, Roeske CL, Mak E, Thorn SL, et al. Deletion of MLIP (muscle-enriched A-type lamin-interacting protein) leads to cardiac hyperactivation of Akt/mammalian target of rapamycin (mTOR) and impaired cardiac adaptation. J Biol Chem. 2015;290(44):26699–714. 10.1074/jbc.M115.678433. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 95. Chen P, Li Z, Nie J, Wang H, Yu B, Wen Z, et al. MYH7B variants cause hypertrophic cardiomyopathy by activating the CaMK-signaling pathway. Sci China Life Sci. 2020;63(9):1347–62. 10.1007/s11427-019-1627-y. [ DOI ] [ PubMed ] [ Google Scholar ] 96. Peter AK, Rossi AC, Buvoli M, Ozeroff CD, Crocini C, Perry AR, et al. Expression of normally repressed Myosin Heavy Chain 7b in the mammalian heart induces dilated cardiomyopathy. J Am Heart Assoc. 2019;8(15):e013318. 10.1161/JAHA.119.013318. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 97. Banks D, Wu M, Higa LA, Gavrilova N, Quan J, Ye T, et al. L2DTL/CDT2 and PCNA interact with p53 and regulate p53 polyubiquitination and protein stability through MDM2 and CUL4A/DDB1 complexes. Cell Cycle. 2006;5(15):1719–29. 10.4161/cc.5.15.3150. [ DOI ] [ PubMed ] [ Google Scholar ] 98. Cuomo ASE, Nathan A, Raychaudhuri S, MacArthur DG, Powell JE. Single-cell genomics meets human genetics. Nat Rev Genet. 2023;24(8):535–49. 10.1038/s41576-023-00599-5. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 99. Hu P, Fabyanic E, Kwon DY, Tang S, Zhou Z, Wu H. Dissecting cell-type composition and activity-dependent transcriptional state in mammalian brains by massively parallel single-nucleus RNA-seq. Molecular cell [Internet]. 2017 [cited 2024 Jun 25];68(5):1006–15. Available from: https://www.cell.com/molecular-cell/pdf/S1097-2765(17)30876-6.pdf . [ DOI ] [ PMC free article ] [ PubMed ] 100. Sloan CA, Chan ET, Davidson JM, Malladi VS, Strattan JS, Hitz BC, et al. ENCODE data at the ENCODE portal. Nucleic Acids Res. 2016;44(D1):D726-732. 10.1093/nar/gkv1160. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 101. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19(1):15. 10.1186/s13059-017-1382-0. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 102. Young MD, Behjati S. SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing data. Gigascience. 2020;9(12):giaa151. 10.1093/gigascience/giaa151. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 103. Fleming SJ, Chaffin MD, Arduini A, Akkad AD, Banks E, Marioni JC, et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat Methods. 2023;20(9):1323–35. 10.1038/s41592-023-01943-7. [ DOI ] [ PubMed ] [ Google Scholar ] 104. Dijk D, Sharma R, Nainys J, Yim K, Kathail P, Carr AJ, et al. Recovering gene interactions from single-cell data using data diffusion. Cell. 2018;174(3):716-729.e27. 10.1016/j.cell.2018.05.061. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 105. Quinlan AR, Hall IM. Bedtools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2. 10.1093/bioinformatics/btq033. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 106. Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31(14):2382–3. 10.1093/bioinformatics/btv145. [ DOI ] [ PubMed ] [ Google Scholar ] 107. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. 10.1093/nar/gkv007. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 108. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc: Ser B (Methodol). 1995;57(1):289–300. 10.1111/j.2517-6161.1995.tb02031.x. [ Google Scholar ] 109. Fang Z, Liu X, Peltz G. GSEApy: a comprehensive package for performing gene set enrichment analysis in Python. Bioinformatics. 2023;39(1):btac757. 10.1093/bioinformatics/btac757. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 110. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50. 10.1073/pnas.0506580102. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 111. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25. 10.1016/j.cels.2015.12.004. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 112. Badia-i-Mompel P, Vélez Santiago J, Braunger J, Geiss C, Dimitrov D, Müller-Dott S, et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinformatics Adv. 2022;2(1):vbac016. 10.1093/bioadv/vbac016. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 113. Schubert M, Klinger B, Klünemann M, Sieber A, Uhlitz F, Sauer S, et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat Commun. 2018;9(1):20. 10.1038/s41467-017-02391-6. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 114. Bravo González-Blas C, De Winter S, Hulselmans G, Hecker N, Matetovici I, Christiaens V, et al. SCENIC+: single-cell multiomic inference of enhancers and gene regulatory networks. Nat Methods. 2023;20(9):9. 10.1038/s41592-023-01938-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 115. Marsh SE, Walker AJ, Kamath T, Dissing-Olesen L, Hammond TR, de Soysa TY, et al. Dissection of artifactual and confounding glial signatures by single-cell sequencing of mouse and human brain. Nat Neurosci. 2022;25(3):306–16. 10.1038/s41593-022-01022-8. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 116. Kleshchevnikov V, Shmatko A, Dann E, Aivazidis A, King HW, Li T, et al. Cell2location maps fine-grained cell types in spatial transcriptomics. Nat Biotechnol. 2022;40(5):661–71. 10.1038/s41587-021-01139-4. [ DOI ] [ PubMed ] [ Google Scholar ] 117. Byrska-Bishop M, Evani US, Zhao X, Basile AO, Abel HJ, Regier AA, et al. High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios. Cell. 2022;185(18):3426-3440.e19. 10.1016/j.cell.2022.08.004PubMedPMID:36055201;PubMedCentralPMCID:PMC9439720. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 118. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6. 10.1038/nbt.1754. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 119. Wick B, Haeussler M. Human multiomic atlas UCSC Trackhub. UCSC Genome Browser. 2026. https://genome.ucsc.edu/cgi-bin/hgTracks?hubUrl=https://cells.ucsc.edu/multiomic-human-heart/hub/hub.txt . 120. Gao W, Hu P, Qiu Q, Bedi K, Sasaki K, Margulies K, et al. Gao et al. 2026 human cardiac datasets. Gene Expression Omnibus. 2026. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE290367 . 121. Gao W. Source Code: An integrative single-nucleus multiomic atlas of the human left ventricle identifies gene regulatory network dynamics across development, aging, and disease. Github. 2026. https://github.com/wgao688/Human-cardiac-multiome-analysis/tree/main . [ DOI ] [ PMC free article ] [ PubMed ] 122. Gao W. Source code: An integrative single-nucleus multiomic atlas of the human left ventricle identifies gene regulatory network dynamics across development, aging, and disease. Zenodo. 2026. https://zenodo.org/records/18061421 . [ DOI ] [ PMC free article ] [ PubMed ] 123. Chaffin M. Single-nuclei profiling of human dilated and hypertrophic cardiomyopathy. Single Cell Portal; 2022. https://singlecell.broadinstitute.org/single_cell/study/SCP1303/single-nuclei-profiling-of-human-dilated-and-hypertrophic-cardiomyopathy#study-download [ DOI ] [ PMC free article ] [ PubMed ] 124. Snyder M. ENCODE v4 left ventricle snRNA-seq. ENCODE. 2025. https://www.encodeproject.org/search/ ? type=MultiomicsSeries&searchTerm=heart%20multiome. 125. Hill M. Integrated multiomic characterization of congenital heart disease. Gene ExpressionOmnibus. 2022. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE203275 . 126. Kanemaru K. Spatially resolved multiomics of human cardiac niches. Heart Cell Atlas. 2023. https://www.heartcellatlas.org . [ DOI ] [ PMC free article ] [ PubMed ] 127. Koenig AL. Cellular Atlas of Human Heart Failure. Gene Expression Omnibus. 2022. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE183852 . 128. Kuppe C. Human Myocardial Infarction atlas. Cellxgene. 2022. https://cellxgene.cziscience.com/collections/8191c283-0816-424b-9b61-c3e1d6258a77 . 129. Litvinukova M. Cells of the human heart. Heart Cell Atlas. 2022. https://www.heartcellatlas.org . 130. Reichart D. Pathogenic variants damage cell composition and single cell transcription in cardiomyopathies. Cellxgene. 2022. https://cellxgene.cziscience.com/collections/e75342a8-0f3b-4ec5-8ee1-245a23e0f7cb . [ DOI ] [ PMC free article ] [ PubMed ] 131. Porrello E, Sim C. Sex-specific control of human heart maturation by the progesterone receptor [snRNAseq_human]. Gene Expression Omnibus. 2021. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE156703 . [ DOI ] [ PMC free article ] [ PubMed ] 132. Simonson B. Single-nucleus RNA sequencing in ischemic cardiomyopathy reveals common transcriptional profile underlying end-stage heart failure. Single Cell Portal. 2023. https://singlecell.broadinstitute.org/single_cell/study/SCP1849/single-nucleus-rna-sequencing-in-ischemic-cardiomyopathy-reveals-common-transcriptional-profile-underlying-end-stage-heart-failure# /. [ DOI ] [ PMC free article ] [ PubMed ] 133. Stine RR, Shapira SN, Lim HW, Ishibashi J, Harms M, Won KJ, et al. EBF2 promotes the recruitment of beige adipocytes in white adipose tissue. Mol Metab. 2016;5(1):57–65. 10.1016/j.molmet.2015.11.001. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 134. Harms MJ, Ishibashi J, Wang W, Lim HW, Goyama S, Sato T, et al. Prdm16 is required for the maintenance of brown adipocyte identity and function in adult mice. Cell Metab. 2014;19(4):593–604. 10.1016/j.cmet.2014.03.007. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 135. Wang J, Zhu X, Wang S, Zhang Y, Hua W, Liu Z, et al. Phosphoproteomic and proteomic profiling in post-infarction chronic heart failure. Front Pharmacol. 2023;14:1181622. 10.3389/fphar.2023.1181622. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 136. Novel Truncated Peptide Derived From circCDYL Exacerbates Cardiac Hypertrophy | Circulation Research [Internet]. [cited 2026 Jan 25]. Available from: 10.1161/CIRCRESAHA.124.325573 [ DOI ] [ PubMed ] 137. van Berlo JH, Elrod JW, van den Hoogenhof MMG, York AJ, Aronow BJ, Duncan SA, et al. The transcription factor GATA-6 regulates pathological cardiac hypertrophy. Circ Res. 2010;107(8):10.1161/CIRCRESAHA.110.220764. 10.1161/CIRCRESAHA.110.220764. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 138. Lan C, Fang G, Li X, Chen X, Chen Y, Hu T, et al. SerpinB1 targeting safeguards against pathological cardiac hypertrophy and remodelling by suppressing cardiomyocyte pyroptosis and inflammation initiation. Cardiovasc Res. 2025;121(1):113–27. 10.1093/cvr/cvae241. [ DOI ] [ PubMed ] [ Google Scholar ] 139. Lu D, Zhang L, Bao D, Lu Y, Zhang X, Liu N, et al. Calponin1 inhibits dilated cardiomyopathy development in mice through the εPKC pathway. Int J Cardiol. 2014;173(2):146–53. 10.1016/j.ijcard.2014.02.032. [ DOI ] [ PubMed ] [ Google Scholar ] 140. Duarte A, Hirashima M, Benedito R, Trindade A, Diniz P, Bekman E, et al. Dosage-sensitive requirement for mouse Dll4 in artery development. Genes Dev. 2004;18(20):2474–8. 10.1101/gad.1239004. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 141. Kitagawa M, Hojo M, Imayoshi I, Goto M, Ando M, Ohtsuka T, et al. Hes1 and Hes5 regulate vascular remodeling and arterial specification of endothelial cells in brain vascular development. Mech Dev. 2013;130(9–10):458–66. 10.1016/j.mod.2013.07.001. [ DOI ] [ PubMed ] [ Google Scholar ] 142. Travisano SI, Harrison MRM, Thornton ME, Grubbs BH, Quertermous T, Lien CL. Single-nuclei multiomic analyses identify human cardiac lymphatic endothelial cells associated with coronary arteries in the epicardium. Cell Rep. 2023;42(9):113106. 10.1016/j.celrep.2023.113106. [ DOI ] [ PubMed ] [ Google Scholar ] 143. Luxán G, D’Amato G, MacGrogan D, de la Pompa JL. Endocardial notch signaling in cardiac development and disease. Circ Res. 2016;118(8):e1-18. 10.1161/CIRCRESAHA.115.305350. [ DOI ] [ PubMed ] [ Google Scholar ] 144. Carmon KS, Lin Q, Gong X, Thomas A, Liu Q. LGR5 interacts and cointernalizes with Wnt receptors to modulate Wnt/β-catenin signaling. Mol Cell Biol. 2012;32(11):2054–64. 10.1128/MCB.00272-12. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 145. Meyer D, Birchmeier C. Multiple essential functions of neuregulin in development. Nature. 1995;378(6555):386–90. 10.1038/378386a0. [ DOI ] [ PubMed ] [ Google Scholar ] 146. Settle S, Marker P, Gurley K, Sinha A, Thacker A, Wang Y, et al. The BMP family member Gdf7 is required for seminal vesicle growth, branching morphogenesis, and cytodifferentiation. Dev Biol. 2001;234(1):138–50. 10.1006/dbio.2001.0244. [ DOI ] [ PubMed ] [ Google Scholar ] 147. Maxwell MA, Muscat GEO. The NR4A subgroup: immediate early response genes with pleiotropic physiological roles. Nucl Recept Signal. 2006;4(8):e002. 10.1621/nrs.04002. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 148. Galea E, Feinstein DL. Regulation of the expression of the inflammatory nitric oxide synthase (NOS2) by cyclic AMP. FASEB J. 1999;13(15):2125–37. 10.1096/fasebj.13.15.2125. [ DOI ] [ PubMed ] [ Google Scholar ] 149. Chehboun S, Labrecque-Carbonneau J, Pasquin S, Meliani Y, Meddah B, Ferlin W, et al. Epstein-barr virus-induced gene 3 (EBI3) can mediate IL-6 trans -signaling. J Biol Chem. 2017;292(21):6644–56. 10.1074/jbc.M116.762021. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 150. Zhang M, Ming S, Gong S, Liang S, Luo Y, Liang Z, et al. Activation-induced cell death of mucosal-associated invariant T cells is amplified by OX40 in type 2 diabetic patients. J Immunol. 2019;203(15):2614–20. 10.4049/jimmunol.1900367. [ DOI ] [ PubMed ] [ Google Scholar ] 151. Schupp JC, Adams TS, Cosme C, Raredon MSB, Yuan Y, Omote N, et al. Integrated single-cell atlas of endothelial cells of the human lung. Circulation. 2021;144(27):286–302. 10.1161/CIRCULATIONAHA.120.052318. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 152. Chandrasekaran P, Negretti NM, Sivakumar A, Liberti DC, Wen H, Peers de Nieuwburgh M, et al. Cxcl12 defines lung endothelial heterogeneity and promotes distal vascular growth. Development. 2022;149(21):dev200909. 10.1242/dev.200909. (PubMed PMID: 36239312; PubMed Central PMCID: PMC9687018). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 153. Kong CSL, M V, Pantaleón-García J, Evans SE, Chen J. Truncated NTRK2 is induced in CAP1 endothelial cells during mouse lung injury-repair. iScience. 2025;28(7):112973. 10.1016/j.isci.2025.112973. (PubMed PMID: 40687816; PubMed Central PMCID: PMC12275058). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 154. Aghajanian H, Kimura T, Rurik JG, Hancock AS, Leibowitz MS, Li L, et al. Targeting cardiac fibrosis with engineered T cells. Nature. 2019;573(7774):430–3. 10.1038/s41586-019-1546-z. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 155. Yuan Q, Tang B, Zhu Y, Wan C, Xie Y, Xie Y, et al. PRDM16 acts as a therapeutic downstream target of TGF-β signaling in chronic kidney disease. JCI Insight. 2025. 10.1172/jci.insight.191458. [ DOI ] [ PMC free article ] [ PubMed ] 156. Patel VB, Zhong JC, Grant MB, Oudit GY. Role of the ACE2/angiotensin 1–7 axis of the renin–angiotensin system in heart failure. Circ Res. 2016;118(8):1313–26. 10.1161/CIRCRESAHA.116.307708. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 157. Sweetwyne MT, Murphy-Ullrich JE. Thrombospondin1 in tissue repair and fibrosis: TGF-β-dependent and independent mechanisms. Matrix Biol. 2012;31(3):178–86. 10.1016/j.matbio.2012.01.006. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 158. Jasso GJ, Jaiswal A, Varma M, Laszewski T, Grauel A, Omar A, et al. Colon stroma mediates an inflammation-driven fibroblastic response controlling matrix remodeling and healing. PLoS Biol. 2022;20(1):e3001532. 10.1371/journal.pbio.3001532. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 159. Bouvrée K, Brunet I, del Toro R, Gordon E, Prahst C, Cristofaro B, et al. Semaphorin3A, Neuropilin-1, and PlexinA1 are required for lymphatic valve formation. Circ Res. 2012. 10.1161/CIRCRESAHA.112.269316. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 160. Sichien D, Scott CL, Martens L, Vanderkerken M, Van Gassen S, Plantinga M, et al. IRF8 transcription factor controls survival and function of terminally differentiated conventional and plasmacytoid dendritic cells, respectively. Immunity. 2016;45(3):626–40. 10.1016/j.immuni.2016.08.013. [ DOI ] [ PubMed ] [ Google Scholar ] 161. Mac Donald A, Guipouy D, Lemieux W, Harvey M, Bordeleau LJ, Guay D, et al. KLRC1 knockout overcomes HLA-E-mediated inhibition and improves NK cell antitumor activity against solid tumors. Front Immunol. 2023;14:1231916. 10.3389/fimmu.2023.1231916. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 162. Kang JB, Raveane A, Nathan A, Soranzo N, Raychaudhuri S. Methods and insights from single-cell expression quantitative trait loci. Annu Rev Genomics Hum Genet. 2023;24(1):277–303. 10.1146/annurev-genom-101422-100437. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 163. Barkal AA, Weiskopf K, Kao KS, Gordon SR, Rosental B, Yiu YY, et al. Engagement of MHC class I by the inhibitory receptor LILRB1 suppresses macrophages and is a target of cancer immunotherapy. Nat Immunol. 2018;19(1):76–84. 10.1038/s41590-017-0004-z. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 164. Bajpai G, Schneider C, Wong N, Bredemeyer A, Hulsmans M, Nahrendorf M, et al. The human heart contains distinct macrophage subsets with divergent origins and functions. Nat Med. 2018;24(8):1234–45. 10.1038/s41591-018-0059-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 165. Rantakari P, Patten DA, Valtonen J, Karikoski M, Gerke H, Dawes H, et al. Stabilin-1 expression defines a subset of macrophages that mediate tissue homeostasis and prevent fibrosis in chronic liver injury. Proc Natl Acad Sci U S A. 2016;113(33):9298–303. 10.1073/pnas.1604780113. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 166. Hu Y, Bringmann H. Tfap2b acts in GABAergic neurons to control sleep in mice. Sci Rep. 2023;13(1):8026. 10.1038/s41598-023-34772-x. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 167. Li XJ, Morgan C, Duy PQ, Evsen L, Hao LT, Machavoine R, et al. The RNA-binding protein TRIM71 is essential for hearing in humans and mice and times auditory sensory organ development. Proc Natl Acad Sci U S A. 2025;122(36):e2505811122. 10.1073/pnas.2505811122. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 168. Okamoto R, Kato T, Mizoguchi A, Takahashi N, Nakakuki T, Mizutani H, et al. Characterization and function of MYPT2, a target subunit of myosin phosphatase in heart. Cell Signal. 2006;18(9):1408–16. 10.1016/j.cellsig.2005.11.001. [ DOI ] [ PubMed ] [ Google Scholar ] 169. Lee KW, Yeo SY, Gong JR, Koo OJ, Sohn I, Lee WY, et al. PRRX1 is a master transcription factor of stromal fibroblasts for myofibroblastic lineage progression. Nat Commun. 2022;13(1):2793. 10.1038/s41467-022-30484-4. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 170. Jahan F, Landry NM, Rattan SG, Dixon IMC, Wigle JT. The functional role of zinc finger E box-binding homeobox 2 (Zeb2) in promoting cardiac fibroblast activation. Int J Mol Sci. 2018;19(10):3207. 10.3390/ijms19103207. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 171. Zhao D, Han X, Huang L, Wang J, Zhang X, Jeon JH, et al. Transcription factor ZFHX3 regulates calcium influx in mammary epithelial cells in part via the TRPV6 calcium channel. Biochem Biophys Res Commun. 2019;519(2):366–71. 10.1016/j.bbrc.2019.08.148. [ DOI ] [ PubMed ] [ Google Scholar ] 172. Zhang J, Liu F, He Y, Zhang W, Ma W, Xing J, et al. Polycystin-1 downregulation induced vascular smooth muscle cells phenotypic alteration and extracellular matrix remodeling in Thoracic aortic dissection. Front Physiol. 2020. 10.3389/fphys.2020.548055. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 173. Read DF, Booth GT, Daza RM, Jackson DL, Gladden RG, Srivatsan SR, et al. Single-cell analysis of chromatin and expression reveals age- and sex-associated alterations in the human heart. bioRxiv; 2022. p. 2022.07.12.496461. Available from: 10.1101/2022.07.12.496461v210.1101/2022.07.12.496461. Cited 2024 Mar 11. [ DOI ] [ PMC free article ] [ PubMed ] 174. Zheng SL, Henry A, Cannie D, Lee M, Miller D, McGurk KA, et al. Genome-wide association analysis provides insights into the molecular etiology of dilated cardiomyopathy. Nat Genet. 2024. 10.1038/s41588-024-01952-y. [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Supplementary Materials 13059_2026_4061_MOESM1_ESM.docx (79.6MB, docx) Additional file 1: Supplemental Figures (Figs. S1-S31) [ 133 – 174 ]. 13059_2026_4061_MOESM2_ESM.xlsx (15.3KB, xlsx) Additional file 2: Table S1. List of data sources for integration. 13059_2026_4061_MOESM3_ESM.xlsx (32.4KB, xlsx) Additional file 3: Table S2. Donor metadata for snRNA-seq datasets. 13059_2026_4061_MOESM4_ESM.xlsx (20.3KB, xlsx) Additional file 4: Table S3. Donor metadata for snATAC-seq datasets. 13059_2026_4061_MOESM5_ESM.xlsx (56.2MB, xlsx) Additional file 5: Table S4. Sex differentially expressed genesand differentially accessible regionsper cell type. 13059_2026_4061_MOESM6_ESM.xlsx (56.2MB, xlsx) Additional file 6: Table S5. Aging DEGs and DARs per cell type. 13059_2026_4061_MOESM7_ESM.xlsx (57.4MB, xlsx) Additional file 7: Table S6. Developmental DEGs and DARs per cell type. 13059_2026_4061_MOESM8_ESM.xlsx (54.1MB, xlsx) Additional file 8: Table S7. Binarized disease DEGs and DARs per cell type. 13059_2026_4061_MOESM9_ESM.xlsx (13MB, xlsx) Additional file 9: Table S8. DCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM10_ESM.xlsx (12.6MB, xlsx) Additional file 10: Table S9. HCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM11_ESM.xlsx (13.3MB, xlsx) Additional file 11: Table S10. HCM vs. non-diseased DEGs and DARs per cell type. 13059_2026_4061_MOESM12_ESM.xlsx (16.8KB, xlsx) Additional file 12: Table S11. Prioritized genes for DCM GWAS loci. 13059_2026_4061_MOESM13_ESM.xlsx (14.3KB, xlsx) Additional file 13: Table S12. Prioritized genes for DCM GWAS loci. Data Availability Statement The processed aggregated snRNA-seq dataset with 2,275,105 nuclei, the aggregated snATAC-seq dataset with 690,044 nuclei, and the newly generated raw expression matrices and raw sequence files generated for this study (labeled as “This study” in Additional file 1: Fig. S1a and Additional file 1: Fig. S6a) will be available on the Gene Expression Omnibus (GEO) under the accession GSE290367 [ 120 ]. All code used for analysis in this study are provided in Github ( https://github.com/wgao688/Human-cardiac-multiome-analysis ) [ 121 ] and are deposited under Zenodo DOI ( https://zenodo.org/uploads/18061421 ) [ 122 ] under the MIT license. The previously published datasets included in this integrated analysis are included in Additional file 2: Table S1. For snRNA-seq, they include Chaffin 2022 [ 123 ], ENCODE v4 [ 124 ], Hill 2022 [ 125 ], Kanemaru 2023 [ 126 ], Koenig 2022 [ 127 ], Kuppe 2022 [ 128 ], Litvinukova 2020 [ 129 ], Reichart 2022 [ 130 ], Sim 2021 [ 131 ], and Simonson 2023 [ 132 ]. For snATAC-seq, they include ENCODE v4 [ 124 ], Kanemaru 2023 [ 126 ], and Kuppe 2022 [ 128 ]. Articles from Genome Biology are provided here courtesy of BMC ACTIONS View on publisher site PDF (6.3 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top