Abstract Recovery of cardiac function is the holy grail of heart failure therapy yet is infrequently observed and remains poorly understood. In this study, we performed single-nucleus RNA sequencing from patients with heart failure who recovered left ventricular systolic function after left ventricular assist device implantation, patients who did not recover and non-diseased donors. We identified cell-specific transcriptional signatures of recovery, most prominently in macrophages and fibroblasts. Within these cell types, inflammatory signatures were negative predictors of recovery, and downregulation of RUNX1 was associated with recovery. In silico perturbation of RUNX1 in macrophages and fibroblasts recapitulated the transcriptional state of recovery. Cardiac recovery mediated by BET inhibition in mice led to decreased macrophage and fibroblast Runx1 expression and diminished chromatin accessibility within a Runx1 intronic peak and acquisition of human recovery signatures. These findings suggest that cardiac recovery is a unique biological state and identify RUNX1 as a possible therapeutic target to facilitate cardiac recovery. __________________________________________________________________ Heart failure (HF) is a growing epidemic that affects over 23 million individuals worldwide and imparts a substantial burden on healthcare systems^[63]1,[64]2. Restoration of cardiac function, a phenomenon referred to as cardiac recovery, represents the holy grail for HF therapies. Cardiac recovery is characterized by restoration of left ventricular (LV) systolic function and reversal of LV dilatation, a process termed LV reverse remodeling. Cardiac recovery is infrequently observed after initiation of goal-directed medical therapy in ambulatory patients and implantation of left ventricular assist devices (LVADs) in patients with advanced HF^[65]3–[66]6. Individuals who experience cardiac recovery have markedly improved survival and quality of life^[67]7,[68]8. At present, much remains to be learned about cardiac recovery, including the underlying cellular and molecular mechanisms, features that distinguish patients with HF who will experience recovery from those who will continue to experience disease progression and whether cardiac recovery represents a reversion to normal or a unique biological state^[69]9. Unraveling the molecular basis of cardiac recovery in humans may help identify new therapeutic targets for HF. The advent of next-generation sequencing has afforded researchers new opportunities to dissect the cellular diversity of human tissue across disease contexts^[70]10,[71]11. Recent studies have used single-nucleus RNA sequencing (snRNA-seq) from healthy and failing human hearts to greatly expand our understanding of human HF and have identified cell-specific transcriptional signatures in genetic, dilated and hypertrophic cardiomyopathies^[72]12–[73]17. Although some studies have explored mechanical unloading of the heart, there is a paucity of information regarding cardiac recovery at the single-cell level. In particular, two studies used gene expression profiling to explore the bulk molecular signature of mechanical unloading after LVAD implantation^[74]18,[75]19. These studies revealed that the transcriptional profile of the unloaded heart differs from the non-failing donor and HF. However, neither study included samples from patients who recovered. More recently, Drakos et al.^[76]20 performed bulk RNA sequencing and proteomics in non-failing hearts and 92 patients with HF at the time of LVAD implantation (25 recovered after a period of mechanical unloading) to characterize features predictive of myocardial recovery. However, this study does not characterize the hearts once they recover. Designing and executing a study focused on cardiac recovery is challenging and must be carefully performed to discriminate cardiac recovery from mechanical unloading. To rigorously identify individuals who exhibit cardiac recovery, patients must undergo extensive phenotyping before and after LVAD implantation, and paired tissue specimens from the same patient need to be collected from similar regions in the LV at the time of LVAD implantation and at the time of LVAD explant. In the present study, we defined the cellular and transcriptional landscape of cardiac recovery using snRNA-seq. We compared non-diseased donor controls to patients with HF who experienced either cardiac recovery or ongoing HF before and after LVAD implantation (n = 40). By performing cell-specific pseudobulk differential expression analysis at the patient level, we uncovered that cardiac recovery represents a unique biological state, and we defined cell-specific signatures of cardiac recovery that were distinct from both healthy and failing hearts. We further showed that transcriptional changes of cardiac recovery were most prominent within macrophages and fibroblasts. We identified a RUNX1 gene regulatory network (GRN) in macrophages and fibroblasts that was associated with, and predictive of, cardiac recovery using a deep neural network. In silico perturbation analysis^[77]21 revealed that RUNX1 activity in macrophages and fibroblasts was predicted to control transition of cell states toward those associated with recovery and away from those associated with progressive HF. Finally, we leveraged a mouse model of cardiac recovery^[78]22 to test the role of the Runx1 GRN in cardiac recovery. These findings uncover the biological state and generate a hypothesis for cell-specific mechanisms that contribute to cardiac recovery, and they identify RUNX1 as a potential therapeutic target to facilitate cardiac recovery. Results snRNA-seq defines the cellular landscape of cardiac recovery We performed snRNA-seq on paired transmural LV specimens obtained from the apical anterior wall of age-matched and sex-matched donor controls (n = 14) and patients with HF at the time of LVAD implant and at the time of LVAD explant. Patients with HF were divided into those who recovered LV systolic function (recovery/reverse remodeled (RR), n = 5) and those who demonstrated persistent reduction in LV ejection fraction (mechanically unloaded (U), n = 8). Echocardiograms after LVAD implant were performed using a predefined protocol where the LVAD speed was reduced to assess intrinsic LV systolic function^[79]23,[80]24 ([81]Fig. 1a). After doublet removal and quality control (QC) ([82]Extended Data Fig. 1), we recovered 185,881 nuclei across 40 patient samples. We then performed dimensional reduction, uniform manifold approximation and projection (UMAP) construction and cell clustering with differential gene expression to annotate cell types ([83]Fig. 1b). We identified 13 transcriptionally distinct cell types marked by canonical gene markers ([84]Fig. 1c, [85]Extended Data Fig. 2a and [86]Supplementary Table 4). We also constructed cell-type-specific gene set scores and detected strong separation across clusters ([87]Extended Data Fig. 2b,[88]c). Cell composition analysis showed a drop in cardiomyocytes and an increase in the stromal cell fraction in HF pre-LVAD and post-LVAD samples relative to controls ([89]Fig. 1d and [90]Extended Data Fig. 2d). Echocardiographic data revealed the maked difference in LV ejection fraction after LVAD implant between the recovered and unloaded groups ([91]Fig. 1e). We then performed cell-type-specific pseudobulk differential gene expression analysis using the following comparisons: pre-LVAD HF (U-pre and RR-pre) versus donor, recovery post-LVAD (RR-post) versus pre-LVAD HF and RR-post versus donor ([92]Extended Data Fig. 3a,[93]b). Fig. 1 |. Study design, global clustering and DE analysis of cardiac remodeling after LVAD implantation. Fig. 1 | [94]Open in a new tab a, Study design. b, Integrated UMAP embedding plot of snRNA-seq data across n = 40 samples and 185,881 nuclei. c, Violin plot for canonical marker genes for cell types. d, Cell cluster composition across conditions colored by cell type. e, Paired ejection fraction measured before and after LVAD implantation split by RR (left) and U (right) with n = 5 and 8 biologically independent samples. Paired two-tailed t-test where *P = 0.0134 (RR) and non-significant P = 0.1169 (U). f, Pseudobulk differential gene expression comparisons in cardiac remodeling after LVAD implantation categories. g, Number of statistically significant differentially expressed genes (adjusted P < 0.05 and log[2] fold change > 0.58 from DEseq2) from f in each cell type as pairwise comparisons where size of dot refers to the sum of axis, and the color refers to the cell type. P values were calculated using the Wald test adjusted for multiple corrections. EF, ejection fraction; NK, natural killer; NS, non-significant. To place this differential expression (DE) analysis in the context of the different facets of cardiac recovery, we built a Venn diagram of overlapping differentially expressed genes ([95]Supplementary Tables 5–[96]7) to identify three key categories relevant to the response to HF and cardiac recovery: (1) genes expressed specifically during recovery (recovery genes); (2) genes associated with HF that were persistently changed (HF persistent genes); and (3) genes associated with HF that returned to normal (HF reversed genes) ([97]Fig. 1f). To identify cell types that may drive cardiac remodeling, we created pairwise scatter plots of the number of differentially expressed genes from the three comparisons ([98]Fig. 1g). Within the majority of cell types, HF persistent and HF recovery genes were found at higher frequency than HF reversed genes, indicating that the recovered state is not simply a reversion to normal but also a unique biological entity. Cell-type-specific signatures of cardiac recovery To dissect which cell types display the strongest signature of recovery, we quantified the number of recovery genes upregulated and downregulated in each cell type. Among the major cell populations, myeloid cells (monocytes, macrophages and dendritic cells) and fibroblasts had the greatest number of recovery genes ([99]Fig. 2a,[100]b and [101]Extended Data Fig. 3c). We then constructed heat maps of the top 25 upregulated and downregulated recovery genes in the cell types with the strongest recovery signatures (myeloid, fibroblast, endothelium, endocardium and cardiomyocyte) and observed consistent and robust enrichment within the RR-post group. We observed modest enrichment of the recovery signature within the U-post, suggesting that mechanical unloading may produce an intermediate phenotype within the continuum of recovery ([102]Fig. 2c). To validate our recovery signature in an independent dataset, we leveraged bulk RNA sequencing from Drakos et al.^[103]20. We trained a random forest (RF) classifier to predict recovery from bulk RNA sequencing profiles taken at the time of LVAD implantation and ranked feature importance based on mean decrease in impurity. Using the top-ranked genes, we calculated the number of overlapping genes with our cell-specific recovery signatures and found that myeloid cells and fibroblasts had the strongest overlap ([104]Fig. 2d). Fig. 2 |. Cell-type specific cardiac recovery. Fig. 2 | [105]Open in a new tab a, Schematic of gene set from DE analysis that marks cardiac recovery. b, Number of statistically significant (adjusted P < 0.05 and log[2] fold change > 0.58 from DESeq2) genes from pseudobulk DE analysis that are upregulated and downregulated in cardiac recovery across cell types. P values were calculated using the Wald test adjusted for multiple corrections. c, Pseudobulk heat maps of top genes upregulated (top) and downregulated (bottom) in cardiac recovery in major cell populations split by donor, HF pre-LVAD, reverse remodeled and unloaded. d, Overlapping genes between recovery predicted genes from bulk RNA sequencing in RR-pre and U-pre^[106]20 and pseudobulk cell-specific recovery genes. e, Number of unique and overlapping cardiac recovery genes in major cell populations. f, Pseudobulk expression of cardiac recovery genes that overlap among the major cell populations. g, Polygenic recovery score of upregulated and downregulated genes in cardiac recovery versus patient ejection fracrtion as a simple linear regression. Dotted line indicates 95% confidence interval; R^2 indicates goodness of fit; and P value indicates whether the slope is significantly non-zero using an F-test. P values were calculated using two-tailed linear regression Wald test with t-distribution. EF, ejection fraction; NK, natural killer; SMC, smooth muscle cell. To assess whether the transcriptional signature of recovery is unique to each cell type, we quantified the number of recovery genes that overlap in more than one cell type across the five major populations and found that most of the identified recovery genes are specific to a given cell type ([107]Fig. 2e). UpSet plots were constructed to assess pairwise overlap of cardiac recovery genes among different permutations of the major cell types ([108]Extended Data Fig. 4). Although the dominant signatures were cell type specific, we found five upregulated (PTPN13, FKBP5, ZBTB16, FGD4 and ZMYND8) and one downregulated (CTD-3252C9.4) recovery genes common to the five major cell types ([109]Fig. 2f). To ascertain whether cell-type-specific cardiac recovery signatures are associated with LV systolic function after LVAD implant, we performed a linear regression analysis using LV ejection fraction as a surrogate for systolic function. We identified strong correlations between the ejection fraction and upregulated and downregulated recovery signatures ([110]Fig. 2g). These findings indicate that transcriptional signatures of cardiac recovery are cell type specific and are reflective of a continuum of cardiac recovery. Polygenic predictors for cardiac recovery can be derived for any of the major cell types ([111]Fig. 2f), which highlights that each of the main cell types enters a distinct state across the continuum of recovery. Cardiomyocytes do not revert to a healthy state in cardiac recovery To dissect the recovery landscape across cell populations, we performed a pseudobulk principal component analysis (PCA) at the sample level ([112]Extended Data Fig. [113]5). We observed a strong separation between donor and HF samples across most cell types. Intriguingly, for cardiomyocytes, post-U and post-RR cardiomyocytes clustered with pre-LVAD HF cardiomyocytes, suggesting a persistence of the HF phenotype during both recovery and mechanical unloading ([114]Fig. 3a). To delineate whether recovery is associated with particular cardiomyocyte states, we performed high-resolution clustering of cardiomyocytes. We detected five distinct cardiomyocyte cell states with unique transcriptional signatures and pathway enrichment ([115]Fig. 3b–[116]e and [117]Supplementary Table 8). Fig. 5 |. RUNX1 is downregulated in fibroblasts in cardiac recovery. Fig. 5 | [118]Open in a new tab a, Fibroblast pseudobulk sample PCA colored by condition. b, UMAP of fibroblast cell states. c, Heat map of marker genes for distinct fibroblast cell states. d, Cell composition of fibroblast cell states across five conditions (left) and nuclei density in four conditions (right). e, Pseudobulk expression heat map of canonical genes upregulated and downregulated in fibroblasts in HF across five conditions. Upregulated (f) and downregulated (h) pseudobulk recovery signature ridge plot split across five conditions. GO biological processes pathways enriched in recovery (g) and pathways down in recovery (i). j, Heat map of pseudobulk DE analysis between RR-pre and U-pre split across five conditions. k, Volcano plot of pathways enriched in unloaded group pre-LVAD implantation where each point is a gene set pathway, and the color of the points represents the degree of statistical significance. l, Transcription factors down in cardiac recovery from ENCODE/ ChEA consensus; x axis is transcription factors, and red dots are transcription factors that are statistically significant. Heat map showing expression of RUNX1 target genes across five conditions. m, PRO-seq coverage in unstimulated and TGF-β-treated in vitro fibroblasts ([119]GSE15582) at the RUNX1 locus. n, Density plots of HF genes in UMAP embedding (left) and linear regression of RUNX1 pseudobulk expression and patient ejection fraction in pre and post cohorts (right). o. RNAscope in situ hybridization of representative ×10 fields from a donor, pre-LVAD HF, U-post and RR-post sample (left); n = 23 biologically independent samples and quantification of number of POSTN^+ cells per ×10 field across the four conditions (right). P values were calculated using unpaired t-test with Welch’s correction. Donor versus pre-LVAD HF (**P = 0.0051). Error bars are mean ± s.e.m. p, Linear regression of RUNX1 pseudobulk expression in fibroblasts and inflammatory signature in macrophages in post-LVAD cohort. Dotted line indicates 95% confidence interval; R^2 indicates goodness of fit; and P value indicates whether the slope is significantly non-zero using an F-test. P values were calculated using two-tailed linear regression Wald test with t-distribution in n and p. P values were calculated with Fisher exact test in g, i, k and l. BP, biological process; EF, ejection fraction; PC, principal component; TF, transcription factor. Fig. 3 |. Cardiomyocytes do not revert to a healthy state in cardiac recovery. Fig. 3 | [120]Open in a new tab a, Cardiomyocyte pseudobulk PCA representation of each patient sample colored by condition. b, UMAP embedding plot of cardiomyocytes. c, Heat map of marker genes for distinct cardiomyocyte cell states. d, Dot plot of gene set z-scores for the top ten genes in each cardiomyocyte cell state. e, Enriched GO pathways for cardiomyocyte cell states using statistically significant marker genes identified using a Wilcoxon rank-sum test (adjusted P < 0.05 and log[2] fold change > 0.58). Dot size refers to gene ratio, and color of dots refers to the adjusted P value. P value was calculated using the Fisher exact test. f, Gaussian kernel density estimation of number of nuclei split by condition. g, Pseudobulk heat map of canonical genes up and down in HF split by condition. h, Fluorescence RNAscope in situ hybridization for MYH6 and NPPA in donor, pre-LVAD unloaded, pre-LVAD reverse remodeled, post-LVAD unloaded and post-LVAD reverse remodeled. Images are at ×10 magnification. i, Quantification of number cells per ×10 field for MYH6 and NPPA across the five conditions. P values were calculated using the Brown–Forsythe and Welch ANOVA tests comparing each condition to donor. For MYH6, n = 38 biologically independent samples and Brown–Forsythe ANOVA test F = 28.91, DFn = 4 and P < 0.0001; U-pre (***P = 0.001), U-post (**P = 0.0011), RR-pre (***P = 0.0004) and RR-post (**P = 0.0019) relative to donor. For NPPA, n = 38 biologically independent samples and Brown–Forsythe ANOVA test F = 6.040, DFn = 4 and P = 0.004; U-pre (***P = 0.0001), U-post (**P = 0.0033), RR-pre (non-significant P = 0.086) and RR-post (non-significant P = 0.1017) relative to donor. Error bars are mean ± s.e.m. j, Density plots of MYH6 and NPPA expression in UMAP embedding. k, Upregulated pseudobulk recovery signature ridge plot split across five conditions. l, Transcription factor protein–protein interactions for genes upregulated in recovery. m, Downregulated pseudobulk recovery signature ridge plot split across five conditions. n, Transcription factor protein–protein interactions for CM2 marker genes downregulated in cardiac recovery. PC, principal component. Gaussian kernel density plots showed shifts in cell states among pre-LVAD HF, U-post and RR-post conditions ([121]Fig. 3f). We then created a pseudobulk heat map of canonical genes upregulated and downregulated in HF^[122]13,[123]14 and found that many HF genes are persistently dysregulated in recovered cardiomyocytes. Genes associated with the non-diseased donor state do not return to normal levels in recovery. Of note, reduced expression of NPPA and ANKRD1 was observed in pre-LVAD and post-LVAD implant in patients who recovered compared to those who did not recover, suggesting that the cardiomyocyte substrate may influence the propensity for recovery ([124]Fig. 3g). We validated key findings using in situ hybridization. Compared to donor controls, the number of MYH6-expressing cardiomyocytes was reduced in all pre-LVAD and post-LVAD HF samples. Consistent with our psudobulk analysis, the number of NPPA-expressing cardiomyocytes was significantly increased in U-pre and U-post HF samples compared to donors. Modest trends were observed in RR-pre and RR-post samples ([125]Fig. 3h–[126]i). MYH6 and NPPA expression was enriched in CM0 and CM1, respectively ([127]Fig. 3j). Using the upregulated and downregulated genes in recovery from our pseudobulk analysis, we created a recovery signature for cardiomyocytes that demonstrated selective enrichment in RR-post group ([128]Fig. 3k,[129]m). To examine transcription factors that may modulate recovery, we used EnrichR to find transcription factors predicted to regulate genes upregulated and downregulated in recovery. Notably, we observed upregulation of a GATA4-associated transcriptional network ([130]Fig. 3l) and downregulation of transcription factors associated with inflammation (STAT3, JUNB, JUND and NFKB1) in recovery ([131]Fig. 3n). The gene signature downregulated in recovery was most enriched in the CM2 subset, which expressed canonical HF genes, including ANKRD1 and NPPA ([132]Extended Data Fig. 6a). Among the genes downregulated in cardiac recovery, ABRA was also expressed in CM2 and found to be reduced in patients with HF who recovered pre-LVAD and post-LVAD implant ([133]Extended Data Fig. 6b,[134]c). To validate that reduced cardiomyocyte ABRA expression is predictive of, and associated with, cardiac recovery, we performed in situ hybridization. Quantification of ABRA-expressing cardiomyocytes confirmed reductions in RR-pre and RR-post compared to U-pre and U-post groups ([135]Extended Data Fig. 6d,[136]e). Collectively, these findings suggest that, in recovery, cardiomyocytes show strong transcriptional changes and enter a new state that is distinct from both healthy and HF states. Inflammatory macrophages preclude functional recovery Pseudobulk PCA identified that RR-post myeloid cells clustered independently from donor or pre-LVAD HF myeloid cells, indicating a strong transcriptional signature of recovery in the myeloid compartment. U-post myeloid cells clustered with both RR-post and pre-LVAD myeloid cells, consistent with the finding that recovery transcriptional signatures identify a spectrum of recovery ([137]Fig. 4a). To decipher elements of cardiac recovery encoded in the monocyte, macrophage and dendritic cell compartment, we subclustered and identified nine distinct cell states with unique transcriptional signatures and pathway enrichment ([138]Fig. 4b,[139]c, [140]Extended Data Fig. 7a,[141]b and [142]Supplementary Table 9). Cell composition analysis showed expansion of Mac1 (SPP1 and TPRG1) and Mac2 (PLAUR and FOSB) and reduction of Mac5 in the U-pre and U-post groups relative to donors and RR-pre and RR-post groups ([143]Fig. 4b and [144]Extended Data Fig. 7d). We next plotted the myeloid recovery signature ([145]Fig. 2) in the UMAP space split by donor, pre-LVAD HF, U-post and RR-post groups and observed a strong enrichment in the RR-post group and modest expression in the U-post group. The recovery signature was evident within all myeloid cell clusters ([146]Fig. 4d). To functionally classify the recovery signature, we performed pathway analysis (EnrichR, WikiPathways database) and detected enrichment for EGF, HGF, androgen receptor, oncostatin M and glucocorticoid receptor signaling ([147]Extended Data Fig. 7c). Fig. 4 |. Pro-inflammatory macrophages and RUNX1 are diminished in reverse remodeling, whereas tissue resident macrophages show signs of recovery. Fig. 4 | [148]Open in a new tab a, Myeloid pseudobulk PCA representation of each patient sample colored by condition. b, UMAP of myeloid cell states (left) and cell state composition across conditions (right). c, Heat map of marker genes for distinct myeloid cell states. d, Recovery upregulated signature in four conditions. e, Ridge plot of CD163 expression split by five conditions. f, RNAscope in situ hybridization of representative ×10 fields across conditions (left) and quantification of number of CD163^+ cells per ×10 field (right) (n = 26). Donor versus pre-LVAD HF (***P = 0.0002), donor versus RR-post (non-significant P = 0.1251), donor versus U-post (**P = 0.0034), pre-LVAD HF versus RR-post (***P = 0.0001), pre-LVAD HF versus U-post (*P = 0.0256) and RR-post versus U-post (*P = 0.0108). g, Linear regression of CD163 pseudobulk expression and patient ejection fraction in pre and post cohorts. h, Nuclei density in unloaded pre-LVAD and reverse remodeled pre-LVAD. i, Mac2 composition in donors and pre-LVAD patients (n = 27). Donor versus U-pre (**P = 0.0067), U-pre versus RR-pre (**P = 0.0027) and donor versus RR-pre (non-significant P > 0.145). j, Inflammation score in U-pre and RR-pre (left) and quantified (right) (n = 27). U-pre versus RR-pre (*P = 0.0321), donor versus U-pre (non-significant P = 0.1006) and donor versus RR-pre (non-significant P = 0.1602). k, FDL plot of myeloid cell states with colors from b. l, Palantir pseudotime in FDL space with terminal states (Mac1, cDC2 and Mono2) (top) and nuclei split by four conditions (bottom). m, CHKA, PLAUR and RUNX1 gene expression along pseudotime (left) (shaded area indicating 1 s.d.) in FDL space (right). n, Pseudobulk RUNX1 expression in pre-LVAD HF, U-post and RR-post (n = 40), pre-LVAD HF versus U-post (non-significant P = 0.1177), donor versus RR-post and pre-LVAD HF versus RR-post (****P < 0.0001) and U-post versus RR-post (**P = 0.0058). o, Linear regression of RUNX1 pseudobulk expression and patient ejection fraction in pre and post cohorts. Dotted line indicates 95% confidence interval; R^2 indicates goodness of fit; and P value indicates whether the slope is significantly non-zero using an F-test. P values were calculated using two-tailed linear regression Wald test with t-distribution in g and o. Error bars are mean ± s.e.m. in f, i and j. P values were calculated using unpaired t-test with Welch’s correction from biologically independent samples in f, i, j and n. EF, ejection fraction; NS, non-significant; PC, principal component. We previously identified a role for cardiac resident macrophages in adaptive remodeling of the heart^[149]25–[150]27. To assess whether cardiac resident macrophages are involved in cardiac recovery, we generated ridge plots and found that CD163 (specific marker of cardiac resident macrophages) expression was decreased pre-LVAD HF groups compared to donor controls. Interestingly, CD163 expression was restored to normal levels in the RR-post group, whereas the U-post group displayed an intermediate phenotype ([151]Fig. 4e). We validated our sequencing findings by performing in situ hybridization for CD163 across patient groups. We observed a marked reduction in the number of CD163^+ cells in pre-LVAD HF samples compared to donors. CD163^+ cells increased to near-normal levels in the RR-post group and modestly increased in the U-post group ([152]Fig. 4f). To delineate whether CD163 expression identifies a continuum of recovery, we performed a linear regression analysis of LV ejection fraction versus pseudobulk CD163 expression at the patient level and found a modest correlation (R^[153]2 = 0.21, P = 0.03) ([154]Fig. 4g). To dissect whether certain transcription cell states of myeloid cells predict recovery before LVAD implantation, we constructed a Gaussian kernel density plot of cell number in the pre-LVAD HF group split by U-pre and RR-pre conditions ([155]Fig. 4h). Intriguingly, we found that the Mac2 cluster was overrepresented in the U-pre condition. Quantification of Mac2 composition across individual patients confirmed expansion of this population in the U-pre relative to donor and RR-pre groups ([156]Fig. 4i and [157]Extended Data Fig. 7d). Mac2 represented an inflammatory population that expressed PLAUR and several chemokines and cytokines ([158]Fig. 4b and [159]Extended Data Fig. 7a,[160]b). To assess differences in inflammatory gene expression between U-pre and RR-pre groups, we constructed and plotted an inflammation gene set score (IL1A, IL1B, TNF, AREG, EREG, CXCL2, CXCL3, CCL3 and CCL4) in UMAP space. The inflammatory signature was localized to the Mac2 population and was present selectively in the U-pre group ([161]Fig. 4j), suggesting that the presence of inflammatory macrophages is a negative predictor of cardiac recovery. To decipher the origin of Mac2, we performed pseudotime analysis using Palantir in the HF samples and identified three terminal monocyte-derived states: Mac2, cDC2 and Mono2 ([162]Fig. 4k,[163]l). Cell density plots showed strong phenotype shifts between each pre-LVAD and post-LVAD group when viewed in a force-directed layout (FDL) space ([164]Fig. 4l). Consistent with the above findings, we observed an enrichment for Mac2 in the U-pre group. Notably, post-LVAD implantation RR-post converged toward a Mac5 phenotype, whereas U-post displayed a divergence toward Mac5 (F13A1 and CD163), Mac2 (PLAUR and FOSB) and Mac4 (IFI44L and MX2) ([165]Fig. 4b,[166]l). In particular, the RR-post group converged toward a phenotype marked by enriched CHKA expression ([167]Fig. 4m). This population largely consisted of CD163^+ cardiac resident macrophages. We next plotted MAGIC-imputed gene expression for CHKA, PLAUR and RUNX1 along pseudotime within the Mac2 lineage. We observed a monotonic increase in PLAUR and RUNX1 expression and decreased CHKA expression as the cell differentiated toward Mac2 ([168]Fig. 4m), highlighting competing differentiation trajectories between HF and recovery. To dissect regulatory changes that may underlie transcriptional signatures of HF and recovery in myeloid cells, we performed transcription factor enrichment analysis with DoRothEA^[169]28 in the U-post and RR-post groups. Transcription factors associated with inflammation (RUNX1, NFKB1, NFKB2, STAT3, ATF2 and JUN) were enriched in myeloid cells from the U-post group ([170]Extended Data Fig. 7e and [171]Supplementary Table 18). Quantification of RUNX1 expression at the patient level via pseudobulk analyses revealed increased RUNX1 expression in the donor, pre-LVAD HF and U-post groups. The expression of RUNX1 was markedly diminished in the RR-post group relative to all other groups ([172]Fig. 4n)—these results suggest that RUNX1 downregulation may modulate macrophage phenotype toward a unique state unlike a healthy or failing heart. We next performed a linear regression for LV ejection fraction versus RUNX1 pseudobulk expression and found a strong negative correlation (R^[173]2 = 0.43, P = 0.001) ([174]Fig. 4o). To assess the interplay between RUNX1 and recovery, we plotted a pseudobulk heat map of RUNX1 target genes that were differentially expressed between U-pre and RR-pre groups and identified numerous pro-inflammatory mediators that were enriched in the U-pre group. Collectively, these analyses highlight the possibility that RUNX1 may prevent recovery by promoting pro-inflammatory gene expression in macrophages ([175]Extended Data Fig. 7f). Cardiac fibroblast remodeling in cardiac recovery Cardiac recovery is associated with reductions in myocardial fibrosis^[176]29,[177]30. To ascertain how cardiac fibroblasts shift during recovery, we first performed a pseudobulk PCA. Donor samples clustered separately from the pre-LVAD and post-LVAD groups. A modest separation was also observed between the RR-post group and the pre-LVAD HF group ([178]Fig. 5a). We then subclustered the fibroblasts into eight cell states ([179]Fig. 5b) marked by distinct transcriptional signatures ([180]Fig. 5c, [181]Extended Data Fig. 8a,[182]b and [183]Supplementary Table 10). Cell composition and density analysis demonstrated that Fib1 (SCN7A) was enriched in donors, whereas Fib3 (POSTN, THBS4 and MEOX1) and Fib7 were enriched in HF (GPC6) ([184]Fig. 5d). Pathway analysis revealed enrichment in genes associated with extracellular matrix remodeling in Fib3 and actin binding in Fib7 ([185]Extended Data Fig. 8c). To assess whether recovered fibroblasts represent a reversion to the normal state, we generated a heat map of genes enriched in donor control and HF fibroblasts at the pseudobulk level. This analysis revealed the presence of both persistently dysregulated genes (SVEP1, FAP, POSTN, GPX3 and APOD) and normalized genes (MEOX1, TGFBR3, ACSM3 and PID1) in recovered fibroblasts ([186]Fig. 5e). Consistent with these findings, in situ hybridization for POSTN revealed persistence of this fibroblast population during recovery ([187]Fig. 5o). We then generated ridge plots of the fibroblast-specific recovery signature and detected an enrichment within the RR-post group ([188]Fig. 5f–[189]h). Gene Ontology (GO) of the top genes upregulated in recovery suggested associations with cytoskeletal organization, glucose homeostasis and receptor tyrosine kinase signaling ([190]Fig. 5g and [191]Supplementary Table 11). Conversely, the top pathways downregulated in recovery included TNF-α/NF-κB and TGF-β signaling ([192]Fig. 5i and [193]Supplementary Table 12). To identify cardiac fibroblast transcriptional changes that predict recovery, we performed pseudobulk DE analysis between the RR-pre and U-pre groups. We detected multiple robust differences between groups, including multiple genes association with inflammation (JUNB, CXCL2, KLF2, KLF4, ISG20 and FOSL2) increased in the U-pre group ([194]Fig. 5j). Pathway analysis of dysregulated genes showed strong upregulation of TNF-α/NF-κB, inflammatory response, TGF-β, IFN-γ and IL-6/ STAT3 pathways ([195]Fig. 5k and [196]Supplementary Table 13). To explore transcriptional programs that may drive recovery, we examined ENCODE/ChEA consensus transcription factors from the ChIP-X database and identified RUNX1 as the most downregulated transcription factor during recovery. We then plotted a heat map of genes predicted to be regulated by RUNX1 during recovery and identified a strong downregulation in the RR-post group ([197]Fig. 5l and [198]Supplementary Table 14). Given the above changes observed in RUNX1 expression, predicted activity and TGF-β enrichment, we leveraged published PRO-seq data from in vitro fibroblasts stimulated with TGF-β to assess the level of active transcription around the RUNX1 locus. We found increased RNA polymerase activity at the RUNX1 locus after TGF-β stimulation relative to an unstimulated control ([199]Fig. 5m), suggesting that TGF-β signaling may regulate RUNX1 expression in cardiac fibroblasts. Next, we sought to delineate RUNX1 expression across fibroblast cell states. Density plots revealed that RUNX1 was expressed in Fib3 along with genes known to contribute to myocardial fibrosis, including POSTN, FAP and MEOX1 (refs. [200]22,[201]31,[202]32). Linear regression of ejection fraction versus pseudobulk RUNX1 expression showed a negative correlation (R^[203]2 = 0.15, P = 0.07) ([204]Fig. 5n). To explore a link between macrophage inflammatory gene expression and associated RUNX1 expression in fibroblasts, we performed a linear regression of fibroblast-specific RUNX1 expression versus the macrophage inflammatory score ([205]Fig. 4n) and observed a correlation at the patient level (R^[206]2 = 0.38, P = 0.04) ([207]Fig. 5p). These results suggest that elevated inflammation in macrophages in U-pre relative to RR-pre may hinder cardiac recovery by modulating RUNX1 expression in fibroblasts. RUNX1 GRN is associated with recovery To decipher the potential contribution of RUNX1 toward promoting the recovery phenotype, we performed the following analyses. First, we validated changes in RUNX1 expression in stromal cells across groups using in situ hybridization. Quantitation of the number of RUNX1^+ stromal cells per ×10 field revealed that there were greater numbers of RUNX1^+ cells in the pre-LVAD HF and U-post groups relative to donors. The RR-post group displayed reduced numbers of RUNX1^+ cells relative to the pre-LVAD HF groups and was similar to donors ([208]Fig. 6a,[209]b). Next, we assessed whether RUNX1 target gene expression ([210]Supplementary Table 15) is predictive of cardiac recovery using deep learning. We split the pre-LVAD HF group (U-pre and RR-pre) into a training set and a test set, built and trained a Keras classifier model (detailed in the [211]Methods section) and applied the model in the test set to predict whether nuclei from recovered samples were derived from the U-pre or RR-pre patients ([212]Fig. 6c). We also trained an RF classifier as a comparator to our Keras model. We found that both the Keras model and the RF classifier predicted recovery. The Keras model had a higher area under the curve (AUC) than the RF classifier in predicting recovery (myeloid: Keras = 0.945, RF = 0.789; fibroblast:Keras = 0.942, RF = 0.836) ([213]Fig. 6d,[214]e). We then tested whether perturbation of RUNX1 could facilitate recovery in human macrophages and fibroblasts. We used CellOracle^[215]21,[216]33 to build cell-specific GRNs and performed an in silico deletion of RUNX1. To define directional changes in cell fate resulting from RUNX1 perturbation, we constructed a vector field and visualized the data UMAP space ([217]Figs. 4b and [218]5b). We leveraged Palantir pseudotime trajectory analysis to infer the baseline flow of cells between states ([219]Extended Data Fig. 9a,[220]b). To identify cell states enriched and depleted after RUNX1 pertubation, we took the inner product between control and experimental vector fields. This analysis indicated that RUNX1 perturbation in myeloid cells resulted in movement of cells away from the Mac1 and Mac2 states and toward the Mac5 and monocyte states, signifying a predicted block in the differentiation of inflammatory macrophages and preservation of cardiac resident macrophages ([221]Fig. 6f–[222]h). RUNX1 perturbation in fibroblasts resulted in movement of cells away from the Fib3 and Fib6 states and toward the Fib0, Fib2 and Fib7 states, suggesting a shift away from states associated with fibrosis and HF ([223]Fig. 6g–[224]i). These findings indicate that RUNX1 expression is reduced during recovery and predicts that downregulation of RUNX1 activity in macrophages and fibroblasts signifies the potential for recovery and drives shifts in macrophages and fibroblasts away from states associated with HF. Fig. 6 |. RUNX1 perturbation in silico and in vivo facilitates cardiac recovery. Fig. 6 | [225]Open in a new tab a, RUNX1 RNAscope in situ hybridization across conditions. b, Quantification of a as proportion of positive interstitial nuclei/total interstitial nuclei per ×10 field. n = 37 biologically independent samples and unpaired t-test with Welch’s correction. *P = 0.0188, **P = 0.007 and ***P < 0.001. Error bars are mean ± s.e.m. c, Schematic of machine learning approach used to predict recovery in macrophages and fibroblasts with Runx1 target genes as features. ROC curves with accuracy metrics for test dataset in predicting recovery in macrophages (d) and fibroblasts (e) using a Keras deep neural net classifier model and an RF classifier. CellOracle in silico Runx1 KO simulation quiver plot of vector field in macrophages (f) and fibroblasts (g). Perturbation score with vector field in macrophages (h) and fibroblasts (i). j, Study design of external validation dataset. k, scRNA-seq and scATAC-seq UMAP embedding with cell labels from RNA sequencing label transferred onto ATAC-seq dataset using publicly available data ([226]GSE15582). l, Mouse Runx1 locus showing from top to bottom: ChIP-seq for BRD4 ([227]GSE46668) and H3K27ac ([228]ENCSR000CDF) in adult mouse heart and scATAC-seq in mouse fibroblasts split by four conditions with called peaks ([229]GSE15582) and peak2gene links from ArchR. Numbers above tract indicate ranges of normalized tag densities. Highlighted area is a Runx1 intronic peak. m, Runx1 expression in macrophages (left) and Postn^− fibroblasts (right) in sham, TAC, TAC + JQ1 and TAC + JQ1 withdrawn. Each dot represents a single cell. P values were adjusted with ordinary one-way ANOVA and Turkey’s multiple comparisons test. Myeloid: sham versus TAC + JQ1 withdrawm (*P = 0.0492), TAC versus TAC + JQ1 (***P = 0.0003) and TAC + JQ1 versus TAC + JQ1 withdrawn (**P = 0.0021). Fibroblast: ****P < 0.0001 for all pairwise comparisons. n, Dot plot of human myeloid and fibroblast recovery signature plotted in mouse macrophages in fibroblasts split by four conditions. NK, natural killer; SMC, smooth muscle cell. To evaluate the biological plausibility of the above predictions, we leveraged published in vivo scRNA-seq and single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq) data obtained from a study investigating the potential of BRD4 inhibition (JQ1 treatment) to promote recovery in a mouse model of pressure overload HF (transverse aortic constriction (TAC))^[230]22,[231]31 ([232]Fig. 6j). We re-processed the scRNA-seq data and annotated major cell populations, independently clustered the scATAC-seq data and performed label transfer to annotate scATAC-seq clusters using the RNA information ([233]Fig. 6k). We selected the fibroblasts from the ATAC-seq dataset to explore changes in chromatin accessibility around the Runx1 locus with JQ1 treatment. We then performed peak-to-gene linkage around the Runx1 locus and identified several regions with increased accessibility in TAC relative to sham fibroblasts that decreased with JQ1 treatment and subsequently increased after withdrawal of JQ1. Using publicly available BRD4 and H3K27ac chromatin immunoprecipitation and sequencing (ChIP-seq) data from mouse left ventricle, we validated many of these peaks as active enhancer marks that were binding sites for BRD4 ([234]Fig. 6l). Next, we quantified Runx1 expression from the scRNA-seq data in macrophages and fibroblasts and found that JQ1 treatment was associated with reduced Runx1 expression ([235]Fig. 6m). To compare how JQ1 treatment compares to a simulated Runx1 deletion, we used CellOracle to perform in silico pertubation of Runx1 in myeloid cells and fibroblasts ([236]Extended Data Figs. 9 and [237]10). With JQ1 treatment, there were fewer Il1b^+ macrophages (cluster 0) and relatively more Lyve1^+ resident macrophages (cluster 2) ([238]Extended Data Fig. 9c). CellOracle-simulated deletion of Runx1 predicted a transition away from Il1b^+ macrophages toward dendritic cells ([239]Extended Data Fig. 9d). Within the fibroblasts, JQ1 treatment resulted in a shift away from Postn^+ fibroblasts (cluster 2) toward Apoe^+ fibroblasts (cluster 0), which were most enriched in the sham condition ([240]Extended Data Fig. 10c). CellOracle-simulated deletion of Runx1 predicted a similar transition away from Postn^+ fibroblasts toward Apoe^+ fibroblasts, albeit to a weaker degree as compared to JQ1 treatment ([241]Extended Data Fig. 10d). It is important to recognize that JQ1 inhibits bromodomain proteins, such as BRD4, which control several transcription factors implicated in fibrosis, including Runx1 ([242]Fig. 6l), Meox1 and others. Thus, JQ1 treatment may have a greater effect than Runx1 deletion alone. Finally, to assess the relationship between JQ1 treatment and human signatures of recovery, we examined the expression of human macrophage and fibroblast genes associated with recovery across the mouse conditions. Within macrophages, the human recovery signature was enriched in the JQ1-treated group and not in the sham group, indicating that BRD4 inhibition drives the acquisition of the recovered state. The transcriptional signature of recovery was evident in both sham and JQ1-treated conditions, highlighting that BRD4 inhibition promotes both the acquisition of the recovered state and reversion to baseline in cardiac fibroblasts ([243]Fig. 6n). Discussion By performing snRNA-seq on donor hearts and samples from patients with HF who either recovered or did not recover LV systolic function after LVAD implantation, we defined the single-cell landscape of human cardiac recovery. This work elucidated that the recovered heart represents a distinct biological entity rather than simply a reversion to the normal or non-failing state. We discovered that the transcriptional signature of cardiac recovery is uniquely encoded across cell types with dominant contributions from cardiac macrophages and fibroblasts. Pro-inflammatory macrophage and fibroblast states were found to be negative predictors of recovery. Construction of GRNs predicted that downregulation of RUNX1 transcriptional activity in cardiac macrophages and fibroblasts may facilitate cardiac recovery, a possibility supported by machine learning and in silico transcriptional factor perturbation techniques and corroborated in a mouse model of cardiac recovery mediated by BRD4 inhibition. For patients suffering from HF, cardiac recovery represents the ideal outcome. Recovery of LV systolic function is observed in select individuals after the initiation of anti-remodeling medications (beta-adrenergic, angiotensin II and neprilysin inhibitors and aldosterone antagonists) and mechanical unloading. The frequency of cardiac recovery declines with the severity of HF. Individuals with end-stage HF who are eligible for LVAD implantation and heart transplantation display the lowest incidence of recovery (1–2%)^[244]3–[245]6. It has been postulated that combined mechanical unloading and pharmacologic treatment may serve as an approach to increase recovery rates in the growing advanced HF population where treatment strategies remain limited by recipient eligibility for heart transplantation and donor organ availability^[246]34–[247]37. An improved understanding of the mechanistic basis by which hemodynamic changes imposed by mechanical unloading translate to molecular and structural changes associated with cardiac recovery^[248]38,[249]39 is essential for the development of new therapeutic strategies to increase the frequency of cardiac recovery. Thus, it is of the utmost clinical importance to delineate the cell types and mechanisms that orchestrate acquisition and maintenance of cardiac recovery. The concept that recovery is not merely a reversion to normal is supported by both clinical and biological data. Clinical trials investigating outcomes after withdrawal of anti-remodeling therapy and mechanical unloading in patients with HF who recovered LV systolic function identified a significant rate of recurrent HF and deterioration of cardiac function. These data suggest that cardiac recovery represents a state of remission and indicate that maintenance of cardiac recovery requires ongoing intervention^[250]40,[251]41. The possibility that the recovered heart is not normal is also supported by physiological and pathological investigations demonstrating abnormal contractile reserve, blunted force frequency relationship, and persistent fibrosis^[252]4,[253]9,[254]42–[255]46. Consistent with our findings, bulk RNA sequencing in a mouse model of cardiac recovery triggered by mechanical unloading identified a transcriptional signature that was distinct from both baseline and failing conditions. Although some HF genes normalized, most remained persistently dysregulated, and new signatures associated with recovery emerged^[256]47. By generating a cell-type-specific transcriptional map of human cardiac recovery, we similarly observed that the recovered human heart retains signatures of HF and, instead, found that recovery was associated with the emergence of cell-specific transcriptional states that were not found in healthy or diseased conditions. Intriguingly, we found that strongest transcriptional signatures of cardiac recovery were predominately encoded outside the cardiac myocyte compartment, with the greatest contribution from cardiac macrophages and fibroblasts. We validated this finding in an independent bulk RNA sequencing dataset^[257]20. Notably, previous studies explored the physiological changes that impact cardiomyocyte remodeling in mechanical unloading and cardiac recovery and found that cardiomyocytes undergo selective remodeling^[258]48; however, there is an incomplete reversal of myocyte contractile function despite signs of physiological recovery^[259]49. Within the myeloid compartment, we detected distinct contributions from monocyte-derived and cardiac resident macrophages, populations with opposing functions and effects on HF pathogenesis^[260]25,[261]50,[262]51. Patients who did not recover displayed an enrichment in monocyte-derived macrophages expressing genes implicated in myocardial inflammation and HF (PLAUR, IL1B, TNF and CCL4)^[263]13. The presence of this population was a negative predictor of recovery. These findings are consistent with previous work identifying an association between the abundant pro-inflammatory monocyte-derived macrophages that expressed CCR2 and failure to recover^[264]50. In contrast, we observed that CD163^+ cardiac resident macrophages are lost in HF and revert to normal levels in those who recover. Given that cardiac resident macrophages strongly express transcriptional signatures associated with cardiac recovery and are known to possess functions important for tissue healing and remodeling^[265]27,[266]52,[267]53, it is likely that this population is a central mediator of the recovery process. Cardiac fibroblast activation is an important mechanism of myocardial fibrosis and HF pathogenesis^[268]32,[269]54. The exact role of fibroblasts in cardiac recovery continues to be debated. Within the recovered heart, cardiac fibroblasts retained signatures of fibroblast activation, including the expression of FAP, MEOX1 and POSTN^[270]13,[271]14,[272]22,[273]32. Interestingly, cardiac recovery was associated with reduced expression of genes associated with inflammatory and TGF-β signaling in cardiac fibroblasts and markers of immune-cell-driven fibroblast activation^[274]22,[275]55. When viewed together with the observation that pro-inflammatory macrophages persist in patients who do not recover, these findings suggest that resolution of inflammatory signaling between macrophages and fibroblasts may be essential for recovery. These findings are consistent with the paradigm that para-inflammation leads to persistent tissue dysfunction, whereas resolution of inflammation restores tissue function^[276]56. Previous studies expanded on the role of the cardiac microenvironment and cellular crosstalk in modulating cardiac recovery^[277]57. To unravel transcriptional networks along the macrophage–fibroblast axis driving recovery, we performed transcription factor enrichment and constructed GRNs associated with recovery. This analysis identified downregulation of RUNX1 activity as a key feature of recovery. RUNX1 is a transcription factor with important roles in inflammatory responses and hematological cancers^[278]58,[279]59. Studies performed in zebrafish suggest a role in fibroblast activation^[280]60,[281]61. To explore a relationship between RUNX1 in macrophages and fibroblasts in cardiac recovery, we applied a deep learning approach^[282]62,[283]63 and found that RUNX1 target gene expression measured at the time of LVAD implantation predicts acquisition of recovered versus non-recovered states in both cell populations. We then used CellOracle^[284]21 to build macrophage-specific and fibroblast-specific GRNs and ascertained the predicted effects of RUNX1 deletion using in silico transcription factor perturbation. Within macrophages, RUNX1 perturbation resulted in predicted loss of pro-inflammatory macrophages and pathogenic fibroblasts with shifts toward states observed in recovered hearts. We then tested our predictions in a clinically relevant mouse model of cardiac recovery mediated by JQ1 treatment, an inhibitor of BRD4 (ref. [285]22). Using publicly available multi-omic data, we showed that BRD4 binds to the Runx1 enhancer in mouse hearts^[286]31, suggesting that inhibition of BRD4 may disturb the Runx1 GRN. We identified enhancer peaks in fibroblasts linked to the Runx1 gene that display increased accessibility in HF that are diminished by JQ1 treatment, similarly to what was observed at the Meox1 locus^[287]22. Consistent with disruption of RUNX1 activity, JQ1 treatment led to reduced Runx1 expression in macrophages and fibroblasts. Finally, we show that JQ1 treatment was sufficient to trigger the emergence of the human cardiac recovery signature in macrophages and fibroblasts. Collectively, these findings highlight the possibility that RUNX1 inhibition may serve as an approach to facilitate cardiac recovery. Our study is not without limitations. Availability of paired (preand post-LVAD implant) myocardial tissue specimens from patients who recover is extremely scarce given the low incidence of cardiac recovery and the small number of clinical centers with robust pipelines to phenotype and capture this unique patient population. Given the small cohort of patients who recovered, it is difficult to extend our conclusions to the global population of patients with HF with differing etiologies, comorbidities and demographics. Furthermore, our study employs an snRNA-seq approach that differs from previous studies that rely on bulk RNA sequencing. Compared to these studies, a strength of our work is the ability to dissect cellular heterogeneity and cell-specific changes that are often masked in bulk profiles where the greatest differences are dominated by the most abundant cell type. Unraveling cell-state-specific changes driving recovery is important to develop targeted therapeutics and establish targeted hypotheses for future studies. A limitation of snRNA-seq is the detection of fewer total transcripts. To mitigate concerns of low transcript numbers, we use pseudobulk DE techniques to increase our power at detecting rarer genes and ensure that differences are consistent across multiple patients. All samples sequenced in this study were obtained from the anterior and apical LV walls, and our findings cannot be generalized to other regions of the heart, as the myocardium is not a homogenous structure. It is important to note that the in silico approach used to perturb Runx1 within macrophages and fibroblasts does not consider cellular crosstalk. Although an important limitation, CellOracle transcription factor perturbation studies have been validated in vitro and in vivo and used by other groups^[288]21,[289]33,[290]64. We recognize that there is no perfect mouse model to recapitulate cardiac recovery^[291]47, and the pressure overload JQ1-mediated recovery dataset has inherent limitations. This dataset is used to test the human-derived hypothesis in conjunction with the human recovery data to offset this issue. Lastly, although we show that in silico perturbation of Runx1 is predicted to facilitate the recovery phenotype, we recognize that conditional deletion of Runx1 in macrophages and fibroblasts models is necessary to establish causality. In conclusion, we provide a comprehensive single-cell transcriptomic map of human cardiac recovery, establish that cardiac recovery is a biological state distinct from healthy and disease and unravel cell-type-specific signatures of recovery. Furthermore, we identify shifts in macrophage and fibroblast phenotypes driven by a RUNX1 GRN that predict the propensity for recovery. Finally, we provide orthogonal sources of evidence to suggest that disruption of RUNX1 in macrophages and fibroblasts may drive the recovered phenotype and that RUNX1 inhibition could be an effective approach to facilitate cardiac recovery. Methods Ethical approval for human specimens This study complies with all ethical regulations associated with human tissue research. Acquisition of donor samples was approved by the Washington University institutional review board (study no. 201104172). All samples were procured with informed consent from patients obtained by Washington University in St. Louis School of Medicine. No compensation was provided for participation. Donor and patient demographical details can be found in [292]Supplementary Table 1. Cardiac phenotyping was performed at the time of LVAD implantation and explant, and relevant cardiac function metrics are provided in [293]Supplementary Table 2. Clinical baseline data were collected for all patients at the time of LVAD implantation and are provided in [294]Supplementary Table 3. Sample selection clinical phenotyping The recovery cohort of patients was selected to match for the following variables to the best extent: ejection fraction at the time of LVAD implant, sex, age and clinical risk profile. Age-matched and sex-matched donors were then pulled from the Washington University in St. Louis School of Medicine biobank repository. Within the LVAD cohort, patients were assigned as ‘reverse remodeled’ or ‘unloaded’. Donors were selected to age and sex match with the LVAD samples. RR and U samples were chosen such that ejection fraction at the time of LVAD implantation was not different between the two groups. Human single-nuclei isolation and library preparation snRNA-seq Cardiac tissues from LVAD cores at the time of LVAD implant (U/RR-pre) and adjacent to core samples at the time of explant (U/RR-post) from paired patients were flash-frozen using liquid nitrogen. Identical regions from the apex of LV from explanted donors were used. Single-nuclei suspensions were generated as previously described. In brief, flash-frozen sections were minced with a razor blade and transferred to a Dounce homogenizer containing 1 ml of lysis buffer (10 mM Tris-HCl pH 7.4, 10 mM NaCl, 3 mM MgCl[2] and 0.1% NP-40 in nuclease-free water) on ice. Samples were homogenized using five strokes; an additional 1 ml of lysis buffer was added; and samples were incubated on ice for 15 minutes. Samples were then filtered with a 40-μm filter, and the filter was rinsed with 1 ml of lysis buffer. The mixture was then centrifuged at 500g for 5 minutes at 4 °C, resuspended in 1 ml of nuclei wash buffer (2% BSA and 0.2 U μl^−1 RNase inhibitor (Thermo Fisher Scientific, AM2694) in 1× PBS) and filtered using a 20-μm pluriStrainer (pluriSelect, SKU43–50020-03). Filtered solution was centrifuged using the above criteria and resuspended in 300 μl of nuclei wash buffer and transferred into a 5-ml tube for flow cytometry. Subsequently, 1 μl of DRAQ5 (5 mM solution; Thermo Fisher Scientific, 62251) was added, and the sample was gently vortexed and allowed to incubate for 5 minutes before sorting. DRAQ5^+ nuclei were sorted into 300 μl of nuclei wash buffer using a BD FACSMelody (BD Biosciences) with a 100 μM nozzle. Sorted nuclei were then centrifuged using the above conditions and resuspended in nuclei wash buffer for a final target concentration of 1,000 nuclei per micro-liter—nuclei were counted on a hemocytometer. Based on the nuclei concentration, 10,000 target nuclei were loaded onto a Chip G for GEM generation using the Chromium Single Cell 5′ Reagent version 1.1 kit from 10x Genomics. Reverse transcription, barcoding, complementary DNA amplification and purification for library preparation were performed as per the Chromium 5′ version 1.1 protocol at the McDonnell Genome Institute. Sequencing was performed on a NovaSeq 6000 platform (Illumina) at a target read depth of 100,000 at the McDonnell Genome Institute. Global transcriptomic map generation Nuclei FASTQ files were aligned to the whole-genome pre-mRNA reference generated from the GRCh38 transcriptome (with the intron flag included) using Cell Ranger version 3 (10x Genomics). Nuclei were filtered to include those with 1,000 < RNA unique molecular identifier (UMI) count < 10,000 and mitochondrial reads less than 5%. After initial QC, scrublet was run on each sample separately in Python with default parameters to score nuclei, and nuclei with a doublet score greater than 0.2 were excluded from downstream analysis. We then leveraged supervised doublet removal as previously described on a per-cell-type basis before combining all objects. In brief, clusters were annotated into major cell populations; each major cell type was subsetted and re-normalized; and PCA, UMAP embedding, clustering and DE analysis were performed. Subclusters that did not express the gene signature of the cell type or had overlapping genes across different cell types were removed. After contamination, all cell type objects were merged to construct a cleaned object. After QC and doublet removal, downstream analysis was performed in Seurat version 4 (ref. [295]65). The cleaned object was normalized using SCTransform^[296]66 with regressing out mitochondrial percent and RNA UMI counts. We then computed the principal components and used these to integrate all samples with Harmony^[297]67. Informed by the elbow plot, we used 80 components to construct the UMAP embedding, found nearest neighbors and clustered the data at multiple resolutions. We then used the FindAllMarkers function in Seurat to perform DE testing and annotated clusters into distinct cell types based on canonical gene markers. To define cell states within each major cell type, we subsetted the major cell populations; re-normalized, re-clustered and re-computed the UMAP; and annotated cell states. Both the global object and each cell type object were saved and used for downstream analysis and plotting in Seurat. All differential gene expression to identify cell types or cell states was performed using the normalized assay with the FindAllMarkers function and the Wilcoxon rank-sum test with min.pct = 0.1 and logfc.threshold = 0.25. A heat map of the top ten genes was made across cell types/ states. To further substantiate our annotations, we used the top genes in each cell type/state to create a gene set z-score to see separation across clusters. Pseudobulk differential gene expression Pseudobulk differential gene expression was performed using the DESeq2 (ref. [298]68) package. After QC, cells were subsetted for each cell type; raw counts were extracted; raw counts were aggregated to the sample level; data were normalized using a regularized log-transform; a pseudobulk PCA was performed; and DE analysis between conditions of interest via DESeq2 was performed. For pseudobulk DE analysis, we made the following comparisons: (1) pre-LVAD HF (U-pre and RR-pre) versus donor; (2) RR-post versus donor; and (3) RR-post versus pre-LVAD HF (U-pre and RR-pre). Genes were deemed statistically significant if adjusted P < 0.05 and absolute (log[2] fold change) > 0.58. Statistically significant differentially expressed genes from comparisons (1)–(3) were then used to compute cardiac recovery, persistent HF and HF reversed genes within R as defined by the Venn diagrams. Additionally, we constructed UpSet plots to look at gene overlap across different cell type permutations. Continuum analysis with ejection fraction To connect transcriptional changes with patient ejection fraction within the LVAD cohort, we performed a simple linear regression of ejection fraction versus (1) pseudobulk expression of genes of interest in specific cell types, such as CD163 and RUNX1, and (2) cell-type-specific pseudobulk gene set score of recovery upregulated and downregulated genes. Pseudobulk aggregation was done at the sample level as defined above. Regressions were performed in Prism 9 (GraphPad Software); R^2 was computed for goodness of fit; and P value indicates whether the slope is non-zero (F-test). Cell composition and density shift calculations We used R to compute cell type composition across conditions. To assess shifts in cell density within both the global object and within individual cell types, we converted the .rds object to an .h5ad file format and used the scanpy.tl.embedding function, which employs a Gaussian kernel density estimation of cell number within the UMAP embedding. Density values are scaled from 0 to 1 within that category. RNAscope in situ hybridization Flash-frozen LV samples were fixed for 24 hours at 4 °C in 10% neutral buffered formalin, washed in 1× PBS and embedded in paraffin. Paraffin-embedded sections were cut at an 8-μm thickness using a microtome. RNA in situ hybridization was performed using the RNAScope^[299]69 Multiplex Fluorescent Reagent Kit version 2 assay and RNAScope 2.5 HD Detection Reagent as per the protocol—RED and RNAScope 2.5 HD Duplex Assay Kits (Advanced Cell Diagnostics) using probes designed by Advanced Cell Diagnostics. Fluorescent images were collected using a Zeiss LSM 700 laser scanning confocal microscope. The following RNAScope probes from Advanced Cell Diagnostics were used: MYH6 (555381), NPPA (531281), CD163 (417061), POSTN (409181), PLAUR and RUNX1. Chromogenic/bright-field/fluorescent images were acquired using a Zeiss Axioscan Z1 automated slide scanner. Image processing was performed using Zen Blue and Zen Black (Zeiss). Positive cells were counted using either of two approaches: (1) for fluorescent images, the number of positive cells counted per ×10 field; or (2) for chromogenic images, the number of positive interstitial cells / total number of interstitial nuclei per ×10 field. For CD163 and POSTN, donor and DCM samples were used from our previously published study and re-quantified as per the above approach to ensure comparability. Pseudotime analysis Palantir^[300]70 was used to perform pseudotime analysis. To dissect monocyte fate specification in HF and cardiac recovery, we subsetted the myeloid object to include U-pre, U-post, RR-pre and RR-post and excluded the proliferating cells. The normalized count matrix was then exported as a .txt file and loaded into Python. Palantir was then used to compute a PCA and a diffusion map using a multiscale low-dimensional embedding computer with 5 eigenvectors. An FDL was then computed for visualization of trajectories, and MAGIC was used to impute data for visualization. The Palantir simulation was executed with classical monocytes as the starting state, no terminal states specified and 600 waypoints used. Generalized additive models within Palantir were then used to compute gene expression trends along the Mac2 lineage using the plot_gene_trend_heatmaps function. Pseudotime, entropy and FDL embedding values were exported, and subsequent plotting for visualization was performed in R/Seurat. Pathway analysis and transcription factor enrichment Statistically significant differentially expressed cardiac recovery genes as defined above from the pseudobulk section were used to perform pathway enrichment analysis. EnrichR ([301]https://maayanlab.cloud/Enrichr/) was used for pathway analysis. Transcription factor enrichment analysis was performed using DoRothEA or the ENCODE and ChEA consensus transcription factors from ChIP-X in EnrichR. Pathway and transcription factor enrichment values were downloaded as .csv files, and plots were generated in Prism. GO enrichment analysis to compare cell states was done in R using the clusterProlifer^[302]71 compareClusters function using only statistically significant marker genes. Deep learning to predict recovery In the macrophages and fibroblasts, we exported the counts matrix for the pre-LVAD HF group (U-pre/RR-pre). The counts matrix was filtered to include only RUNX1 target genes. The target genes were obtained from low-throughput and high-throughput functional studies ([303]https://maayanlab.cloud/Harmonizome/gene_set/RUNX1/CHEA+Transcrip tion+Factor+Targets). Next, we used the sklearn.model_selection function to split our data into a train set and a test set (with test.size = 0.3). We then used the StandardScaler function from the sklearn library to scale the train and test feature matrices. We used Keras ([304]https://github.com/keras-team/keras) from TensorFlow ([305]https://www.tensorflow.org/about/bib) with a dropout l2 regularization to train a deep neural network to classify nuclei into RR or U category. The following network structure was used in Python3, and all code has been uploaded to our GitHub: model = keras.Sequential([ keras.layers.Flatten(input_shape=(len(df.columns) - 1,)), keras.layers.Dense(16, activation=tf.nn.relu, kernel_regularizer=regularizers.l2(0.001)), layers.Dropout(0.5), keras.layers.Dense(16, activation=tf.nn.relu, kernel_regularizer=regularizers.l2(0.001)), layers.Dropout(0.5), keras.layers.Dense(1, activation=tf.nn.sigmoid),]) The model was compiled using an ‘Adam’ optimizer and a binary cross-entropy loss function. AUC metrics were calculated. To provide an alternative comparison, we also trained an RF classifier using the RandomForestClassifier function from the sklearn.ensemble library ([306]https://scikit-learn.org/stable/about.html#citing-scikit-learn). Both models were applied on the test dataset, and receiver operating characteristic (ROC) curves were plotted to compare the two approaches. Deep learning was repeated several times, and the accuracy metrics were consistent across simulation runs; in each run, a different permutation of train/test was used. GRN construction CellOracle (version 0.10.5)^[307]21 was used to perform GRN analysis in macrophages and fibroblasts separately. Processed data were imported into scanpy as .h5ad files. The dataset was then downsampled to 20,000 nuclei to limit computational memory usage. Highly variable genes (3,000) were kept for GRN construction. A promoter DNA sequences base GRN was initialized for humans. Unscaled counts were used to generate the CellOracle object, and k-nearest neighbor (KNN) imputation was performed. A cluster-specific GRN was then constructed, and the processed GRN was saved for subsequent transcription factor knockout (KO) analysis ([308]Supplementary Tables 16, [309]17, [310]20 and [311]22). In silico transcription factor perturbation For in silico RUNX1 KO simulation in macrophages and fibroblasts, the processed Oracle object and inferred GRNs were loaded. The GRN was then fit using ridge regression models for the simulation. To simulate a KO, the RUNX1 expression was set to 0. Next, we calculated cell state transition probabilities, which are visualized as vectors on a digitized grid. To establish a baseline developmental flow field, we used Palantir as described above to calculate pseudotime values ([312]Supplementary Tables 19 and [313]21). Next, the Gradient_calculator function from the CellOracle library was used to calculate a developmental vector field on the digitized grid. As described in CellOracle, an inner product was calculated between the baseline vector field and the post-KO simulation vector field to assess perturbations in different cell states, to understand which populations are enriched and depleted after RUNX1 KO. All visualization parameters were used as per CellOracle recommendations. External in vivo dataset analysis scRNA-seq analysis. We used publicly available data with annotations from the original manuscript. Specifically, we extracted the fibroblasts (POSTN^+ and POSTN^−) and myeloid cells. Using the human pseudobulk cardiac recovery genes, we calculated a gene signature for these genes in the mouse dataset and compared sham, TAC, TAC + JQ1 and TAC + JQ1 withdrawn. For the purposes of reference mapping, we extracted the raw filtered matrices for all samples, processed the data using the same pipeline that we applied for the human data and built a single-cell map that we annotated using canonical gene markers. This object was then used to annotate the ATAC-seq data. scATAC-seq analysis. Raw fragment files were used to construct arrow files in ArchR^[314]72. In brief, we used ArchR to calculate doublet scores; we filtered doublets; and we only kept cells with a transcription state site (TSS) enrichment score > 10 and number of unique fragments > 10,000. Next, we computed an iterative latent semantic indexing (LSI), integrated data with Harmony across samples, added clusters and computed a UMAP embedding. To annotate cells from the ATAC-seq data, we used the addGeneIntegrationMatrix function with the above processed and annotated scRNA-seq reference. Using these cluster annotations, we made pseudobulk replicates with the addGroupCoverages function. Then, we used the addReproduciblePeakSet function and MACS2 to call peaks. We then subsetted the fibroblasts, used ArchR to compute peak-to-gene links and used browser plots to visualize peak-to-gene links around the Runx1 locus. We split the browser track by condition to assess chromatin accessibility with peaks linked to the Runx1 gene. Within this same locus, we also leveraged publicly available BRD4 and H3K27ac ChIP-seq data^[315]22—.bed files were downloaded, and tracts were plotted using IGV (version 2.13.1). Statistics and reproducibility Statistical significance was calculated in Prism 9. For RNAscope image quantification, five fields were quantified per patient in areas of maximal staining, and representative images from those areas are shown in [316]Figs. 3–[317]6 and [318]Extended Data Figs. 6–[319]8. Hematoxylin and eosin staining was performed on adjacent areas (for which tissue was available) used for snRNA-seq, and representative images are shown in [320]Suplementary Fig. 2. For snRNA-seq analysis, data integration with Harmony was performed to minimize batch effects across libraries. Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Extended Data Extended Data Fig. 1 |. Quality control metrics. Extended Data Fig. 1 | [321]Open in a new tab nCount_RNA, nFeature_RNA, percent.mt, and scrublet doublet score split by (A) condition and (B) cell type. Extended Data Extended Data Fig. 2 |. Global clustering. Extended Data Fig. 2 | [322]Open in a new tab (A) Heatmap of top10 marker gens for each cell type identified via DE analysis. (B) DotPlot for cell type gene set scores from (A) where the x-axis is cell type gene signature and y-axis is the cluster. (C) Gene set z-scores for top gene markers for each cell type plotted in the UMAP embedding. (D) Cell type composition for each of the patient samples. Extended Data Extended Data Fig. 3 |. Pseudobulk DE analysis to unravel cardiac recovery. Extended Data Fig. 3 | [323]Open in a new tab (A) Pseudobulk DE analysis in each cell type in 3 comparison groups: pre-LVAD HF vs donor, RR-post vs donor, and RR-post vs pre-LVAD HF. Red dots indicate statistically significant genes (adjusted p-value < 0.05). (B) Total number of statistically significant (adjusted p-value < 0.05 and log2FC > 0.58) per cell type in comparisons from (A). (C) Number of overlapping genes in five major cell populations which are up and down in the comparisons from (A). Red number is the number of cardiac recovery genes. P-values calculated using Wald test adjusted for multiple corrections. Extended Data Extended Data Fig. 4 |. Cardiac recovery overlap amongst cell types. Extended Data Fig. 4 | [324]Open in a new tab UpSet plot showing overlap in cardiac recovery genes from ([325]Fig. 2) in five major cell populations which are (A) up and (B) down in cardiac recovery. Extended Data Extended Data Fig. 5 |. Cell-specific pseudobulk analysis. Extended Data Fig. 5 | [326]Open in a new tab Pseudobulk PCA analysis in each cell type colored by five conditions (donor, U-pre, U-post, RR-pre, and RR-post). Extended Data Extended Data Fig. 6 |. ABRA expression enriched in unloaded group. Extended Data Fig. 6 | [327]Open in a new tab (A) DotPlot of cardiomyocyte specific recovery up- and down signature grouped by CM cell states. (B) Density plot of ABRA expression in UMAP embedding. (C) DotPlot of ABRA expression in cardiomyocytes grouped by condition. (D) RNAscope images of ABRA in 5 conditions and scale bar is 100 um. (E) RNAscope images quantified across an array of patients. N = 37 biologically independent samples and p-values calculated using Wald test adjusted for multiple corrections; donor vs U-pre (*P = 0.023), U-pre vs RR-pre (***P < 0.0001), U-pre vs RR-post (***P = 0.0003), U-post vs RR-pre (***P = 0.0007), U-post vs RR-post (*P = 0.0188), and RR-pre vs RR-post (***P = 0.0007). Extended Data Extended Data Fig. 7 |. Macrophage diversity in recovery. Extended Data Fig. 7 | [328]Open in a new tab (A) Gene set z-scores for top gene markers for each cell state plotted in the UMAP embedding. (B) Enrich GO using compareclusters from cluster Prolifer across macrophage cell states. P-value calculated using Fisher exact test. (C) WikiPathways enriched in cardiac recovery. P-value calculated using Fisher exact test. (D) Paired comparison of Mac 2 cluster composition at patient level split by U and RR group from biologically independent samples. (E) DoRothEA TF enrichment analysis in U-post and RR-post zoomed in on some key differentially enriched TFs. (F) Overlap between Runx1 target genes and DE genes between U-pre and RR-pre with heatmap of respective genes split by condition. Extended Data Extended Data Fig. 8 |. Fibroblast diversity in recovery. Extended Data Fig. 8 | [329]Open in a new tab (A) Gene set z-scores for top gene markers for each cell state plotted in the UMAP embedding. (B) DotPlot for cell type gene set scores from (A) where the x-axis is cell type gene signature and y-axis is the cluster. (C) Enrich GO using compareclusters from cluster profiler across fibroblast cell states. P-value calculated using Fisher exact test. Extended Data Extended Data Fig. 9 |. CellOracle simulation in myeloid cells in TAC. Extended Data Fig. 9 | [330]Open in a new tab (A) Myeloid cell states, (B) Marker genes for cell states, (C) Cell state composition and cell density plots in TAC and TAC + JQ1, and (D) Cell oracle Runx1 KO perturbation score with vector field. Extended Data Extended Data Fig. 10 |. CellOracle simulation in fibroblasts in TAC. Extended Data Fig. 10 | [331]Open in a new tab (A) Fibroblast cell states, (B) Marker genes for cell states, (C) Cell state composition and cell density plots in TAC and TAC + JQ1, and (D) Cell oracle Runx1 KO perturbation score with vector field. Supplementary Material supplementary figures [332]NIHMS1915768-supplement-supplementary_figures.pdf^ (151.9MB, pdf) supplementary tables [333]NIHMS1915768-supplement-supplementary_tables.xlsx^ (8.3MB, xlsx) Acknowledgements