We collected five independent HF datasets, including GSE135055, GSE203160, GSE48166, GSE198945, and GSE116250, all containing bulk RNA-seq data from cardiac tissues of human HF patients and healthy controls. In addition, we selected our previously obtained human cardiac bulk RNA-seq datasets GSE165303 for validation. Additionally, GSE235897 was included, providing bulk RNA-seq data of human monocyte-derived macrophages (hMDMs) under trained immunity conditions. All the above datasets were obtained from the GEO public database (https://www.ncbi.nlm.nih.gov/geo/). For single-cell transcriptomic data, we accessed the SCP1303 project from the Broad Institute (https://singlecell.broadinstitute.org/single_cell), which includes raw scRNA-seq data from failing human hearts with dilated and hypertrophic cardiomyopathy. During data preprocessing, probe IDs were mapped to gene symbols. For genes associated with multiple probes, the probe with the highest average expression was retained as the representative. Detailed information on sample composition, disease etiology, sequencing platforms, sample sizes, and direct accession links for each dataset, is provided in the Supplementary Methods.

The Fragments Per Kilobase of transcript per million mapped reads (FPKM) values from the GSE135055 and GSE235897 datasets were log₂(X + 1) transformed, and all subsequent analyses were performed using R version 4.4.2. Differentially expressed genes (DEGs) between groups were identified using the “limma” package (version 3.62.2), with significance thresholds set at p < 0.05 and |log2FC| > 0.585 (corresponding to approximately a 1.5-fold change). The results were visualized using volcano plots and hierarchical clustering heatmaps, which were generated with “ggplot2” package (version 3.5.2) and “pheatmap” package (version 1.0.12).

To explore the potential biological functions and pathways implicated by the differentially expressed genes identified from the GSE135055 dataset, Kyoto encyclopedia of genes and genomes (KEGG) and Gene ontology (GO) enrichment analyses were performed on the DEGs using the “clusterProfiler” package (version 4.14.6), based on the “org.Hs.eg.db annotation” package (version 3.20.0), with significance defined as p < 0.05. The enrichment results were visualized using the “enrichplot” package (version 1.26.6) and the “ggplot2” package (version 3.5.2). To further evaluate innate immune-related transcriptional alterations at the pathway level, GSEA was performed on the GSE135055 dataset to assess pathway-level alterations between heart failure and healthy control samples. A pre-ranked gene list was generated based on gene expression between the two groups. Gene sets were obtained from the Molecular signatures database (MSigDB), including REACTOME_Innate_Immune_System from the Reactome pathway collection (C2: CP: Reactome) and GOBP_Innate_Immune_Response from the Gene Ontology Biological Process collection (C5: GO: BP). GSEA was conducted using the GSEA (version 4.4.0) following the pre-standard analysis workflow (8). Statistical significance was evaluated using normalized enrichment score (NES), along with nominal p values, false discovery rate (FDR) q values, and family-wise error rate (FWER) p values, as provided by the GSEA software.

To characterize the immune cellular landscape and coordinated immune interactions associated with heart failure, we identified immune cell subpopulations by mapping sequencing data to a validated immune cell signature matrix (LM22) representing 22 distinct immune cell types. The CIBERSORT algorithm was employed to quantify the relative abundance of these immune subpopulations and to compare compositional differences among samples (9). The resulting CIBERSORT-derived cell proportion matrix was extracted, converted to numeric format, and immune cell types exhibiting zero variance across samples were excluded to avoid spurious correlations. Spearman correlation coefficients were then calculated across immune cell types to characterize immune cell-cell association patterns. For visualization, the correlation matrix was rendered as an upper-triangular bubble plot (with the diagonal set to 1), with circle size proportional to p value and color indicating correlation direction; the correlation matrix was also exported for downstream analyses. To quantify coordinated shifts between immune populations, we additionally computed pairwise immune cell ratios from the CIBERSORT-derived fractions. For each sample and any two immune cell types “A” and “B”, the “A”/”B” ratio was defined as: (f_A + ϵ)/(f_B + ϵ), where f denotes the estimated cell fraction and ϵ = 1 × 10^-6 is a small pseudocount added to prevent division-by-zero artifacts when a cell type fraction approached zero. For visualization and group comparisons, ratios were log₂-transformed as log₂[(f_A + ϵ)/(f_B + ϵ)]. To further explore the correlations among immune functions, we selected ten representative innate immune-related subfunctions. Detailed descriptions of the selected pathways and the gene set-based scoring method are provided in the Supplementary Methods. Briefly, a Spearman correlation analysis was conducted to construct an immune functional interaction network, where correlations with a coefficient ≥ 0.6 were considered indicative of strong synergistic interactions. The network was visualized using Cytoscape (version 3.10.1), with node degree values used to highlight key functional interactions.

To identify core gene modules and key genes driving the transcriptional alterations observed in the GSE135055 dataset, WGCNA was performed using the “WGCNA” package (version 1.72-5) to identify gene co-expression modules associated with clinical traits. Genes were filtered based on expression variance, retaining the top 25% of the most variable genes. A soft-thresholding power (β) was selected to approximate scale-free topology (R² ≥ 0.90), and the resulting adjacency matrix was transformed into a topological overlap matrix (TOM) to measure network connectivity. Modules were identified via hierarchical clustering with dynamic tree cutting (minimum module size = 30), and similar modules were merged using a cut height of 0.25. Pearson correlation analysis was used to assess relationships between module eigengenes and clinical traits, and the two modules showing the strongest associations were selected for further analysis. All network construction and visualization were performed in R, with additional support from the “dplyr” package (version 1.1.4).

To determine which candidate gene plays more prominent and influential roles in heart failure, potential candidate genes were first screened by integrative analysis of heart failure-associated differentially expressed genes, high-weight genes derived from WGCNA, and trained immunity-related genes using a Venn diagram-based approach. A machine learning-based prioritization framework was then established in R to systematically rank these candidates together with identification of the core genes. The pre-selected genes from GSE135055 were used as input features, and gene expression data were merged with phenotypic group labels to construct a classification dataset. Stratified random sampling was applied to divide the samples into a training set (70%) and a test set (30%). Prior to model training, all gene expression features were centered and scaled based on the training data to ensure comparability across features. Six widely used machine learning models were trained using the “caret” package (version 7.0-1), including random forest (RF), support vector machine with a radial basis kernel (SVM), logistic regression (GLM), least absolute shrinkage and selection operator (LASSO) regression, k-nearest neighbors (KNN), and neural networks (NNET), with five-fold cross-validation employed to evaluate model performance. Model performance was initially evaluated using cross-validation-based classification metrics, and model generalizability was further assessed on the independent test set using the “DALEX” package (version 2.4.3). Predicted probabilities for the positive class were used to calculate probability residuals and root mean square error (RMSE) to characterize prediction uncertainty and calibration. Partial dependence plots were generated to visualize the marginal effect of individual genes on model predictions, while reverse cumulative distribution plots and boxplots of residuals were used to evaluate predictive accuracy. In addition, permutation-based variable importance analysis was performed to quantify and visualize the contribution of each gene across different models. To further screen candidate genes, we incorporated the mean variable-importance of different genes derived from machine learning models into principal component analysis (PCA). Specifically, for each candidate gene, a variable-importance matrix obtained from six models (RF, SVM, GLM, LASSO, KNN, and NNET). The importance values were z-score standardized across genes within each model (centered and scaled), and subsequently re-weighted to integrate model contributions based on their relative ranking derived from the preceding machine learning evaluation. A 3D PCA plot was generated based on these values, and genes that are farthest along each of the first three principal components (PC1, PC2, and PC3) were selected as potential validation targets. All figures were generated and exported using the “ggplot2” and “fs” packages (version 1.6.6).

To evaluate the classification performance and cross-dataset robustness of the selected genes, receiver operating characteristic (ROC) analysis was performed using the “pROC” package (version 1.18.5). ROC curves were generated in 4 independent heart failure cohorts for each gene, including GSE203160, GSE116250, GSE198945, and GSE48166. For each dataset, a gene-signature score was calculated by aggregating the expression values of the selected genes, and ROC curves were constructed based on the signature scores and corresponding group labels. The area under the curve (AUC) was calculated for each dataset and visualized in a combined ROC plot. Additionally, meta-analysis was performed using the “meta” package (version 8.1-0) to synthesize gene level standardized mean differences (SMDs) across independent datasets. Forest plots were generated to visualize dataset-specific and pooled effect sizes. Combined effect estimates from random-effects models, together with heterogeneity statistics, were used to assess the consistency and stability of gene expression patterns across datasets.

Single-gene GSEA was performed using a single-gene-based, pre-ranked approach to characterize biological pathways associated with the expression level of the gene of interest. Briefly, genes were ranked according to their Spearman correlation coefficients with the target gene, calculated across samples, and ranked in descending order. To avoid undefined values, minimal expression values were truncated prior to log₂ transformation. The resulting ranked gene list was used as input for single-gene GSEA analysis. GO gene sets were obtained from the Molecular signatures database (MSigDB) collection (c5.go.v2024.1.Hs.symbols.gmt) and analyzed using the “clusterProfiler” package (version 4.14.6). All GO terms were initially retained (pvalueCutoff = 1), and significantly enriched terms were defined based on a false discovery rate (FDR) q value < 0.25. For visualization, both enrichment curves and ridge plots were generated using the “enrichplot” package (version 1.26.6) to highlight the top enriched GO terms.

For single-cell transcriptomic analyses, we accessed the SCP1303 project from the Broad Institute Single Cell Portal. which contains raw scRNA-seq data derived from failing human hearts diagnosed with dilated cardiomyopathy (DCM) and hypertrophic cardiomyopathy (HCM), along with non-failing donor controls. The dataset includes gene expression profiles generated using droplet-based scRNA-seq technology, with detailed cell-type annotations provided by the original study. Raw scRNA-seq data were processed in Python (Scanpy version 1.9.3). Cells with low quality were excluded based on standard quality control criteria, including abnormal total UMI counts, gene detection rates, and elevated mitochondrial gene expression. Gene expression matrices were normalized and log-transformed where appropriate, and highly variable genes were identified for downstream dimensionality reduction. PCA was performed, followed by construction of a neighborhood graph and nonlinear dimensionality reduction using Uniform manifold approximation and projection (UMAP) for visualization. Based on curated cell-type annotations from the original dataset, macrophage populations were extracted for focused downstream analyses. To improve computational efficiency and ensure compatibility with trajectory inference, the data were further streamlined by retaining essential expression layers and metadata, generating an intermediate AnnData object for subsequent analysis. The data were then converted in R using zellkonverter (version 1.10.0) into a CellDataSet (CDS) compatible with monocle3 (version 1.3.1). Trajectory inference was conducted using the standard monocle3 workflow, including preprocessing, dimensionality reduction, cell clustering, principal graph learning, and pseudotime ordering. The root cell was manually specified to represent macrophages at the inferred early stage of the trajectory. To characterize transcriptional dynamics along pseudotime, MTURN expression was visualized across the inferred trajectory. Cells were subsequently stratified based on MTURN expression levels, with the top 40% classified as MTURN-high and the bottom 40% as MTURN-low, while intermediate cells were excluded to enhance group separation. Differential gene expression analyses were performed between MTURN-high and MTURN-low macrophages separately within DCM and HCM samples, enabling disease-specific characterization of MTURN-associated transcriptional programs.

Pairwise correlations were computed to assess associations between MTURN expression and heart-failure marker genes using normalized bulk RNA-seq expression matrices derived from human cardiac tissue samples in our previous dataset (GSE165303). Expression tables were processed in R using the “data.table” package (version 1.17.8), and correlation analyses were performed across all samples. Two-sided Spearman correlation coefficients and corresponding nominal p values were calculated for each marker, with p < 0.05 considered statistically significant. Results were visualized using scatter plots generated with the “ggplot2” package.

Based on published methods (10–12), we established a macrophage trained immunity model using THP-1-derived macrophages, complementing the Al(OH)₃/PHAD/mannan signature used in the transcriptomic discovery with an independent β-glucan-based approach. Briefly, THP-1 cells were differentiated into macrophages with phorbol 12-myristate 13-acetate (PMA; 100 nM, 24h), then “trained” with β-glucan (5 µg/mL, 24h). After a 24h rest in complete medium, cells were restimulated with lipopolysaccharide (LPS; 10 ng/mL, 24h); the resulting supernatant was collected as conditioned medium (CM) and applied to AC16 human cardiomyocytes for 24h. Controls used CM was generated from non-trained THP-1-derived macrophages that were restimulated with LPS (10 ng/mL, 24h) and applied to AC16 cells for 24h. THP-1 cells were maintained in RPMI-1640 supplemented with 10% FBS, 1% penicillin-streptomycin (P/S), and 0.05 mM β-mercaptoethanol; THP-1-derived macrophages were cultured in DMEM with 10% FBS and 1% P/S; and AC16 cells were maintained in DMEM/F12 (1:1) supplemented with 10% FBS and 1% P/S. All cells were cultured at 37 °C in a humidified atmosphere containing 5% CO₂.

Total RNA was extracted from THP-1-derived macrophages and AC16 cardiomyocytes using the RNeasy Kit (QIAGEN, #217004) in accordance with the manufacturer’s instructions. First-strand cDNA was generated using SuperScript II Reverse Transcriptase (Invitrogen, #18080044). Then qPCR was performed with AZURE CIELO3 system using SYBR Green Master Mix (Alkali Scientific, #QS1001). The mRNA expression levels were normalized to GAPDH and determined using the 2^-ΔΔCt method. All primer sequences for RT-qPCR are listed in Supplementary Table 1.

All statistical analyses were conducted using GraphPad Prism (version 8.02). Results are expressed as the mean ± standard deviation (SD). Comparisons between two groups were performed using a two-tailed unpaired Student’s T-test. For comparisons involving multiple groups, one-way analysis of variance (ANOVA) followed by Tukey’s post hoc test was applied. A p-value less than 0.05 was considered statistically significant (* denotes p < 0.05; ** denotes p < 0.01).

To investigate the transcriptomic alterations associated with heart failure, we first analyzed bulk RNA-seq data from the GSE135055 dataset, including 21 failing hearts and 9 healthy heart tissue samples. Compared with healthy controls, 336 genes were significantly upregulated in the cardiac tissue of heart failure patients, including NPPA and NPPB, which encode atrial and B-type natriuretic peptides and are well-established markers of cardiac stress and heart failure severity (13), as well as HBA1, HBA2, and HBB, which encode hemoglobin subunits and may reflect hypoxia-related or erythrocyte-associated transcriptional alterations in failing hearts (14). Conversely, 146 genes were significantly downregulated, including CA14, which encodes carbonic anhydrase 14 involved in pH regulation (15); ADAM11, a member of the disintegrin and metalloproteinase family implicated in cell-cell interactions (16); CPNE5, a calcium-dependent phospholipid-binding protein that we have previously shown to play a critical role in maintaining aortic integrity during early sepsis (17); and LSAMP, a cell-adhesion molecule linked to tissue structural maintenance (18). Collectively, these suppressed transcripts consistently point to a coordinated attenuation of signaling and structural integrity programs in failing myocardium. (Figures 1A, B, Supplementary Table 2). To investigate the biological relevance of these DEGs, we performed functional enrichment analyses. KEGG pathway analysis revealed significantly enriched pathways, among which cytoskeleton in muscle cells and the PI3K signaling pathway ranked among the top hits (Figure 1C). These genes can be broadly categorized into heart failure-associated stress and remodeling markers (ANKRD1 and FHL1) (19, 20), core sarcomeric and contractile components (MYOM1, MYL3, TNNI1, and TPM3) (21–23), cytoskeleton remodeling and mechanotransduction regulators (XIRP1, XIRP2, and DIAPH1) (24, 25), metabolic enzymes (CKM and ENO2) (26), and extracellular matrix-related genes with immune-modulatory potential (BGN, VCAN, FN1, THBS1, and multiple collagen family members) (27, 28). Such an enrichment pattern is reasonable, as the GSE135055 dataset contains a substantial number of heart failure-associated structural and contractile genes that are well represented in pathways related to muscle cytoskeleton organization and PI3K signaling. To identify immune-specific transcriptional alterations beyond tissue structure-dominated pathways, we next focused on three innate immune-related pathways including malaria, leishmaniasis, and the complement/coagulation cascades. Our analysis showed that these three pathways were enriched for genes involved in complement activation (C3, CFD, CFH/CFHR1, and MASP1) (29), coagulation-inflammation crosstalk (F2R/F2RL2, THBD, PROS1, SERPINE1/2, and VTN) (30), Fc receptor-mediated immune responses (FCGR3A and FCGR3B) (31), as well as endothelial adhesion molecules and key signaling regulators such as VCAM1, JAK2, NFKBIA, and TGFB2 (32, 33). Collectively, these analysis results suggest the presence of a pronounced innate immune activation signature in heart failure. Interestingly, GO enrichment analysis further identified pathways related to tissue remodeling and repair, including collagen fibril organization, collagen-containing extracellular matrix, as well as extracellular matrix structural constituent (Figure 1D). Such pathways likely reflect the complex structural and signaling adaptations occurring in the failing myocardium and provide a contextual framework within which immune-related processes may operate. Furthermore, we performed targeted gene set enrichment analysis using two innate immune-related gene sets: REACTOME Innate Immune System (NES = -1.12, FDR q-value = 0.031, and FWER p-Value = 0.031) and GOBP Innate Immune Response (NES = -1.13, FDR q-value = 0.038, and FWER p-Value = 0.038) (Figures 1E, F). Both gene sets were significantly enriched in heart failure samples. Together, these findings indicate that dysregulation of innate immune-related functions may contribute to the pathogenesis and progression of heart failure.

Next, to investigate the relationship between HF-associated transcriptional alterations and innate immune activation, we quantified immune cell infiltration in the left ventricular cohort (GSE135055) (Figure 2A), and found that HF samples exhibited a distinct immune remodeling compared with healthy controls (Figure 2B). Specifically, the proportions of M2 macrophages and follicular helper T cells were markedly decreased in HF, whereas CD8⁺ T cells were significantly increased. Of note, a similar reduction in M2-like macrophage signatures has been reported in several independent transcriptomic analyses of failing human myocardium (34–37). Importantly, the deconvolution-inferred “M0/M1/M2” macrophage fractions were not interpreted as definitive evidence of canonical polarization programs, but could be as phenotype-associated transcriptional signatures reflecting immune reprogramming at the tissue level. To figure out the potential relationship underlying these changes across innate immune subfunctions, we constructed an immune functional interaction network using the immune infiltration profiles from GSE135055, and revealed that there was a high interconnectivity between trained immunity- and macrophage polarization-related modules (Figure 2D), suggesting that these processes may represent central axes of immunological remodeling in HF. In addition, we examined immune cell-cell correlation patterns within HF samples, which revealed coordinated shifts across myeloid and lymphoid compartments rather than isolated changes (Figure 2C). In particular, monocytes were positively correlated with activated CD4⁺ memory T cells, and regulatory T cells (Tregs) were positively associated with M0 macrophages signatures. In contrast, Tregs showed a negative correlation with activated NK cells, and M1 macrophages were inversely correlated with memory B cells. It is important to note here, M0 macrophages in deconvolution frameworks represent transcriptionally inferred, relatively non-polarized myeloid states rather than functionally inert cells; therefore, these correlations could be interpreted as reflecting immune microenvironmental balance and myeloid differentiation potential, rather than direct cell-cell interactions (38, 39). Within this framework, the positive association between Tregs and M0-like macrophage signatures may indicate a relatively immunoregulatory milieu in which macrophage polarization is restrained or delayed, consistent with a transitional myeloid reprogramming state. Accordingly, such coordinated immune remodeling is compatible with sustained innate immune activation and reprogramming, a hallmark of trained immunity described in chronic inflammatory settings (7). Together, these results suggest the presence of trained immunity-like immune reprogramming in HF. Moreover, analysis of immune cell ratios revealed that heart failure patients had significantly elevated ratios of monocytes to M0 and M2 macrophages, along with an increased ratio of M1/M2 compared to healthy controls. (Figures 2E-G). These ratio shifts could be interpreted as functional readouts, consistent with trained immunity-associated macrophage reprogramming rather than macrophage polarization. Therefore, dynamic remodeling of macrophage-associated signatures and trained immunity-related immune reprogramming may contribute to HF development and progression.

At present, comprehensive studies summarizing gene expression changes associated with trained immunity, particularly in macrophages, are still lacking. Hence, we selected GSE235897 database in which human monocytes were differentiated into macrophages and subsequently trained with a combination of Al(OH)₃, PHAD, and mannan (Figure 3A), as this model provides a basis for investigating the transcriptional changes associated with macrophage-trained immunity. With differential expression analysis, we identified 147 significantly upregulated and 171 significantly downregulated genes in macrophages following training (Figures 3B, C, Supplementary Table 3). These DEGs could be compiled into a gene set representing the transcriptional signature of trained immunity.

Given that the DEGs derived from the heart-failure dataset (GSE135055) were extensive, we next applied WGCNA to identify the core genes associated with the HF phenotype. After confirming the absence of outlier samples through hierarchical clustering, we selected a soft-thresholding power of 9 to satisfy the scale-free topology criterion (Figures 3D-G). Subsequently, multiple gene modules were identified using dynamic tree cutting. Among them, the black module (MEblack, r = 0.92, p < 1e-12) was positively correlated with heart failure, while the green module (MEgreen, r = - 0.74, p < 1e-6) was negatively correlated. Together, these two modules encompassed a total of 697 genes (Figure 3H, Supplementary Table 4). Finally, Venn integration of WGCNA module genes, HF DEGs, and the trained-immunity signature identified seven candidate hub genes—PTGDS, TNFRSF12A, MLLT11, RND3, MTURN, AQP3, and COL23A1—that may represent key molecular links between macrophage trained immunity and HF pathogenesis (Figure 3I).

Building on the integrative identification of seven candidate genes with distinct expression patterns in heart failure, we next applied multiple machine learning algorithms to assess their diagnostic significance (Figure 4A). Specifically, generalized linear model (GLM), support vector machine (SVM), random forest (RF), and neural network-based models were implemented in parallel to capture complementary predictive characteristics across algorithms. Residual boxplots showed that GLM and SVM had more tightly clustered error distributions, with reduced variance and fewer extreme residuals, indicating greater model stability (Figure 4B). Meanwhile, reverse cumulative distribution plots of residuals revealed that GLM and NNET yielded the lowest prediction errors across the majority of samples, suggesting superior predictive performance (Figure 4C). Furthermore, feature importance rankings further demonstrated that GLM, SVM, and RF maintained stable performance under variable perturbations, indicating enhanced model robustness and reduced sensitivity to noise or sampling variation (Figure 4D). To more accurately identify top-performing genes, we conducted a weighted PCA integrating outputs from multiple models to minimize bias introduced by any single algorithm (Figure 4E). The PCA results revealed substantial variation in the spatial distribution of genes, with PTGDS, MTURN, and TNFRSF12A located furthest from the centroid in three-dimensional space, reflecting strong and consistent contributions across models. Therefore, these genes exhibit strong discriminative capacity, with MTURN in particular highlighting a close association with macrophage-mediated trained immunity-related transcriptional programs in heart failure.

To validate the machine learning-derived performance of the candidate genes, especially MTURN, we evaluated their expression and diagnostic value in four independent HF datasets. Importantly, ROC curve analysis demonstrated that MTURN consistently achieved high AUC values across all datasets (maximum AUC = 0.992, outperforming PTGDS and TNFRSF12A (Figures 5A-C, G). Furthermore, after normalizing expression data, we assessed the differential expression patterns of these genes and found that MTURN remained consistently upregulated in HF samples (Figures 5D-F, H). A subsequent meta-analysis of the combined datasets revealed that MTURN had the highest SMD of 1.57 (95% CI: 0.70-2.45), indicating robust and consistent upregulation in HF samples (Figure 5I). In addition, single-gene GSEA revealed that MTURN-associated genes were significantly enriched in immune activation pathways, including leukocyte-mediated immunity, regulation of adhesion molecules, and immune response remodeling (Figures 5J, K). Next, to explore the relevance of the trained immunity transcriptional signature to heart failure, we performed correlation analysis between MTURN expression, and the dysregulated genes identified in the trained immunity model. Notably, MTURN showed strong positive correlations with multiple pro-inflammatory and immune activation-related genes (Supplementary Table 5), including CXCL8, CEACAM3, CYGB, and COL23A1 (40–42). In contrast, MTURN expression was negatively correlated with genes involved in immune regulation and metabolic homeostasis, such as IL1RN, CD36, and ALDH1A2 (43, 44) (Supplementary Figure S1). Collectively, these findings suggest that the elevated expression of MTURN in HF is a consistent and widespread phenomenon, highlighting its strong association with macrophage-mediated immune regulation and underscoring its potential as a promising molecular target for future investigation.

To address the potential role of MTURN in macrophages during the development of heart failure, we selected a scRNA-seq dataset (SCP1303 project from failing human hearts with DCM and HCM) for analysis and UMAP visualization revealed the apparent enrichment in adipocytes, endothelial cells, and macrophages. Further dot-plot analysis results showed that higher expression levels of MTURN were exhibited in adipocytes, followed by endothelial cells, and macrophages, within human hearts (Figure 6A, Supplementary Table 6). More interestingly, the expression levels of MTURN were increased by 2.1-fold in cardiac macrophages only, but not in either adipocytes or endothelial cells within failing hearts, compared to healthy donors (Figure 6B). To explore the dynamic regulation of MTURN during macrophage state transitions, we performed pseudotime trajectory analysis. Cells were ordered along an inferred temporal axis and revealed progressive increase in MTURN expression during macrophage differentiation and activation, followed by stabilization at later pseudotime stages (Figures 6C-E). Next, we separated the MTURN-expressing macrophages into high and low expression groups and performed parallel differential expression analysis in HCM and DCM samples (Figure 6F, Supplementary Tables 7, 8). This analysis identified several immunity related genes, such as PIEZO1, IGSF6, and EIF3, that were consistently upregulated in MTURN-high macrophages across both disease conditions. Given the prominent enrichment of PIEZO1 in MTURN-high macrophages, we further investigated the dynamic behavior of PIEZO1 at the single-cell level. Pseudotime trajectory analysis revealed that PIEZO1 exhibited a clear pseudotime-dependent expression pattern, characterized by early induction during macrophage state transitions, preceding the upregulation of MTURN, followed by a gradual decline at later stages (Supplementary Figure S2A). Although PIEZO1 and MTURN displayed partially distinct temporal dynamics, both genes showed coordinated downregulation toward the terminal pseudotime states, suggesting their involvement in stage-dependent macrophage activation programs rather than terminal differentiation. In addition, macrophages were separated into PIEZO1-high and PIEZO1-low subsets, and comparative analysis demonstrated that MTURN expression was significantly higher in PIEZO1-high group than in PIEZO1-low group (Supplementary Figure S2B). Together, these findings suggest a close functional association between PIEZO1 and MTURN during macrophage activation and immune training in heart failure.

To validate our bioinformatics findings, we first revisited our previous bulk RNA-seq dataset (GSE165303) comprising human cardiac samples from healthy controls and patients with DCM (Figure 7A). Consistent with our computational analysis, we observed that the expressions of PIEZO1 and MTURN were significantly increased in the DCM group (Figures 7B, C). We next examined whether MTURN expression was associated with established heart failure markers at the tissue level. Spearman correlation analysis revealed that MTURN expression was positively correlated with the heart-failure biomarkers NPPA (r = 0.55; p = 2.33 × 10^-9), NPPB (r = 0.32; p = 1.04 × 10^-3), and TNNI3 (r = 0.32; p = 1.27 × 10^-3) (Figures 7D-F). In contrast, direct correlation analysis between PIEZO1 and MTURN expression in bulk cardiac tissue did not reveal a significant linear association (Spearman r = 0.07, p = 0.499; as shown in Supplementary Figure S3). These results suggest that, although PIEZO1 and MTURN are concurrently upregulated in DCM, their expression levels are not tightly coupled at the whole-tissue level. To further investigate their regulation in a cell-type-relevant context, we established an in vitro trained immunity model using human THP-1-derived macrophages and AC16 cardiomyocytes (Figure 8A). In trained immunity THP-1-derived macrophages, the mRNA levels of glycolysis-related genes (LDHA, SLC2A1, and HK2) were significantly increased (Figure 8B). Meanwhile, the expression of inflammation-related genes (TNF, IL-6, and IL-1β) was also markedly upregulated (Figure 8C), indicating successful induction of trained immunity. Notably, trained immunity further upregulated MTURN and PIEZO1 expression in THP-1-derived macrophages (Figure 8D), indicating that these two genes are co-responsive under macrophage trained immunity conditions. Finally, conditioned medium collected from these macrophages markedly increased the expression of heart-failure markers (NPPA, NPPB, and TNNI3) in AC16 cardiomyocytes (Figure 8E). Collectively, these results validate our bioinformatic analysis and suggest that trained immunity-associated macrophage activation is accompanied by coordinated upregulation of MTURN and PIEZO1. While their expression is not directly correlated at the bulk tissue level, their concordant induction in macrophages supports a potential functional convergence in the context of heart failure-associated immune reprogramming.