Abstract
Heart failure with preserved ejection fraction (HFpEF) as a high
heterogeneity clinical syndrome, is commonly associated with diastolic
dysfunction, and has no effective therapy, which is obviously distinct
from Heart failure with reduced ejection fraction (HFrEF). Currently
the differences of cell type heterogeneity between HFpEF and HFrEF
remain largely unknown. Here we illustrate an atlas consisting of
21,747 cardiac cells, including both HFpEF and HFrEF. Cell-cell
communication analysis reveals cardiomyocytes rather than endothelium
or fibroblasts were dominant communication “hub” in HFpEF. The subtypes
of cardiomyocytes are highly heterogeneous between HFpEF and HFrEF.
Notably, a specific subtype of cardiomyocytes shows significant gene
expression associated with the metabolism of fatty acids. Additionally,
regulon analysis reveals that Ppargc1a, Atf6, E2f6, and Mitf exhibited
specific elevated regulation in the subtype of cardiomyocytes of HFpEF.
Furthermore, we have identified 210 HF susceptibility genes from
HF-associated GWAS data. After integrating scRNA-seq, GWAS, and eQTL
data, the genetic susceptibility underlying HFpEF and HFrEF were
discussed. In conclusion, this study not only comprehensively
characterizes the differences of cardiomyocytes changes but also
provides insights into potential targets for cell type- and
subtype-specific molecules between HFpEF and HFrEF.
Subject terms: Heart failure, Gene regulatory networks
__________________________________________________________________
Single-cell RNA sequencing reveals that cardiomyocyte heterogeneity
drives the divergence between HFpEF and HFrEF, identifying metabolic
subtypes and regulatory factors. Integrated analysis of GWAS and
scRNA-seq data highlights cell-type and subtype specific susceptibility
genes
Introduction
Heart failure (HF) is a complex and serious heterogeneous clinical
condition^[40]1 and a major cause of cardiovascular
hospitalization^[41]2, impacting over 64 million individuals
worldwide^[42]3. In 2012, the global economic burden of HF was
estimated to be $108 billion annually across 197 countries^[43]4. HF is
classified into several subgroups based on left ventricular ejection
fraction (EF), including HF with reduced EF (HFrEF; EF < 40%), HF with
preserved EF (HFpEF; EF ≥ 50%), along with two additional forms of
HF^[44]1. Epidemiological studies have reported that HFrEF or HFpEF are
the most common subtypes, making up more than 85% of all HF
cases^[45]5. In recent years, there have been growing recognition of
the pathological mechanism underling HFrEF, particularly the activation
of the sympathetic nervous system and the renin-angiotensin-aldosterone
system^[46]6,[47]7, which contribute to the pathologic left ventricular
remodeling and cardiac fibrosis^[48]8. The primary classes of
medications used to treat HFrEF have been shown to effectively target
these compensatory pathways above. For a long time, HFpEF was
considered as an early stage of HFrEF, characterized by cardiac
hypertrophy, myocardial fibrosis, and hypertension^[49]9–[50]11.
However, numerous clinical trials assessing the HFrEF-targeted
medications in patients with HFpEF, such as TOPCAT trial (aldosterone
inhibitors)^[51]12, PARAGON-HF trial (angiotensin receptor
blockers)^[52]13, and PRESERVE-HR trial (beta blockers)^[53]14, have
yielded unfavorable outcomes. Researchers now widely recognize HFpEF as
a distinct clinical entity with considerable unmet medical need.
Therefore, it is crucial to explore the differences underlying
mechanisms between HFrEF and HFpEF.
HFpEF is currently regarded as a multisystem disorder, as it
simultaneously affects the heart, vasculature, lung, skeletal muscle,
liver, and kidney^[54]15. Autopsy findings from HFpEF patients have
identified the amyloid depositions, alterations in the stiffness of
myofilament protein, and cardiometabolic disorders linked to diastolic
dysfunction in the cardiac muscle^[55]16–[56]18. However, in contrast
to HFrEF, heart transplants are rarely performed for patients with
HFpEF due to associated complications, and very few medical centers
globally conducting in situ heart biopsies^[57]9. Consequently, there
is extremely limited patients’ tissue data available, particularly from
live cardiomyocytes (CMs). A promising, attractive alternative is to
develop animal models to mimic pathogenesis and explore the molecular
mechanisms related to HFpEF. One study indicated that combining a
high-fat diet with the inhibition of constitutive nitric oxide synthase
could exacerbate concomitant metabolic and hypertensive stress in HFpEF
mice, effectively replicating the multisystem, multiorgan, and
multicellular cardiovascular characteristics observed in human
HFpEF^[58]19. Despite these advancements, the cellular and molecular
mechanisms that underlie HFpEF have remained largely unclear.
Single-cell RNA sequencing (scRNA-seq) has become a robust method for
exploring gene expression at the individual cell level, thereby
deepening our understanding of cellular diversity and gene regulation.
Guo and colleagues utilized single-cell transcriptional profiling of
HFpEF mouse hearts to shed light on the role of Angptl4 in facilitating
interactions between fibroblasts (FBs) and angiogenesis^[59]20.
However, this study excluded CMs, which typically exceed 100 μm in
length. As a primary cellular element in the heart, CMs are crucial for
controlling cardiac contraction and relaxation. Currently, sequencing
methods focus on extracting cell nuclei to tackle the challenge posed
by the inability to encapsulate cardiomyocytes in droplets^[60]21.
However, this method results in the loss of most cytoplasmic mRNA data,
leaving the multi-nucleus information and cytoplasmic content of CMs
unverified^[61]22,[62]23. Given the size, complexity, and importance of
CMs, using nanogrid single-cell selection, nanocell imaging system, and
sequencing platform provides a promising avenue for investigating CMs
heterogeneity and target exploration^[63]24,[64]25.
In this study, to surmount the challenges, we isolated CMs and NCMs in
Langendorff-perfusion heart from high-fat diet and L-name induced
HFpEF, and freshly profiled 10,255 cells by scRNA-seq. By integrating
with the published scRNA-seq dataset of HFrEF, we established a
comprehensive mouse cardiac cell atlas of 21,747 cells across two major
HF subtypes. This atlas facilitated the characterization of
heterogeneity between HFpEF and HFrEF at multiple levels and thereby
provided candidates of HF subtype and cell type-specific intervention
targets for potentially therapeutic applications. Additionally, the
genetic susceptibility underlying HF subtypes was deliberated by
combining HF-associated GWAS and eQTL information.
Results
Overview of the single-cell atlas across different HF murine models
Murine adult CMs are sensitive to mechanical stress, pH changes, and
ionic fluctuations, leading them quickly transition into a pathological
state and ultimately undergo cell death once removed from their native
oxygen-rich environment. To comprehensively analyze the cellular makeup
of the murine heart, we utilized bespoke single-cell platforms that
were designed to handle larger cell sizes and select individual viable
cells via an imaging system. Specifically, the nanogrid single-cell
selection platform features a 125 μm large-aperture nozzle that allows
for extraction of intact adult CMs from cardiac cell suspension. The
integrated imaging unit facilitated the concurrent evaluation of the
viability of sorted CMs at a microscopic level, thereby improving the
quality of single cell sequencing data (Fig. [65]1A). Consequently, we
established the HFpEF murine model by applying consistent metabolic and
hypertensive stress through a high-fat diet and nitric oxide synthase
inhibitor NG-Nitroarginine methyl ester hydrochloride (L-name), which
effectively replicated the conditions of hypertension, obesity, glucose
intolerance, exercise intolerance, and other characteristics observed
in HFpEF patients (Supplementary Fig. [66]1). We then applied the
Langendorff perfusion method to extract both CMs and noncardiomyocytes
from ventricular region of mouse hearts and utilized sequencing
nanogrid platform to generate scRNA-seq profiles for a control group
(N = 3) and a HFpEF group (N = 3), with the HFpEF challenges lasting 15
weeks (Fig. [67]1A). We successfully isolated 87.84% CMs and 99.14%
non-CMs of all murine heart cells as individual viable cells
(Supplementary Data [68]1). After filtering out cells with low number
of detected genes and UMIs, we obtained 10,255 single cells exhibiting
high-quality gene expression profiles from HFpEF case-control settings,
with a median of 833,382 reads depth per cell, a median alignment rate
of 66.49% per cell, and a median of 3014 detected genes per cell
(Supplementary Fig. [69]2A, Supplementary Data [70]2). Furthermore,
profiling over 1470 cells per group enabled the identification of a
putative cell type present at 5% significance threshold with 99%
confidence (Supplementary Fig. [71]2B).
Fig. 1. Single-cell analysis of heart failure in HFpEF and HFrEF murine
models.
[72]Fig. 1
[73]Open in a new tab
A Schematic of the study design and workflow. Cardiomyocytes and
noncardiomyocytes isolated from murine heart in HFpEF study were
stained with Hoechst 33,342 and propidium iodide before selection.
Single live cell (Hoechst^+; PI^−) were identified using the Nanowell
imaging system and selected for subsequent experiments. In
collaboration with a publicly available HFrEF scRNA-seq dataset (GEO
accession number: [74]GSE120064), a comprehensive mapping of the mouse
cardiac cell atlas was performed. B Uniform Manifold Approximation and
Projection (UMAP) showing 21,747 single cells from both HFpEF and
HFrEF. Each dot represents a single cell. Cell populations were
identified by the expression of known marker genes. C Dot plot shows
differential marker genes in each cell type. CM cardiomyocyte,
EC endothelial cell, FB fibroblast, GN granulocyte, MP macrophage, T T
cell, B B cell, RBC red blood cell, W week, L-name
Nω-nitro-Largininemethyl ester, TAC transverse aortic constriction, NCM
noncardiomyocyte.
To obtain comprehensive insights into the cellular composition and
molecular underpinnings of HF in various murine models, we integrated
the scRNA-seq data from 11,492 cells in the HFrEF-model study by Ren et
al.^[75]25 (GEO accession number: [76]GSE120064). Canonical correlation
analysis, along with the exclusion of pan-housekeeping genes that
exhibited consistent expression across cells, was performed to mitigate
potential batch or individual bias in the sequencing data
(Supplementary Fig. [77]2C–E, [78]3A, B, Supplementary Data [79]3).
Finally, our murine single-cell HF atlas included 21,747 cells
representing nine major cell types, including CMs, endothelial cells
(ECs), FBs, macrophages (MPs), granulocytes (GNs), endocardial
endothelium cells (EcCs), T cells, B cells, and red blood cells (RBCs),
identified by their specific molecular markers (Fig. [80]1B, C,
Supplementary Data [81]4). The murine heart atlas contained 49.7% CMs,
15.2% FBs, 12.1% ECs, 14.7% immune cells (MPs, GNs, T, and B cells),
and 1.9% EcCs (Supplementary Fig. [82]3C, D, Supplementary Data [83]5).
The Human Heart Cell Atlas study found that the ventricular region of
the heart consists of 49.2% CMs and 15.5% FBs, proportions of which
align closely with our murine atlas findings^[84]26.
Global transcriptomic profiling highlights distinct transcriptional
landscapes in HFpEF and HFrEF
To investigate the transcriptomic differences between HFpEF and HFrEF
in cardiac tissue, we aggregated the scRNA-seq data by sample groups
(HFpEF, HFpEF_Control, HFrEF, HFrEF_Control) and created pseudobulk
expression profiles. Differential expression analysis revealed 1,355
upregulated and 636 downregulated genes when comparing HFpEF to
HFpEF_Control (Supplementary Data [85]6). Pathway enrichment analysis
of these upregulated expressed genes indicated significant enrichment
for striated muscle cell development and fatty acid metabolic processes
in HFpEF pseudobulk profiles (Supplementary Data [86]7), in contrast to
the canonical association of HFrEF with fibrosis and inflammation
pathways reported in previous studies^[87]25.
Given the critical role of mitochondrial function in cardiac
metabolism, we analyzed mitochondrial gene expression across HF
subtypes. In HFpEF, 13 mitochondrial genes showed significant
dysregulation, including upregulation of mt-Nd1 (1.8-fold upregulation)
and downregulation of mt-Rnr2 (0.8-fold downregulation). Conversely,
HFrEF exhibited alterations in 15 mitochondrial genes, such as mt-Ts2
(1.9-fold upregulation) and mt-Te (0.5-fold downregulation)
(Supplementary Data [88]8). Functional enrichment analysis also
revealed subtype-specific pathway perturbations. In HFpEF,
mitochondrial related genes were enriched in lipid metabolism
(NES = 2.56, adjusted P = 0.017), fatty acid metabolism (NES = 2.50,
adjusted P = 0.016), and carbohydrate metabolism (NES = 2.15, adjusted
P = 0.016), whereas HFrEF-associated mitochondrial genes predominantly
mapped to OXPHOS (NES = 2.01, adjusted P = 0.027) and metal ion
homeostasis pathways (NES = 1.85, adjusted P = 0.037) (Supplementary
Fig. [89]4, Supplementary Data [90]9). These findings suggested
divergent mitochondrial adaptations underlie the metabolic
heterogeneity between HFpEF and HFrEF.
CMs played a pivotal role as the dominant communication hub in pathological
progression
To investigate cellular essential function within a single-cell atlas,
Wang et al. developed a cell-cell interaction map and revealed that ECs
engaged in the highest counts of communications with other cell types,
underscoring their central role in cellular crosstalk^[91]24.
Therefore, it was vital to comprehend the intricate network of
cell-cell communication in the heart to identify the main hub that
orchestrated various signaling pathways associated with different forms
of HF. Initially, we utilized CellChat based on mass action models to
assess cell-cell interactions in scRNA-seq datasets obtained from HF
tissues in murine studies of HFrEF and HFpEF. Consistent with previous
finding^[92]27, our cell-cell communication analysis revealed that FB,
as a source cell type, exhibited the most significant alterations in
interaction counts (FB: N = 92) among all cell types (Fig. [93]2A).
Importantly, top ligands secreted by FB were associated with collagen
synthesis and extracellular matrix organization in HFrEF, including
COLLAGEN and FN1 (Supplementary Fig. [94]5A). Moreover, the frequency
and intensity of interaction between CMs and other cell types,
including JAM, GAS FGF, showed significant reductions (Supplementary
Fig. [95]5B, C). In addition, CMs emerged as the dominant communication
“hub” in HFpEF, which exhibited the highest differential upregulated
interaction counts (CM: 82) and interaction strength (CM: 0.0193)
(Fig. [96]2B), suggesting that CMs were highly interconnected and
played a crucial role in signal transmitting and influencing
neighboring cells more than other cell types in HFpEF. CMs were found
to enhance the activation of the VCAM, EPHA, and HSPG signaling
pathways (Supplementary Fig. [97]6), suggesting that HFpEF may
exacerbate chronic inflammation. These findings highlighted the
importance of understanding the pivotal role of CMs in unraveling the
complexities associated with HFpEF.
Fig. 2. Characterization of cardiomyocytes at different heart failures both
HFpEF and HFrEF.
[98]Fig. 2
[99]Open in a new tab
A Total differential interaction numbers of cell-to-cell signals in
cell types in the HFrEF group versus HFrEF control group. The cell type
indicated on the Y-axis and the X-axis acts as the numbers of the
differential signaling pathways. B Differential interaction numbers
(left) and strength (right) of cell-to-cell signals between the HFpEF
control group and HFpEF group. The color blue signifies a reduction in
the number of interactions or strength within the model group, while
red denotes an increase in the number of interactions or strength
within the model group. The cell type indicated on the Y-axis serves as
the initiator of the signaling pathway, while the cell type indicated
on the X-axis acts as the recipient of the signaling pathway. The
histogram depicts the total variation in numbers or strength of
interaction for each cell type involved in sending or receiving
signals. C UMAP projection displaying 10,744 cardiomyocytes isolated
from hearts in both HFrEF and HFpEF studies, with cells labeled
according to their respective cluster numbers. D Dot plot showing
differential marker genes among CM clusters. Pchd7: Protocadherin 7;
Ndufa4: NDUFA4 Mitochondrial complex associated; Nmrk2: Nicotinamide
riboside kinase 2; Tmem176b: Transmembrane protein 176B; Dlc1: DLC1 Rho
GTPase activating protein; Ndufa5: Ubiquinone oxidoreductase subunit
A5; Xirp2: Xin actin binding repeat containing 2; Nppa: Natriuretic
peptide A. E The gene ontology (GO) analysis was conducted in the top
50 specifically expressed genes within CM1, CM2, CM4, CM5 cluster. The
selected top four categories are shown. F The gene ontology (GO)
analysis was conducted in the top 50 specifically expressed genes
within CM6, CM7, CM8 cluster. The selected top four categories are
shown. G The gene ontology (GO) analysis was conducted in the top 50
specifically expressed genes within CM3 cluster. The selected top four
categories are shown. H UMAP plot shows depicting fatty acid metabolism
among CM clusters. The color intensity of the fatty acid metabolism
score indicates the level of cellular fatty acid metabolic activity,
with the red darker colors indicating higher activity.
The subtypes of CMs involved in HFpEF and HFrEF were heterogeneous
By further elucidating the molecular disparities between CMs in HFrEF
and HFpEF, we employed unsupervised clustering to categorize all 10,036
CMs (Fig. [100]2C), which ultimately identified eight distinct
sub-clusters. The eight subtypes comprised CMs derived from various
mice, indicating a strong reliability and reproducibility in our data,
without significant influence from batch effect or individual
variability (Supplementary Fig. [101]7A). Among these sub-clusters, CM1
constituted the largest group of 18.8%, while CM8 represented the
smallest of 4.6%. (Supplementary Data [102]10).
Within these shared sub-clusters, CM1 displayed a remarkable abundance
of genes encoding contractile force-generating sarcomere protein (Ttn),
calcium-dependent strong adhesive molecule (Pcdh7) and calcium-mediated
process (Ryr2). The top 50 enriched marker pathways in CM1 included
those associated with heart contraction and muscle system process
(Fig. [103]2D, E, Supplementary Fig. [104]7B). The CM2 sub-cluster was
enriched by genes encoding components of the troponin complex (Tnni3),
light chain of myosin protein (Myl3) and mitochondrial respiration
chain complex members (Ndufa4). Pathway enrichment in CM2 involved in
generation of precursor metabolites and energy (Fig. [105]2D, E,
Supplementary Fig. [106]7B); Both CM4 and CM5 sub-clusters showed
highly expression of typical endothelial or fibroblast markers, such as
Pecam1 and Cdh5 or Dcn and Vim, related to angiogenesis or
extracellular organization (Fig. [107]2D, E, Supplementary
Fig. [108]7B). Additionally, noncanonical CMs had been identified, with
experimental validation conducted in several studies to confirm their
existence and elucidated their specific roles within the
heart^[109]24,[110]26.
To explore the disparities between CMs in HFrEF and HFpEF, we performed
a comparative analysis utilizing these marker genes alongside existing
literature on single-cell sub-clustering. The cell subtypes common to
both models exhibited marker gene expression profiles closely aligned
with those documented in the literature regarding public CM scRNA-seq
subtypes^[111]26,[112]28. We identified three CMs (CM6, CM7, CM8)
sub-clusters that were specifically elevated in the context of HFrEF
pathology. The CM6 population was significantly more abundant in the
HFrEF control group and was consistent with previous findings
suggesting that CM6 might represent age-dependent ventricular CMs,
characterized by high expression of genes related to cytokine and
chemokine signaling (such as Ccl2, Ccl4, Cxcl2, Ccl7) (Fig. [113]2F,
Supplementary Fig. [114]7B)^[115]24. CM8 constituted 14.5% of CMs in
HFrEF but only 0.8% in HFpEF CMs (Supplementary Data [116]10). Among
CMs, CM7 exhibited the highest expression of the myosin gene Myh7
(TPM-Like value mean: 1.996), cardiomyopathy-associated protein Xirp2
and matrix metallopeptidase 9 (Supplementary Fig. [117]7B), which
played a role in regulating the remodeling of actin cytoskeleton and
extracellular matrix. Top enriched marker pathways in CM7 were also
related to cardiac muscle tissue morphogenesis (Fig. [118]2F). However,
the disparities in myosin expression between CM7 and CM8 were
characterized by the high levels of Myl4 and Myl7 in CM8 (Supplementary
Fig. [119]7B), indicating a different gene program involved in
structure remodeling of myocardial myosin. Additionally, CM8 displayed
increased expression of Nppa (TPM-Like value mean: 1.052), which
encoded atrial natriuretic peptide precursor, indicating cardiac
hypertrophy (Fig. [120]2D). The top 50 enriched marker pathways in CM8
also included those associated with interferon signaling pathway
(Fig. [121]2F), indicating inflammation abnormalities.
CM3 was specifically found in HFpEF with robust expression of
nicotinamide riboside kinase 2 (Nmrk2). Nmrk2 had been observed to be
present in minimal quantities in healthy myocardium, while levels of
Nmrk2 were significantly elevated in cardiac tissue under pathological
condition^[122]29. The result suggested the homeostasis of nicotinamide
adenine dinucleotide was disrupted in HFpEF. In addition, CM3 showed
highly expression of cell structure genes (Tuba4a, Acta2, Actc1) and
fatty acid metabolism-related genes (Acaa2, Fabp3, Ech1) (Supplementary
Fig. [123]7B), indicating significantly alterations in the
contractile-diastolic function and fatty acid metabolism of CM3 in
HFpEF. Notably, through pathway enrichment analysis of 50
differentially expressed genes in CM3, we observed that CM3 was active
in processes such as fatty acid metabolism, lipid oxidation and cardiac
muscle systolic function (Fig. [124]2G). To confirm the HFpEF-specific
role of the CM3 cardiomyocyte subpopulation, we validated its key
markers Nmrk2 and Acaa2 via immunofluorescence. In HFpEF hearts, Nmrk2
signal intensity was significantly higher than in both controls and
HFrEF. Similarly, Acaa2 expression increased in HFpEF compared to both
controls and HFrEF (Supplementary Fig. [125]8A). These findings
corroborate the scRNA-seq-based identification of CM3 as a
HFpEF-enriched subpopulation with distinct metabolic features.
Moreover, to investigate how Nmrk2 upregulation contributes to HFpEF
pathology, we simulated its knockdown using scTenifoldKnk. Pathway
analysis revealed enrichment in ventricular cardiac muscle tissue
morphogenesis, cardiac muscle contraction and response to metal ions
(Supplementary Fig. [126]8B), suggesting Nmrk2 may modulate both
metabolic and contractile functions in CM3 CMs. Furthermore, to
quantify fatty acid metabolic activity at single-cell level, we
calculated that the scores for all active metabolic pathways in CMs
using vision algorithm. Subclusters were then ranked according to their
average metabolic scores. Remarkably, CM3 exhibited the highest
activity in fatty acid-related metabolic process (mean activity value:
0.10) among all CMs, indicating a high energy expenditure in lipid
oxidation (Fig. [127]2H), which might be linked to their specific role
at metastatic sites. These results implicated CM3 in HFpEF-associated
metabolic dysregulation, raising the possibility that fatty acid
metabolism modulation might have therapeutic potential. Further studies
are required to test this hypothesis and evaluate its clinical
translatability.
Regulatory heterogeneity among various cell types
To investigate the transcriptional regulatory networks underlying the
integrated dataset, we applied a gene co-expression and motif-based
approach to identify the active transcriptional factor (TF) and their
associated downstream regulated gene sets, referred to as regulons’.
The binary activity of each regulon (“on” or “off”) was defined based
on the area under the curve (AUC) at individual cell level. We
identified a total of 370 regulons, and each regulon was characterized
by a single predominant TF, with a median of 171 regulated genes
(Supplementary Data [128]11).
In total, seventy-two unique regulons with top regulon specificity
scores (RSS) were prioritized as cell-type typical regulons (cell-type
TRs) (examples in Fig. [129]3A, full list in Supplementary
Data [130]12). Notably, four out of eight cell-type TRs identified in
CMs were reported to be associated with metabolic traits (Ppargc1a,
Nfe2l1, Rxrg, Atf6), while three were correlated with myocardial
hypertrophy or inflammation (Mitf, E2f6, Hmgb3)^[131]30–[132]38. One
prime example was Ppargc1a (78 g), a regulon with second highest RSS in
CMs (Fig. [133]3B). CMs were found to colocalize with cells exhibiting
active Ppargc1a (78 g) AUC (non-zero AUC values); Additionally, CMs
were observed to colocalize with Ppargc1a (78 g) “on” cells (AUC
threshold was 0.12, Fig. [134]3C), 94.91% of those being CMs
(Supplementary Fig. [135]9A). We further confirmed that the target gene
expressions (Ech1, Hspa9, Ldhb, Mfn2, Ndufa4, Sdhb) of Ppargc1a (78 g)
were higher in CMs compared to other cell types (non-CMs) using qPCR
(Supplementary Fig. [136]9B). Besides its predominant expression in CMs
(Fig. [137]3D), Ppargc1a was consistently reported as a regulator of
mitochondrial energy and lipid metabolism primarily in CMs^[138]39. For
other cell types, FBs were found to colocalize with cells exhibiting
active Tcf21 (95 g), and also colocalize with Tcf21 (95 g) “on” cells
(AUC threshold was 0.05) (Supplementary Fig. [139]10). Highly expressed
in FBs, Tcf21 has been identified essential for FBs
development^[140]40, along with other FB TRs were associated with
differentiation, and glycolysis of cardiac FBs^[141]41,[142]42. ECs
were found to colocalize with cells exhibiting active Elk3 (36 g), and
also colocalize with Elk3 (36 g) “on” cells (AUC threshold was 0.06)
(Supplementary Fig. [143]11). With predominant expression in ECs, Elk3
was reported to regulate preservation of endothelial barrier
function^[144]43, along with other EC TRs were associated with hypoxia
adaptation^[145]44. In addition, six additional colocalizations were
documented, including Foxc1 (12 g) in EcCs, Mxd1 (43 g) in GNs, Irf5
(66 g) in MPs, Pax5 (25 g) in B, Eomes (18 g) in T and Gata1 (22 g) in
RBCs^[146]45–[147]51. Our findings underscored the importance of
identifying cell-type TRs to enhance our understanding of the functions
and regulatory roles of different cell types in the HF development,
particularly in CMs and metabolic traits.
Fig. 3. Cell-type and subtype typical regulons (TRs).
[148]Fig. 3
[149]Open in a new tab
A Cell-type typical regulons (TRs) with top regulon specificity score
(RSS) were displayed. X-axis represented nine cell types: cardiomyocyte
(CM), endothelial cell (EC), endocardial cell (EcC), fibroblast (FB),
macrophage (MP), granulocyte (GN), B lymphocyte (B), T lymphocyte (T),
red blood cell (RBC). Y axis represented regulons (with total number of
genes in bracket). Each dot represented a regulon. Dot size represented
RSS, color-coded by Z scores (scaled RSS). B The RSS of all regulons in
CM (Y axis) were displayed with the descending rank by RSS (X axis).
Cell-type TRs were colored by red and Ppargc1a (78 g) with RSS of 0.65
was highlighted. C A representative cell-type TR in CM, Ppargc1a
(78 g), was showcased. The upper left part showed the UMAP distribution
of CM (orange) and other cells (gray). The upper right part displayed
the UMAP distribution of area under curve (AUC) score of Ppargc1a
(78 g) across all cells (higher AUC referred to higher regulon
activity). The lower left part showed the criteria of binary activity
(on/off) of Ppargc1a (78 g), determined by the elbow-based approach
(AUC cutoff of 0.12). The lower right part displayed UMAP distribution
of Ppargc1a (78 g) binary activity (on/off). D The mRNA expression
level of Ppargc1a across all cell types were displayed with mean and
standard error (error bar). E Six elevated TRs in HFpEF CMs. F
Transcriptional regulation network in HFpEF CM. For the left part, four
circles in the inner represented 4 elevated-regulated subtype TRs of
HFpEF in CM: Atf6 (157 g), Ppargc1a (78 g), Mitf (50 g) and E2f6
(4,759 g). They were connected by arrows with direction (arrow pointed
to transcription factor being regulated). For the right part, LSD
enrichment proportion of Ppargc1a (78 g) was 14.10%
(P = 2.94 × 10^−18), Atf6 (157 g) was 14.10% (P = 1.67 × 10^−16, E2f6
(4,759 g) was 1.70% (P = 1.06×10^−65), Mitf (50 g) was 4.00%
(P = 3.78 × 10^−3). Genes aligned with the outer circle were target
genes colored consistent to the four transcription factors in the
middle. LSDs indicated lipid metabolism/systolic & diastolic disorder
genes. *** represented P < 0.001, ** represented P < 0.01. G Four
attenuated TRs in HFrEF CMs. H Transcriptional regulation network in
HFrEF CM. For the left part, four circles in the inner represented 4
attenuated-regulated subtype TRs of HFrEF in CM: Rxrg (512 g), Ppargc1a
(78 g), Mitf (50 g) and Uqcrb (1437 g). They were connected by arrows
with direction (arrow pointed to transcription factor being regulated).
For the right part, LSD enrichment proportion of Ppargc1a (78 g) was
14.10% (P = 2.94×10^−18), Rxrg (157 g) was 9.18% (P = 5.71 × 10^−68,
Uqcrb (1,437 g) was 2.16% (P = 9.91 × 10^−25), Mitf (50 g) was 4.00%
(P = 3.78 × 10^−3). Genes aligned with the outer circle were target
genes colored consistent to the 4 transcription factors in the middle.
I Regulons with top proportion of “on” binary activity in CM3 cluster
of HFpEF.
Regulatory heterogeneity across different models
To identify the regulatory networks associated with HF, we defined 186
and 79 subtype typical regulons (subtype TRs) for HFpEF and HFrEF,
respectively, based on regulon binary activity (Supplementary
Data [150]13, Methods). Given the significance of lipid metabolism and
cardiac systolic & diastolic dysfunction in HF, we defined 140 genes
related to lipid metabolism-systolic & diastolic dysfunction (LSD) and
28 LSD-enriched regulons through LSD enrichment analysis (Supplementary
Data [151]14, Methods).
Among six elevated TRs in HFpEF CMs, four contribute to the formation
of a transcriptional regulation network (Fig. [152]3E, F). With highest
“on” proportion among four groups (Supplementary Fig. [153]12), Atf6
(157 g) regulated Ppargc1a (78 g) and Mitf (50 g), and the latter two
regulated each other (Fig. [154]3F). This regulatory interaction was
further corroborated by previous studies^[155]52–[156]55. In HFrEF CMs,
four attenuated TRs also formed a regulatory network (Fig. [157]3G, H).
In the HFpEF CM network, Ppargc1a and Mitf were target genes of Atf6
(157 g) and E2f6 (4,759 g) (Fig. [158]3F), whereas in the HFrEF
network, they were target genes of Rxrg (512 g) and Uqcrb (1,437 g)
(Fig. [159]3H). Notably, the binary activity of Ppargc1a (78 g) and
Mitf (50 g) exhibited contrasting trends in HFpEF and HFrEF CMs,
indicating contrasting transcriptional regulation driven by their
distinct upstream transcription factors (Supplementary Data [160]15).
Additionally, Ppargc1a (78 g), Atf6 (157 g), E2f6 (4,759 g), Mitf
(50 g), Rxrg (512 g) and Uqcrb (1,437 g) were all LSD-enriched regulons
(P = 2.94 × 10^−18, 1.67 × 10^−16, 1.06 × 10^−65, 3.78 × 10^−3,
5.71 × 10^−68, 9.91 × 10^−25; LSD gene proportion = 14.10, 7.64, 1.70,
4.00, 9.18, 2.16% respectively); The top regulons with the highest “on”
proportion in CM3 were all LSD-enriched regulons, with Ppargc1a (78 g)
was also CM3 TR (RSS was 0.28) (Fig. [161]3I, Supplementary
Fig. [162]13). These findings indicated a strong association between
HFpEF CMs and active lipid metabolism alongside systolic & diastolic
dysfunction, particularly highlighting CM3.
In addition to CMs, an average of 36 subtype TRs were defined for all
cell types (Supplementary Data [163]16), providing a framework for
exploring the regulatory networks in various cell types underlying
HFpEF and HFrEF.
Integrative analysis of HF-associated loci, susceptibility genes, and
metabolism
Integrative analysis of genome-wide association studies (GWAS) and
scRNA-seq deepened our understanding towards HF from a population
perspective. Initially, eight HF-associated GWAS studies were included,
and a total of 210 HF-susceptibility genes related to 91 HF-associated
loci were finally identified through a tailored workflow (Fig. [164]4A,
Supplementary Data [165]17-[166]19) (Methods). HF-associated loci (with
a median of three HF-susceptibility genes) were dispersed across
chromosomes, with chromosomes 1, 2, and 6 housing the majority of these
loci (Loci counts in Fig. [167]4B). While 14 LSD-enriched regulons were
significantly enriched in 22 HF-associated loci (Loci involved / loci
enriched in Fig. [168]4B, Supplementary Data [169]20). The discovery
that one LSD-enriched regulon was enriched in multiple HF-associated
loci suggested these loci might synergistically influence metabolism
traits in HF development.
Fig. 4. Cell-type and subtype HF-susceptibility typical genes (HF-TGs) and
loci.
[170]Fig. 4
[171]Open in a new tab
A Identification workflow of HF-susceptibility genes and loci. B
LSD-related loci enrichment. X axis represented chromosome 1 to
chromosome 22. Y axis in pink represented HF-GWAS loci counts in each
chromosome. Y axis in brown represented percentage of HF-GWAS loci with
genes being members of LSD-enriched regulons. The Y axis in blue
represented percentage of HF-GWAS loci with genes significantly
enriched in LSD-enriched regulons. C Identification workflow of
cell-type and subtype HF-susceptibility typical genes (HF-TGs). D
Overview of cell-type and subtype HF-TGs. Cell-type HF-TGs in CM were
Actn2, Hspb7, Ttn, Mlip and Mlf1; Cell-type HF-TG in MP was Rgs10;
Cell-type HF-TG in FB was Fkbp7; Cell-type HF-TG in GN was Hp; Subtype
HF-TGs in HFpEF were Actn2, Hspb7 and Ttn; Subtype HF-TGs in HFrEF were
Mtss1, Nfia and Mcmbp. CM indicates cardiomyocyte, EC endothelial cell,
EcC endocardial cell, FB fibroblast, GN granulocyte, MP macrophage, B B
lymphocyte, T T lymphocyte, RBC red blood cell. E External validation
of cell-type HF-TGs in human heart. X axis represented six cell types:
VAM (Ventricular Cardiomyocyte), ACM (Atrial Cardiomyocyte),
Endothelial, Fibroblast, Lymphoid (Lymphocyte), Myeloid (Myelocyte). Y
axis presented mRNA expression of ACTN2, HSPB7, TTN, MLIP and MLF1. F
Colocalization of GWAS and eQTL of ACTN2. For the left part, eQTL
signals of ACTN2 colocalized with those of GWAS (PPH[4] = 99.90% for
rs12724121). Pearson correlation (r = 0.53) was shown between
log-transformed P of eQTL (Y axis) and GWAS (X axis). SNPs were
color-coded by R² of linkage disequilibrium (1000 Genomes, EUR, EAS,
AFR) with the lead SNP (red dots) [blue dots represented R^2 < 0.2,
pale blue represented R^2 < 0.4, green represented R^2 < 0.6, yellow
represented R^2 < 0.8, pink represented R^2 < 1.0]. For the right part,
regional association SNPs of GWAS (upper) and eQTL (lower) within ±
300 kb of rs12724121 (red dot) were presented. The horizontal line
indicated genome-wide significant P for GWAS (5.00 × 10^−8) and
genome-wide empirical P for eQTL of ACTN2 from Heart-Atrial Appendage
(6.40 × 10^−5).
To elucidate cell types and subtypes of locations where
HF-susceptibility genes predominately functioned, we prioritized the
overlap of cell-type typical genes (cell-type TGs) and subtype typical
genes (subtype TGs) overlapped with the HF-susceptibility genes for
further analysis (Fig. [172]4C, Methods). Finally, we delineated 8
HF-susceptibility cell-type TGs and 51 HF-susceptibility subtype TGs
(examples in Fig. [173]4D, full list in Supplementary Data [174]21).
Notably, CMs exhibited the highest number of HF-susceptibility
cell-type TGs (N = 5) among all cell types (Fig. [175]4D), indicating
pivotal role of CMs in HF development. The mRNA expression levels of
ACTN2, HSPB7, TTN, MLIP, MLF1 (cell-type TGs in CMs) were consistently
highest levels in CMs of human heart^[176]56, confirming the
reliability of our analysis (Fig. [177]4E, Methods). A similar trend
had been found in other HF-susceptibility cell-type TGs (Supplementary
Fig. [178]14). Furthermore, HFpEF was characterized by an elevation of
HF-susceptibility subtype TGs (85.71% elevation), while HFrEF showed
the opposite trend (80.00% attenuation), suggesting markedly different
expression patterns between HFpEF and HFrEF (Supplementary
Fig. [179]15, Supplementary Data [180]21).
We employed a Bayesian method for the colocalization of GWAS and eQTL
(Methods) to highlight the potential genetic regulatory role of
HF-associated loci. Colocalization analysis found that six
HF-associated loci were colocalized with one or more eQTLs in
heart-associated tissues (Supplementary Data [181]22). For example,
rs12724121 served as an eQTL for ACTN2 (cell-type TG in CMs) in
heart-atrial appendage tissue, which colocalized with HF GWAS signals
(the posterior probability with one shared causal variant
(PPH[4]) = 99.90%) (Fig. [182]4F). This suggested that rs12724121 could
partially explain the effects of the ACTN2 expression and contributed
to the molecular mechanisms underlying HF. For another example,
rs10497529, rs142556838, and rs2220127 were identified as eQTLs for
FKBP7 in FBs, which colocalized with HF GWAS signals (PPH[4] = 99.90,
98.80, and 96.50%, respectively, Supplementary Data [183]22). This
suggested that these SNPs had the potential to regulate FKBP7
expression and susceptibility to HF. Notably, previous studies
nominated ACTN2 as a key factor in HF^[184]57, and our analysis made
cell-type and subtype location clearer where ACTN2 functioned.
Intriguingly, Actn2 was found to be regulated by LSD-enriched regulon
Uqcrb (1437 g), suggesting an alternative metabolic mechanism of Actn2
in CMs of HFpEF. Furthermore, we incorporated the functional study
specifically evaluating the role of Actn2 in H9C2 CMs using
Actn2-targeting shRNA. This intervention induced significant molecular
changes characteristic of diastolic dysfunction, including marked
upregulation of established biomarkers Ryr2, Trpc1, Scn5a, and Trpc3
(Supplementary Fig. [185]16). These results provide functional evidence
supporting ACTN2’s critical regulatory role in pathways underlying
diastolic abnormalities.
Discussion
The identification of disease-associated studies specific to cell types
is crucial to gain insights into pathogenesis and develop novel
therapies for HFrEF. However, the intricate nature of HFpEF presents
challenges for current single-cell RNA sequencing datasets, which often
suffer from complexity and insufficient sample sizes necessary for a
comprehensive understanding the pathogenesis of HFpEF^[186]9.
Therefore, there is a pressing need to investigate HFpEF using resolved
single-cell transcriptional data^[187]58. In this study, we performed
scRNA-seq on cardiac cells isolated from murine hearts exhibiting
concomitant metabolic and hypertensive stress-induced HFpEF, which
predominantly manifests in the male mouse phenotype and faithfully
recapitulates numerous systemic and cardiovascular characteristics of
HFpEF observed in humans^[188]19. Furthermore, we comprehensively
analyzed the scRNA-seq data of HFpEF and HFrEF to elucidate the
specific pathogenesis.
CMs, a pivotal cell type in the heart, are widely acknowledged for
their intricate involvement in both normal and pathological cardiac
functions. In our study, we identified common subclusters (CM1, CM2,
CM4, and CM5) presented in both HFrEF and HFpEF scRNA-seq datasets.
This finding suggests that there are shared molecular characteristics
contributing to the development and progression of HF, regardless of
its etiology. Additionally, the presence of distinct subclusters (CM3,
CM6, CM7, and CM8) indicates potential variations in disease
pathogenesis or treatment responses among different HF subtypes. The
observed segregation between HFrEF and HFpEF CMs further supports the
idea that they represent discrete clinical entities rather than a
continuum within the spectrum of HF. Notably, CM3, the exclusive subset
of CMs found specifically in HFpEF, exhibited impaired fatty acid
metabolism alongside systolic and diastolic dysfunction. This
reinforces the notion that metabolic dysfunction is a pivotal
pathological manifestation of HFpEF^[189]59. It is truly exhilarating
to witness the remarkable potential demonstrated by drugs regulating
metabolism, such as empagliflozin and semaglutide, in the treatment of
HFpEF^[190]60–[191]62. Therefore, focusing on cardiomyocyte metabolism,
particularly with respect to fatty acid metabolism, may indeed be the
primary therapeutic strategy for addressing HFpEF.
Dysregulated transcription factors represent a unique category of
therapeutic targets that influence the dysregulation of gene
expression, particularly the activation of genes related to cardiac
hypertrophy and fibrosis, hallmark features of HF^[192]63. Our data
analysis revealed both common and specific TFs across different cell
types and different HF models. Among them, Atf6 has been demonstrated
to exert a pivotal adaptive role in cardiac hypertrophy by mitigating
protein misfolding^[193]64. Moreover, the contrasting trends observed
in the binary activity of Ppargc1a (78 g) and Mitf (50 g) in HFpEF and
HFrEF CM networks is intriguing. This observation coincides with
elevated activity of upstream Atf6 (157 g) and E2f6 (4759 g) in HFpEF,
while upstream Rxrg (512 g) and Uqcrb (1437 g) exhibited attenuated
activity in HFrEF. These findings indicate the need for further precise
interventions targeting dominant transcription factors altered in
different HF subtypes. We further performed a functional profile of
these regulons within HFpEF and HFrEF CM networks (full list of
top-ranking pathways in Supplementary Data [194]23). Specifically,
Ppargc1a (78 g) was significantly involved in cytoskeleton-related
pathways^[195]65, with their upstream regulons Atf6 (157 g) and E2f6
(4759 g) were associated with adrenergic and insulin-related
pathways^[196]66 in HFpEF CM network. The literature review was also
consistent with this (Supplementary Data [197]24). The use of
TF-related inhibitors has also shown promising therapeutic effects in
HF treatment. Therefore, employing targeted approaches towards TFs,
such as regulating their intrinsically disordered region, utilizing
auto-response inhibitors and developing proteolytic targeting chimeras,
has the potential to enhance treatment strategies for various subtypes
of HF.
The GWAS study on HF provides a novel perspective for discovering and
identifying genetic variants associated with HF. Multiple research
teams have successfully identified risk variants of significant genomic
impact associated with HF in extensive cohorts of individuals with HF
pathology, as well as in the control subjects, thereby enhancing the
scientific rigor and academic merit of this study^[198]67,[199]68. By
integrating the findings from GWAS and our single-cell analysis, we
have successfully identified 30 as a susceptibility gene specific to
HFrEF, while 21 has been pinpointed as a susceptibility gene specific
to HFpEF. Notably, the Alexis Battle research team previously reported
that rs580698 can modulate the enhancer of ACTN2 gene, thereby
contributing to the progression of HF^[200]57. In our study, we have
discovered another variant, rs12724121, which may also play a role in
HFpEF progression by regulating ACTN2 expression levels in CMs. For
non-CMs, rs10497529, rs142556838, and rs2220127 were identified as
eQTLs for FKBP7 in FBs, which colocalized with HF GWAS signals. By
reviewing findings in cell-cell communication analysis, strong
interaction signals between FBs and CMs were observed in HFrEF and
HFpEF (Fig. [201]2A, B). Fkbp family stabilized the cardiac ryanodine
receptor, potentially mitigating pathological calcium leakage from CMs
during HF^[202]69,[203]70. It was hypothesized that under the potential
regulation of HF-associated SNPs - rs10497529, rs142556838, and
rs2220127, elevated Fkbp7 expression in FBs might influence calcium
regulation by strengthening interactions with CMs and influencing
sarcoplasmic reticulum calcium release in the context of HF. However,
further investigations were required to refine this hypothesis. In all,
the identification of HF subtypes and susceptibility genes specific to
certain cells may offer new avenues for precise targeted therapy in the
treatment of HF.
While our study contributes to the understanding of the pathogenesis of
HF, it is essential to recognize its inherent limitations. Firstly,
male mice were chosen for analysis due to the model’s enhanced clarity
in this gender, as the underlying reasons why gender acts as a risk
factor for HFpEF remain elusive. However, we recognize that this
approach limits the generalizability of our findings, as female
individuals often exhibit distinct pathophysiological features in
HFpEF, such as a higher prevalence of obesity-related inflammation and
preserved cardiac output. Future studies should incorporate both sexes
to explore gender-specific differences in HF pathogenesis and
therapeutic responses. Secondly, given that HF is a multifaceted and
progressive chronic ailment, it is imperative to comprehend its
development across multiple time points. Moreover, our analysis solely
focused on alterations in mRNA levels at the single-cell transcriptome
level. However, considering the intricate and diverse regulatory
mechanisms within organisms, it is crucial to integrate other omics
approaches to comprehensively explore the pathogenesis of HF from
various perspectives. For instance, our analysis of fatty acid
metabolism pathways in CM3 subclusters may not fully capture
post-transcriptional regulation, protein expression, or epigenetic
changes, which are essential for understanding the complete regulatory
network underlying HF. Lastly, it is important to note that the GWAS
information utilized primarily derived from blood samples, which may
not accurately reflect actual cardiac expression patterns.
In conclusion, we conducted single-cell sequencing analysis on the
hearts of HFpEF mice, integrating publicly available HFrEF data. Our
study provides a comprehensive characterization of the transcriptional
profile and regulatory factors in different subtypes of HF.
Additionally, we screened for potential susceptibility genes and
variants related to HF utilizing population-based GWAS data and mapping
to specific cell types and HF subtypes based on mouse single-cell data.
The clinical implications of our findings suggest that the pathogenesis
of HFpEF extends beyond endothelial dysfunction and inflammation,
incorporating disorders of fatty acid metabolism and aberrant
transcription factor activation. Collectively, our findings enhance the
understanding of the transcriptional and molecular underpinnings of CMs
in distinct subtypes of HF, thereby unveiling potential therapeutic
targets for diverse HF subtypes.
Methods
Experimental animals
All animal experiments were conducted in accordance with the Guidelines
for the Care and Use of Laboratory Animals. Official approval has been
obtained from the Animal Care and Use Committee of Peking Union Medical
College (Approval Number: ACUC-A02-2024-012), ensuring compliance with
ethical regulations. Wild-type mice (C57BL/6 J) were procured from the
Beijing Vital River Laboratory Animal Technology Co., Ltd (Beijing,
China). Adult male mice aged 8 weeks were utilized. The mice were
maintained on a 12 h light/dark cycle from 6 AM to 6 PM with
unrestricted access to food (D12450J, Research Diet for the chow groups
and D12492, Research Diet for the HFD groups) as well as water. During
the construction of the HFpEF model, L-name (0.5 g/L, Sigma-Aldrich)
was administered in the drinking water, with the pH adjusted to 7.4.
All the data from animal studies were collected and analyzed by an
independent technician who was blinded to the study design and group
assignment.
Echocardiography and Doppler imaging
Transthoracic echocardiography was performed using a VisualSonics Vevo
2100 system (Visual Sonics, Toronto, CA) as previously
described^[204]19. Echocardiography and Doppler imaging were conducted
at 15 weeks after chow or high-fat diet initiation. Anesthesia was
induced by administering 3% isoflurane and then confirmed by the
irresponsiveness to the firm pressure on one of the hind paws. Vital
signs, including body temperature, respiratory rate, and
electrocardiograph, were continuously monitored. The inhalational flow
rate of isoflurane was adjusted to achieve mouse anesthesia while
maintaining their heart rate within the range of 400–480 beats per
minute or 305–390 beats per minute. Systolic function was assessed
using short-axis M-mode scans at the midventricular level while
maintaining a heart rate of 400–480 beats per minute. Diastolic
function measurements were obtained from apical four-chamber views
using pulsed-wave and tissue Doppler imaging at the level of the mitral
valve, with a heart rate maintained at 305–390 beats per minute. Three
cardiac function parameters were calculated: heart rate, left
ventricular EF, and mitral E wave and E′ wave. All mice recovered from
anesthesia without complications. All data were measured at least three
times.
Speckle-tracking echocardiography and strain analysis
Speckle-tracking echocardiography and strain analysis were conducted
using VevoStrain software (VisualSonics, Toronto, CA) to calculate
global strain in B-mode longitudinal dimensions obtained from the
parasternal long-axis view. The selection of B-mode images was guided
by high frame rates and the capability to visualize both endocardium
and epicardial left ventricular wall borders. Average peak global
strain values were calculated from six distinct anatomical segments of
the left ventricle.
Tail-cuff blood pressure measurements
Blood pressure was assessed using the tail-cuff method (CODA, USA) as
previously described^[205]71. Briefly, mice were acclimated to
short-term restraint prior to testing. For noninvasive measurements, an
individual mouse was placed in holders on a temperature-controlled
platform at 37 °C, ensuring steady-state conditions during recordings.
Blood pressure readings were obtained for a minimum of three
consecutive days and averaged from at least 10–15 measurements per
session.
Intraperitoneal glucose-tolerance test
The intraperitoneal glucose tolerance test was conducted by
administering a glucose injection (2 g/kg in saline) following a 6 h
fasting period. Tail blood glucose levels (mg/dl) were measured using a
glucometer at seven time points: before administration (0 min), and
then at 15-, 30-, 45-, 60-, 90-, and 120 min post-glucose infusion.
Exercise exhaustion test
After a three-day acclimatization period to treadmill exercise, an
exhaustion test was conducted on both HFpEF and control groups of mice.
The mice were subjected to uphill running (20°) on the treadmill
(Columbus Instruments), starting at a warm-up speed of 5 m minute^−1
for 4 min, followed by an increase in speed to 14 m minute^−1 for
2 min. Subsequently, every 2 min, the speed was incremented by 2 m
minute^−1 until the mouse reached exhaustion. Exhaustion was defined as
the failure to resume running within 10 s after direct contact with an
electric-stimulus grid. The distance covered during running was
measured.
Langendorff isolation of cardiomyocytes (CMs) and non-cardiomyocytes (NCMs)
in adult mice
CMs and NCMs were isolated from HFpEF and control group mice, following
established protocols utilizing a standard Langendorff device^[206]72.
Briefly, mice were administered 5 IU/g via intraperitoneal injections
to prevent the blood-clot-induced obstruction in coronary arteries.
After fifteen minutes, the mice were anesthetized intraperitoneally
with Tribromoethanol of 0.175 mg/g body weight. Subsequently, the heart
is removed and soaked in a room-temperature calcium-free separation
buffer to promote natural blood clearance. It is then quickly
transferred to a cold calcium-free separation buffer. All blood
vessels, connective tissues, and the surrounding aortic arch outside
the heart were meticulously severed. Then, the heart was securely
fastened to an aortic cannula filled with perfusion fluid. It is
crucial to ensure that the depth of cannula insertion does not surpass
the aortic valve; hence, cannula, and aorta were firmly affixed using a
nylon thread. The heart was then subjected to perfusion with
calcium-free isolation solution (composed of 137 mM NaCl, 5.4 mM KCl,
1 mM NaH[2]PO[4], 1 mM MgCl[2], 20 mM HEPES, 20 mM taurine, 2 mM sodium
pyruvate, 10 mM BDM and glucose adjusted to pH 7.4) at a controlled
rate of approximately 2 ml/min for five cycles. Next, the heart was
perfused with an enzyme buffer for approximately 15 min at a consistent
flow rate of 37 °C. The enzyme buffer consisted of 0.67 mg/ml type II
collagenase (Gibco, 17101015) and 10 mg/15 ml BSA (Sigma,
V900933-100G). Once softened, the heart tissue with the septum was
carefully detached from the cannula and the left ventricle. The heart
tissue was then isolated and finely minced. The digested tissue
underwent thorough pipetting up and down multiple times, ensuring
proper dissociation of cells. Subsequently, CMs were pelleted by
centrifuging the dissociated cells for 1 min at a speed of 100 × g at a
temperature of 4 °C. The resulting supernatant was transferred to a
fresh 15 ml centrifuge tube and subjected to additional round of
centrifugation for 3 min at a speed of 300 × g to collect NCMs.
Live/dead cell staining
CMs and NCMs were stained using the LIVE/DEAD Viability/Cytotoxicity
Kit (Invitrogen) according to the manufacturer’s instructions to assess
cell viability. The cells were suspended in 2 ml of calcium-free
isolation buffer with 10% FBS, and then added with 4 µl of EthD-1
(2 mM, Component B, for staining dead cells) and 1 µl of calcein AM
(4 mM, Component A, for staining live cells). After incubating at 37 °C
for 15 min, the staining process was terminated by adding an equal
volume of PBS. Then, the cells were settled naturally for a duration of
10 min before removing the supernatant. CMs and NCMs were resuspended
in 2 ml of calcium-free isolation buffer with 10% FBS and stained with
Hoechst 33,342 and propidium iodide at 37 °C for 20 min. Cells were
centrifuged again at 100 g for 3 min at room temperature and
resuspended in the required volume of pre-warmed 1× PBS (without Ca^2+
or Mg^2+, pH 7.4). The resuspended cells were transferred into 12-well
culture plates at a density of 80,000 cells per well. Images were
captured immediately using a Zeiss fluorescence microscope (Axio
observer D1; Zeiss).
Nanowell dispensing, imaging, and selection of single CM or NCM
Single cell CMs and NCMs suspensions were dispensing into the ICELL8
chip (Takara Bio USA) resulting in the nanowell with cells. After cell
dispensing, the chip was centrifuged at 300 × g for 5 min. Then, the
nanogrid wells were automatically imaged using Micro-Manager.
CellSelect^TM software filtered nanowells with multiple cells or no
cells and identified single cell based on imaging automatically. Every
nanowell contains an adapter sequence with barcode (16nt) and unique
molecular identifier (UMI, 14nt). Finally, a file containing position
information and cell identities of the nanowell for each chip was
generated.
Reverse transcription (RT) and next-generation sequencing
The sequencing library was prepared as previously described^[207]25.
Briefly, 100 nanograms RNA was used for generating the sequencing
library by the KAPA mRNA HyperPrep Kits (KK8580, KapaBiosystems). Then,
the dNTPs, RT buffer, RT oligo, RT enzyme, and D-Rase-free water were
mixed and automatically distributed to each nanopore in the chip. The
cDNA synthesis was performed according to manufacturer’s instructions,
followed with the purification and quality evaluation of the product.
After quantification, cDNA was selected for library generation.
Finally, all libraries were sequenced using Illumina NextSeq 500
sequencer (FC-404-2005, Illumina) following the 35nt paired-end
sequencing protocol.
Single-cell data pre-processing and quality control
The raw reads were processed using the Perl pipeline script provided by
WaferGen^[208]25, following several main steps: (1) The validity of
reads was verified. The reads with barcode/unique molecular identifier
(UMI)-contained read 1 were retained. (2) The reads were processed to
filter the adapters using cutadapt (v1.8.1) with the following
parameters: -m 20, --trim-n, --max-n 0.7, and -q 20. (3) Reads were
mapped to genomes of mouse, yeast, E. coli, and adapter sequences using
bowtie2 (v2.2.4), and contaminants were removed with fastq screen
(v0.5.1.4). Only the reads mapped to mouse genome were retained. (4)
Reads were aligned to UCSC mm10 genome by STAR (v2.5.2b) and assigned
to Ensembl gene (GRCm38.75) using featureCounts (subread-1.4.6-p1
command line). Cells were included at a minimum threshold of 500
detected genes and 10,000 UMIs per cell; unique read alignment rate >
50%; and mitochondrial read percentage < 65% (mean + 1.5 s.d.). To
mitigate the potential impact of mitochondrial genes on downstream
analysis, reads mapped to 37 mitochondrial genes (gene symbols starting
with “Mt-”) were excluded. The UMI count for each gene was normalized
by using the following formula:
[MATH:
NVg=log2UgC×10000UtC :MATH]
Here
[MATH:
NVg
:MATH]
represents the normalized value of the gene.
[MATH:
UgC
:MATH]
was the UMI count for each gene in a cell and
[MATH:
UtC
:MATH]
was the total number of UMIs in the corresponding cell. Finally, the
value was subsequently subjected to natural logarithm transformation.
Data normalization and integration
The processed data from both HFpEF and HFrEF studies were normalized
using the logarithm-transformation-based normalization by the
‘NormalizeData’ function with the scale factor of 10,000. Top 2000
variable features were selected based on the
‘SelectIntegrationFeatures’ function. A set of anchors between HFpEF
and HFrEF studies was identified based on the top 2000 variable
features and neighbor-based algorithm by ‘FindIntegrationAnchors’
function. Finally, the pre-computed anchors were utilized to integrate
HFpEF and HFrEF studies using the ‘IntegrateData’ function. To
rigorously address batch effects when integrating newly generated HFpEF
data with publicly available HFrEF datasets, we employed Harmony
(v1.0), a probabilistic batch correction algorithm that preserves
biological variation while aligning cells across batches (Korsunsky et
al., 2019). The integration workflow included the following steps: (1)
Log-normalization of raw counts from all datasets; (2) Selection of
highly variable genes (top 2000) using variance-stabilizing
transformation; (3) PCA dimensionality reduction (n = 50 PCs); (4)
Harmony integration with default parameters using ‘disease subtype’ and
‘batch source’ as covariates.
Single-cell scale, dimension reduction, and unsupervised clustering
The clustering was conducted via Seurat R packag. First, the integrated
dataset was scaled and centered features with using a liner model for
the regression by using ‘ScaleData’ function. Second, we computed 30
principal components by running a principal component analysis (PCA)
dimensionality reduction. Third, we computed the k.param shared nearest
neighbors for the integration study with 30 dimensions of reduction by
‘FindNeighbors’ function. The final resolution value was set as 0.8.
Fourth, we applied the Uniform Manifold Approximation and Projection
(UMAP) to reduce dimensions with 30 neighboring points used in local
approximations of manifold structure via the uwot R package (v0.1.16).
Nine cell types, including CMs, endothelial cells (ECs), FBs, MPs, GNs,
EcCs, T cells, B cells, and red blood cells (RBCs) were determined
according to the expression of known markers from published studies
(Supplementary Data [209]4). For instance, cell populations with
significantly higher expression of Tnnc1, Myh6, and Actn2 were
annotated as “CMs”, whereas those highly expressing Col1a2 and Col3a1
were considered as “FBs”. For CMs, the CMs were extracted and
re-clustered separately. The PCA components were based on the
aforementioned computed PCA components. Differentially expressed genes
for each cell type in HFpEF compared to HFpEF_Control are listed in
Supplementary Data [210]25. Subpopulations were identified for CMs
using the ‘FindSubCluster’ function with default parameters. Highly
correlated sub-clusters were merged and identified 8 CM sub-clusters.
Also, the markers of each sub-cluster were identified by FindAllMarkers
function to check the expression heatmap of sub-clusters markers,
identifying of distinct patterns to effectively distinguish between
sub-clusters.
Mitochondrial gene analysis
The mitochondrial gene analysis was conducted performed using by
MitoCarta3.0 database and clusterProfiler R package (v.4.14.6).
Mitochondrial genes were identified using gene symbols starting with
“Mt-”. Differential expression analysis was conducted separately for
HFpEF and HFrEF using FindMarkers, with adjusted p-values < 0.05
considered significant. Gene Set Enrichment Analysis (GSEA) was further
performed on mitochondrial-related pathways using the MitoCarta3.0
database. Activated pathways (adjusted p < 0.05), including lipid
metabolism, fatty acid metabolism, carbohydrate metabolism, OXPHOS and
metals and cofactors were analysis.
Cell-cell interaction network
The cell-cell interaction analysis was conducted performed using by
CellChat R package (v1.6.1)^[211]73. First, we loaded the CellChatDB
signaling molecule interaction database
([212]http://www.cellchat.org/cellchatdb/), which includes the known
structural composition of ligand-receptor interactions, encompassing
‘Secreted Signaling’, ‘ECM-Receptor’, and ‘Cell-Cell Contact’. Second,
we identified the over-expressed signaling genes associated with each
cell group and over-expressed ligand-receptor interactions by the
function of ‘identifyOverExpressedGenes’ and
‘identifyOverExpressedInterac-tions’. We calculated cell-cell signaling
interaction probability and strength using mass action models between
any interacting cell groups. This step is performed independently in
HFpEF and HFrEF studies. Fourth, the “hub” was defined as the cell type
with the highest differential number of interactions.
Single-cell metabolic activity analysis
The inference of metabolic activity at individual cell level was
performed using the scMetabolism (v0.2.1)^[213]74. Two metabolic
reference datasets were used, 85 metabolism pathways from the KEGG
database and 82 metabolism pathways from REACTOME database. First, we
performed the homology-based conversion of murine genes to human genes
by the nichenetr (v2.0.1, ‘convert_mouse_to_human_symbols’ command
line). Second, the metabolic activity assay of individual cells was
calculated by the AUC algorithm based on the expression matrix of
genes. Finally, single-cell metabolic activity of specific pathways
(fatty acid biosynthesis, fatty acid elongation, fatty acid
degradation, fatty acid metabolism) was analyzed.
Computational Simulation of Gene Knockdown
To explore the functional impact of key genes identified in CM3 CMs, we
applied scTenifoldpy (v0.1.3), a single-cell virtual knockdown
framework. Briefly, CM subpopulation’s gene expression matrix was
subsetted, and a gene regulatory network was constructed using
principal component regression. Virtual knockdown of Nmrk2 was
simulated by removing its expression values and recalculating
downstream gene perturbations.
Lipid-systolic dysfunction (LSD) gene selection
First, we identified genes related to fatty acid metabolism and cardiac
contraction and diastole function in GO database. Second, to maximize
the screening of target genes, all significantly different genes in the
HFpEF group and HFpEF Control group (with P < 0.05 and fold change > 0)
were identified. Third, genes overlapped with the genes selected in the
first step were chosen. Finally, 140 Genes related to fatty acid
metabolism and cardiac systolic dysfunction were selected as LSD
related differential genes.
Gene regulatory network inference
The gene regulatory network was inferred using on a gene co-expression
and motif-based assemble pipeline SCENIC (v1.3.1) following 4four main
steps. First, the RNA counts were used as the input. Genes with average
expression > 0.03 and detected in at least 1% of cells were retained.
Second, the potential TF-gene pairs were calculated by a co-expression
analysis using the GENIE3 (v1.24.0). Third, the regulons were
identified by both considering the potential TF-gene pairs and reported
TF binding information from a motif database RcisTarget (v1.19.2).
Fourth, the activity of each regulon in each cell was calculated by the
area under the recovery curve (AUC) and integrated the expression ranks
across all genes in a regulon, using the AUCell (v1.24.0).
LSD associated enrichment
To identify LSD-enriched regulons, odds ratios (OR) were calculated
based on the actual and expected number of genes in both LSD and
regulons (OR > 2). Two-tailed Chi-square test was used to assess the
significance (P < 0.05). To identify LSD-enriched loci, odds ratios
(OR) were calculated based on the actual and expected number of genes
in both HF-associated loci and LSD-enriched regulons (OR > 2). A
two-tailed Fisher’s exact test was employed to evaluate significance
(P < 0.05), as the expected frequencies of genes were less than 5.
Regulon specificity score (RSS) calculation
To quantify cell type typical regulon, regulon specificity score (RSS)
was calculated based on an entropy-based strategy^[214]75 following
three steps:
First, the distribution of regulon activity score (RAS) and normalized
RAS was calculated. For each active regulon, a vector was used to
represent the distribution of RAS in the cell population (
[MATH: n :MATH]
is the total number of the cells, R is regulon):
[MATH:
PR=p1<
mrow>R,⋯,pnR
:MATH]
RAS are normalized so that:
[MATH: ∑i=1
np
iR=1
:MATH]
A vector was used to indicate whether a cell belongs to a specific cell
type (C is cell type):
[MATH:
PC=p1<
mrow>C,⋯,pnC
:MATH]
where
[MATH:
PiC
=1,cellbelongstothespecificcelltype<
/mtr>0,otherwise<
/mtable> :MATH]
The vector is also normalized so that:
[MATH: ∑i=1
np
iC=1
:MATH]
Second, the Jensen-Shannon divergence (JSD), a metric commonly used for
quantifying the difference between two probability distributions, was
evaluated. JSD was defined as:
[MATH: JSDPR,PC=HPR+PC
2−HPR+HPC2 :MATH]
where
[MATH: HP=−∑pil<
mi>ogpi :MATH]
represents the Shannon entropy of a probability distribution
[MATH: P :MATH]
. The range of JSD values is between 0 and 1, where 0 means identical
distribution and 1 means extreme difference.
Third, the regulon specificity score (RSS) is defined by converting JSD
to a similarity score:
[MATH: RSSR,C=1−JSDPR,PC :MATH]
For each cell type, top eight unique regulons with highest RSS were
prioritized as cell-type typical regulons (cell-type TRs).
Defining subtype typical regulons (subtype TRs)
The odds ratios (OR) were calculated based on the actual and expected
number of cells being “on” binary activity in each regulon. A two-sided
Chi-square test was employed, and significance was assessed with
P < 0.05. HFpEF-elevated TRs were defined to meet the criterion that
regulon “on” proportion in HFpEF was simultaneously higher (OR > 1.5)
than in both HFpEF control and HFrEF. HFpEF-attenuated TRs were defined
to meet the criterion that regulon “on” proportion in HFpEF was
simultaneously lower (OR < 0.67) than in both HFpEF control and HFrEF.
HFrEF-elevated TRs were defined to meet the criterion that regulon “on”
proportion in HFrEF was simultaneously higher (OR > 1.5) than in both
HFrEF control and HFpEF. HFrEF-attenuated TRs were defined to meet the
criterion that regulon “on” proportion in HFrEF was simultaneously
lower (OR < 0.67) than in both HFrEF control and HFpEF.
Workflow of HF-susceptibility genes and loci identification
The workflow contains following three steps:
Step 1: HF-associated lead variants identification.
By searching GWAS Catalog and PubMed using key word “Heart failure” or
“HF” (“GWAS Catalog & Published studies” in Fig.[215]4A), 114
HF-associated variants from eight studies with P < 5.00×10^−8 were
originally acquired (“114 variants & 8 studies” in Fig. [216]4A).
Linkage disequilibrium (LD) R² < 0.2 was used to acquire independent
lead variants, and was calculated according to European (EUR), East
Asian (EAS) and African (AFR) population using LDlinkR package
(v1.3.0). The loci of HF-associated lead variants were sourced from the
UCSC Genome Browser database^[217]76. Finally, 110 independent lead
variants & 93 loci were identified to be HF-associated lead variants
and loci (“110 variants & 93 loci” in Fig.[218]4A).
Step 2: HF-candidate variants identification.
HF-candidate variants were calculated based on HF-associated lead
variants above, considering that the functional variants in HF might
not always be the lead variants but variants near them. LD R² > 0.8 was
used to define HF-candidate variants which had potential to be the
functional variants near HF-associated lead variants. Finally, a total
of 1,001 HF-candidate variants & 93 loci were identified to be
HF-candidate variants and loci (“1001 variants & 93 loci” in
Fig.[219]4A).
Step 3: HF-susceptibility genes identification.
To identify HF-susceptibility genes, a ‘various-to-gene’ (V2G)
machine-learning model on Open Target Genetics
([220]https://genetics.opentargets.org/) was utilized. The V2G model
integrated functional genomics information including variant-gene
distance, expression quantitative trait locus (eQTL), and transcription
initiation site (TSS). Top three genes with highest V2G scores for each
HF-candidate variant were identified to be HF-susceptibility genes
(N = 210). Two loci without matched genes were filtered. A total of 91
loci were finally preserved, with 210 HF-susceptibility genes
potentially being regulated (“1,001 variants & 91 loci & 210
HF-susceptibility genes” in Fig.[221]4A).
Defining HF-susceptibility cell-type typical genes (cell-type TGs)
First, cell-type TGs were defined to meet the criteria that
log[2]FC > 2 and detected in at least 50% cells in the group. A
two-sided Wilcoxon rank-sum test was employed, and significance was
assessed using the Bonferroni correction with an adjusted P < 0.05.
Second, by performing intersection with HF-susceptibility genes from
GWAS, overlapped genes were defined as HF-susceptibility cell-type TGs.
Murine genes were converted to their human counterparts using
homology-based methods (Homologene R package, v1.4.68.19.3.27),
followed by manual verification.
Defining HF-susceptibility subtype typical genes (subtype TGs)
First, subtype TGs were identified. A two-sided Wilcoxon test was
employed, and significance was assessed using the Bonferroni correction
with an adjusted P < 0.05. HFpEF-elevated TGs were defined to meet the
criterion that gene expression in HFpEF was simultaneously higher
(log[2]FC > 0) than in both HFpEF control and HFrEF. HFpEF-attenuated
TGs were defined to meet the criterion that gene expression in HFpEF
was simultaneously lower (log[2]FC < 0) than in both HFpEF control and
HFrEF. HFrEF-elevated TGs were defined to meet the criterion that gene
expression in HFrEF was simultaneously higher (log[2]FC > 0) than in
both HFrEF control and HFpEF. HFrEF-attenuated TGs were defined to meet
the criterion that gene expression in HFrEF was simultaneously lower
(log[2]FC < 0) than in both HFrEF control and HFpEF.
Second, HF-susceptibility subtype TGs were identified. By intersecting
subtype TGs above with HF-susceptibility genes from GWAS, overlapped
genes were prioritized to be HF-susceptibility subtype TGs. Murine
genes were converted to their human counterparts using homology-based
methods (Homologene R package, v1.4.68.19.3.27), followed by manual
verification.
Validation in the human heart cell atlas
The human heart scRNA-seq data were downloaded from Human Heart Cell
Atlas database
([222]https://www.heartcellatlas.org/#team/global_raw.h5ad) and were
converted to Seurat file format (SeuratDisk, v0.0.0.9020). Given the
discrepancy of cell-type annotation, cell types of Human Heart Cell
Atlas database were mapped to our dataset based on the following rules:
Ventricular_Cardiomyocyte and Atrial_Cardiomyocyte to CM; Endothelial
to EC; Fibroblast to FB; Lymphoid to B and T cells; Myeloid to GN. RNA
raw counts of genes were fetched and shown.
Colocalization analysis between eQTL and GWAS signals
Colocalization analysis between eQTL and GWAS was performed using a
Bayesian-based framework (Coloc, v5.2.3). The summary statistics for
GWAS of all-cause HF were obtained from the GWAS Catalog database under
accession code GCST90162626^[223]68, containing 115,150 cases and
1,550,331 controls of diverse genetic ancestry. Expression quantitative
trait loci (eQTL) data was obtained from Genotype-Tissue Expression
(GTEx) v8^[224]77 by ezQTL tool
([225]https://analysistools.cancer.gov/ezqtl/#/qtls). “Artery-Aorta”,
“Atery-Cornary”, “Heart-Atrial-Appendage”, “Heart-Left-Ventricle”,
“Thyroid”, and “Whole-Blood” were selected to be HF-associated tissues
with eQTL. Four posterior probabilities were estimated: association
with GWAS only (PPH[1]), association with eQTL only (PPH[2]),
association with GWAS and eQTL but two independent variants (PPH[3]),
and association with GWAS and eQTL having one shared variant (PPH[4]).
Linkage disequilibrium (LD) and minor allele frequency (MAF)
information from African, East Asian and European (AFR, EAS, EUR) were
obtained from the 1000 Genomes database for colocalization. Evidence
for colocalization was determined based on PPH[4] > 0.8 (300 kb region
centered on the lead variants).
Immunofluorescence staining
Mouse heart paraffin sections were deparaffinized and submitted to
antigen retrieval in sodium citrate buffer as previously described.
Afterwards, sections were incubated with Dako serum-free blocking
solution (Dako) for 1 hour at RT. Sections were incubated with the
primary antibodies against Nmrk2 and Acaa2 overnight at 4°C. Some
sections were incubated with isotype controls or only stained with
secondary antibodies. After washing three times with PBS, sections were
treated with secondary antibodies for 1 hour in the dark. Then DAPI was
used to stain the nuclei. Moreover, heart sections were processed with
the TrueVIEW Autofluorescence Quenching Kit (Vectorlabs) according to
the manufacturer’s instructions to remove the cardiac autofluorescence.
After washing, sections were mounted with FluoroshieldTM Aqueous
Mounting medium. Images were acquired in a confocal microscope or
microscope at 20x magnification (Leica Stellaris 5 Confocal
Microscope).
Cell culture
H9C2 rat CMs (Procell) were cultured in high-glucose Dulbecco’s
Modified Eagle Medium (DMEM, Gibco) supplemented with 10% fetal bovine
serum (FBS) and 1% penicillin/streptomycin. Cells were maintained in a
humidified incubator under 5% CO[2] atmosphere at 37 °C, with medium
refreshed every 48 h. Subculturing was performed using 0.25%
trypsin-EDTA upon reaching 80–90% confluency. HEK293T cells were
maintained in high-glucose DMEM supplemented with 10% FBS and 1%
penicillin/streptomycin at 37 °C/5% CO[2]. Cells at 70–80% confluency
in 10 cm dishes were used for transfection.
Rat ACTN2 shRNA design and lentiviral vector construction
Two candidate shRNA sequences (19-21 nt) targeting rat Gene ACTN2 mRNA
(NCBI RefSeq: [226]NM_001170325) were designed using the NCBI RNAi
Design Tool. Sequences were filtered for specificity (BLASTn against
rat transcriptome) and GC content (30–55%). A scramble shRNA with no
homology to mammalian genes served as negative control. Forward and
reverse oligonucleotides containing: shACTN2-1:
5’-CTATGTTCTCGATCTGGGTGC-3’; shACTN2-2: 5’-CAATGTCCGACACCATCTTGC-3’.
pSIH-H1 shRNA cloning vector was linearized with BamHⅠ and EcoRI (NEB).
Annealed shRNA duplexes were ligated using T4 DNA Ligase at 16 °C for
16 h. Ligation products were transformed into E. coli via heat shock.
Positive clones were selected on LB agar plates with ampicillin
(100 μg/mL). Plasmid DNA was extracted and sequenced using the H1
promoter primer. Validated constructs were stored at −80 °C.
Lentivirus production and transduction
For each dish, cells were transfected using Lipofectamine 2000
(Invitrogen) with a plasmid mixture containing 10 μg shRNA construct
(validated pSIH-H1 vector), 7.5 μg psPAX2 packaging plasmid, and 2.5 μg
pMD2.G envelope plasmid. DNA was diluted in 1.25 mL Opti-MEM I Reduced
Serum Medium, combined with 50 μL Lipofectamine 2000 in 1.25 mL
Opti-MEM, incubated for 20 min at room temperature, then added dropwise
to cells. After 6 h incubation, medium was replaced with high-glucose
DMEM supplemented with 10% FBS. Viral supernatant was harvested at 48
and 72 h post-transfection, pooled, filtered through 0.45 μm PVDF
membranes. Concentrated virus was obtained by ultracentrifugation at
120,000 × g for 2 h at 4 °C, with pellets resuspended in 100 μL
ice-cold PBS and stored at −80 °C. For H9C2 transduction, H9C2 cells
were seeded in 6-well plates at 1.5 × 10⁵ cells/well 24 h prior to
infection. Thawed concentrated virus was combined with polybrene
(8 μg/mL final concentration) in DMEM. Cells were transduced with virus
(final MOI = 20). Knockdown efficiency was validated 72 h
post-transduction through quantitative real time-PCR and western
Blotting.
RNA isolation and quantitative real time-PCR
RNA was isolated from H9C2 cells using TRIzol reagent. After
homogenization, chloroform was added, followed by centrifugation
(12,000 × g, 15 min, 4 °C). The aqueous phase was mixed with
isopropanol, centrifuged, and the RNA pellet washed with 75% ethanol
(DEPC-treated water). RNA was dissolved in DEPC-treated water and
quantified via NanoDrop 2000. For cDNA synthesis, 1 μg RNA was
reverse-transcribed using the FastKing cDNA Kit (Tiangen). Quantitative
real-time PCR (qPCR) was performed on the Bio-Rad iQ5 system with
TransStart Green qPCR SuperMix. Reactions (20 μL) contained 10 μL 2×
SuperMix, 3 μL cDNA, 2 μL primers (1 μM each), and 5 μL H₂O. The
thermal profile included 95 °C for 2 min, followed by 40 cycles of
95 °C (15 s) and 60 °C (30 s). Gene expression was normalized to Actin
using the ΔΔCt method.
Western blotting analysis
Total proteins were extracted from H9C2 cells using a lysis buffer
(50 mM Tris pH 7.4, 500 mM NaCl, 1% NP40, 20% glycerol, 5 mM EDTA, 1 mM
PMSF) supplemented with phosphatase/proteasome inhibitors (Cell
Signaling Technology). Insoluble debris was removed by centrifugation
(14,000 × g, 10 min, 4 °C), and supernatants were quantified via BCA
assay (Thermo Scientific). Samples were mixed with 5× loading buffer
and denatured (95 °C for 5 min). Proteins were separated by SDS-PAGE,
transferred to PVDF membranes (Millipore), and incubated overnight at
4 °C with primary antibodies. Membranes were washed (TBST ×3) and
incubated with HRP-conjugated secondary antibodies (goat
anti-mouse/rabbit IgG, 1:10,000) for 1 h at RT. Bands were visualized
using SuperSignal West Pico substrate and quantified with Image Pro
Plus 6.0 software.
Software availability
We performed the analyses using the following software packages:
Cutadapt (v1.8.1), [227]https://github.com/marcelm/cutadapt/; Bowtie 2
(v2.2.4), [228]https://bowtie-bio.sourceforge.net/ bowtie2/index.shtml;
FastQ Screen (v0.5.1.4),
[229]https://github.com/StevenWingett/FastQ-Screen; STAR,
[230]https://github.com/alexdobin/STAR; Rsubread,
[231]https://bioconductor.org/packages/release/bioc/html/
Rsubread.html; Seurat (v5.1.0), [232]https://satijalab.org/seurat/;
clustree, [233]https://github.com/lazappi/clustree; CellChat (v1.6.1),
[234]https://github.com/sqjin/CellChat; MitoCarta3.0,
[235]http://www.broadinstitute.org/mitocarta； CellChatDB,
[236]http://www.cellchat.org/cellchatdb/; scTenifoldpy (v0.1.3),
[237]https://pypi.org/project/scTenifoldpy/; scMetabolism (v0.2.1),
[238]https://github.com/wu-yc/scMetabolism; Nichenetr (v2.0.1),
[239]https://github.com/saeyslab/nichenetr; SCENIC (v1.3.1),
[240]https://github.com/aertslab/SCENIC; RcisTarget (v1.19.2),
[241]https://github.com/aertslab/RcisTarget; GENIE3 (v1.24.0),
[242]https://github.com/aertslab/GENIE3; AUCell (v1.24.0),
[243]https://github.com/aertslab/AUCell; LDmatrix Tool,
[244]https://ldlink.nih.gov/; Open Targets Genetics,
[245]https://genetics.opentargets.org/; LDlinkR (v1.3.0),
[246]https://github.com/CBIIT/LDlinkR; Homologene (v1.4.68.19.3.27),
[247]https://github.com/oganm/homologene; Human Heart Cell Atlas
database, [248]https://www.heartcellatlas.org/#team/global_raw.h5ad;
SeuratDisk (v0.0.0.9020), [249]https://github.com/
mojaveazure/seurat-disk; Coloc (v5.2.3),
[250]https://github.com/chr1swallace/coloc.
Statistics and reproducibility
Data was analyzed by GraphPad Prism 10 (GraphPad Software, LLC, San
Diego, CA) and R software. All results were expressed as
mean ± standard error of the mean (mean ± SEM). P value was calculated
by two-tailed Chi-square test, two-tailed Fisher’s exact test,
two-tailed Wilcoxon test as indicated in supplemental materials;
Adjusted P was calculated based on Bonferroni correction; * represented
P < 0.05, ** represented P < 0.01, *** represented P < 0.001. P < 0.05
was considered statistically significant.
Reporting summary
Further information on research design is available in the [251]Nature
Portfolio Reporting Summary linked to this article.
Supplementary information
[252]Supplementary Information^ (13.5MB, pdf)
[253]42003_2025_8827_MOESM2_ESM.pdf^ (111.4KB, pdf)
Description of Additional Supplementary Materials
[254]Supplementary Data 1-25^ (6.7MB, xlsx)
[255]Reporting Summary^ (2MB, pdf)
Acknowledgements