The data for this study were sourced from publicly available datasets. Bulk RNA-seq analysis incorporated four independent cohorts from the GEO database: GSE57338 (136 controls vs. 177 heart failure samples), GSE5406 (14 vs. 194), GSE141910 (166 vs. 200), GSE116250 (14 vs. 50) and GSE79962 (11 vs. 20) (10–14). For single-cell resolution, three scRNA-seq datasets were analyzed: GSE145154 (105,533 cells from 6 samples), GSE183852 (49,723 cells from 8 samples, including 2 controls and 5 HF cases), and GSE121893 (4,933 cells from 12 samples, including 2 controls and 10 HF cases) (15–17). Genes associated with histidine metabolism were identified using the corresponding pathway (map00340) from the KEGG database, with the full list provided in Supplementary Table 1.

For bulk RNA-seq data, batch effects across different datasets were removed using the sva R package prior to machine learning analysis (18). Normalised expression matrices were employed, comprising FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values for sequencing data and log2-transformed values for microarray data. For single-cell data, analysis was performed using the Seurat package (v5.30) (19). During initial quality control, we excluded cells with >5,000 or <200 detected genes, >10% mitochondrial gene content, and >5% ribosomal or haemoglobin gene content. The remaining cells were subsequently normalised using the SCTransform method. Principal Component Analysis (PCA) was performed based on the top 3,000 highly variable genes. Clustering was conducted using a resolution of 0.8, and cell types were manually annotated based on canonical marker genes. The Harmony algorithm was employed to integrate samples from different scRNA-seq datasets (20).

To evaluate pathway activity in individual cells, six different gene set enrichment approaches were applied: AddModuleScore, AUCell, UCell, singscore, ssGSEA, and GSVA. Using the Seurat platform, AddModuleScore computes a module score by contrasting the mean expression of a gene set against a randomly selected control set (with 100 genes). The AUCell tool (version 1.30.1) determines an Area Under the Curve value from the cumulative distribution of expression ranks within the gene set, which represents the global measure of pathway activity. In addition, it is specifically optimized for the high sparsity and dropout effect of single-cell RNA-seq data. Therefore, we use the enrichment scores calculated by this algorithm to perform subsequent analyses. UCell (version 2.12.0) follows a comparable logic but employs the Mann–Whitney U statistic to improve computational speed. With the singscore package (version 1.28.0), a score is produced that incorporates both gene ranking and expression level. It is worth noting that, due to the characteristic sparsity of single-cell RNA-seq data, the top highly variable genes were supplied as input for the ssGSEA and GSVA algorithms (as implemented in GSVA version 2.2.0) (21). For all subsequent analyses, the histidine metabolism score generated by AUCell was adopted as the main measure, while results from the remaining five techniques served for additional verification. The correlation analysis between AUCell scores and scores of other algorithms is shown in the Supplementary Table 2.

To identify functionally enriched categories, Over-Representation Analysis (ORA) was carried out with reference to the Gene Ontology (GO) database. The clusterProfiler R package (version 4.16.0) was employed to detect statistically over-represented terms associated with biological processes, molecular functions, and cellular components (22). Visualisation of the significant enrichment results was accomplished using ggplot2 (version 3.5.02). For the bulk RNA-seq data, differential expression analysis was performed with the limma package (version 3.62.1) (23). Genes meeting the thresholds of an absolute log2 fold change (|log2FC|) greater than 0.585 (corresponding to a fold change >1.5) and a p-value below 0.05 in comparisons between HF and Control samples were classified as differentially expressed.

After batch effect correction, the integrated dataset was divided for model construction. Given its large sample size, the GSE141910 cohort was assigned to the training set. The remaining three cohorts (GSE79962, GSE5406, and GSE57338) were set aside as external validation sets to evaluate the generalizability of the models. A diverse set of 15 machine learning algorithms was utilized, comprising: Neural Network, Logistic Regression, Linear Discriminant Analysis, Quadratic Discriminant Analysis, k-Nearest Neighbors, Decision Tree, Random Forest, Extreme Gradient Boosting, Ridge Regression, LASSO, Elastic Net, Support Vector Machine, Gradient Boosting Machine, Stepwise Logistic Regression, and Naïve Bayes. Following systematic optimization of hyperparameters, a total of 117 unique model variants were developed and assessed. Model performance was measured using four criteria: F1 Score, Area Under the Curve, Accuracy, and Recall. The resulting performance outcomes were visualized with the ComplexHeatmap package in R (version 2.24.1) (24).

Global cell-cell communication networks were inferred using the CellChat R package (version 2.1.2) (25). For analyses within macrophage and fibroblast subpopulations, CellChat objects were created separately for the “high histidine” and “low histidine” groups (defined by AUCell scores). CellChat infers communication probability based on a curated ligand-receptor (L-R) interaction database and the average gene expression within each cell cluster. The overall interaction strength and count were quantified and compared between the two metabolic states using the compareInteractions function. Differential signalling pathways and key “hub” populations (identified via netAnalysis_computeCentrality) were analysed to characterise the communication patterns specific to the high-histidine microenvironment. MultiNicheNet Analysis: To investigate specific, target-gene-level interactions between FRZB⁺ MSCs and CRYM⁺ cardiomyocytes, the MultiNicheNet R package (v2.0.1) was employed (26). This algorithm predicts L-R interactions correlated with downstream gene expression changes in the receiver cell. The analysis was conducted bidirectionally (MSC-to-cardiomyocytes and cardiomyocytes-to-MSC) comparing the HF versus Control conditions. MultiNicheNet prioritises active L-R pairs by assessing their correlation with the expression of predicted downstream target genes. Subsequently, the target genes regulated by these prioritised L-R pathways were analysed against the GO database to infer the biological consequences of the communication.

To quantify the abundance of key cell populations within the bulk tissue of the GSE141910 cohort, we employed the IOBR R package (version 0.99.8) (27). We specifically selected two complementary deconvolution algorithms encapsulated within IOBR to robustly assess the stromal origin of our candidate genes (e.g., FRZB): MCPcounter: This algorithm was chosen for its ability to provide abundance scores for specific cell types based on curated marker gene sets. Our analysis focused on the “Fibroblasts” score (Fibroblasts_MCPcounter) to test the direct correlation between our candidate genes and a specific mesenchymal cell population. ESTIMATE: In parallel, the ESTIMATE algorithm was used to calculate a global StromalScore. This score quantifies the overall presence of the entire stromal compartment (including, but not limited to, fibroblasts) within the tissue microenvironment. The rationale for this dual-method approach was to establish a more robust validation: by demonstrating that our candidate gene expression correlated positively with both a specific cell-type score (MCPcounter) and a global compartment score (ESTIMATE), we could more confidently confirm their dominant expression within the cardiac stromal lineage.

All statistical analyses were conducted using the R (version 4.5.1) and Python (version 3.10) programming languages. Data comparisons across groups were performed using either an independent-samples t-test or a Kruskal–Wallis test. A p-value of less than 0.05 was regarded as statistically significant.

This study examined remodeling of the cardiac microenvironment in heart failure (HF) through analysis of scRNA-seq dataset GSE145154. After implementing stringent quality control steps—such as filtering out low-quality cells (e.g., those exhibiting elevated expression of mitochondrial, erythrocyte, or ribosomal genes) and removing duplicates (Supplementary Figures S1A–C)—and integrating data from multiple samples (Figure 1A), all cells were grouped into 19 clusters via unsupervised clustering (Figure 1B) and assigned to 10 major cell types using established marker genes (Figure 1C). Comparative assessment of cellular composition uncovered notable shifts between HF and control samples (Figure 1D). Control samples were characterized by a higher abundance of stromal cells, including fibroblasts and endothelial cells, while HF specimens displayed pronounced accumulation of immune cells such as T cells and B cells, indicating substantial immune infiltration and tissue remodeling in HF. To confirm annotation accuracy, we visualized canonical marker gene expression on UMAP projections, which confirmed cluster-specific expression patterns (Figure 1E). Differentially expressed genes were also detected across individual cell types (Figure 1F).

To quantitatively evaluate compositional variation and intergroup heterogeneity, we applied three computational approaches. MiloR analysis (Figure 1G) identified significant enrichment of T cells, NK cells, pericytes, B cells, and SMCs in HF tissues. Conversely, myeloid cells, fibroblasts, endothelial cells, and cardiomyocytes were broadly distributed across both conditions, implying considerable heterogeneity within these populations. Augur analysis (Figure 1H), which evaluates transcriptomic state changes, highlighted fibroblasts and myeloid cells as the populations with the largest transcriptional shift between groups. Finally, subpopulation preference (RO/E) analysis (Figure 1I) further reinforced that fibroblasts and myeloid cells were more characteristic of control tissues, whereas lymphocyte subsets were preferentially associated with HF. Together, these complementary analytical approaches consistently identified fibroblasts and myeloid cells as the most heterogeneous and differentially distributed populations within the failing cardiac microenvironment. Based on these findings, we focused subsequent investigations on these two cell types.

Based on the preceding findings, we next explored heterogeneity in specific biological functions within the heart failure (HF) microenvironment. Activity of the histidine metabolism pathway at single-cell resolution was evaluated by integrating six gene set enrichment algorithms—ssGSEA, GSVA, UCell, AUCell, singscore, and AddModuleScore. Among all cell types, cardiomyocytes, fibroblasts, mast cells, and myeloid cells displayed notably elevated histidine metabolism scores (Figure 2A). For a global characterization of the metabolic landscape, cells were stratified into high- and low-activity groups based on quartile values. This broad classification allowed for a comprehensive observation of the heterogeneous distribution pattern, where regions enriched for high-activity cells predominantly coincided with tissue zones associated with fibroblasts and myeloid cells (Figures 1I, 2B). To corroborate these findings at the tissue level, we analyzed four independent transcriptomic datasets (GSE5406, GSE57338, GSE141910, GSE116250). A complex dysregulated expression profile of histidine metabolism was detected in cardiac tissues from HF patients relative to controls (Figures 2D–G). Notably, this activation pattern was observed across HF patients with varying etiologies, including dilated cardiomyopathy (DCM) and ischemic cardiomyopathy (ICM), indicating that enhanced histidine metabolism represents a common biological hallmark of heart failure. Finally, we examined expression changes of key genes involved in this pathway. Rather than uniform upregulation, a dysregulated expression profile was observed (Figures 2H–K): several genes (e.g., ALDH2, HNMT, HDC) were significantly elevated in HF samples, whereas other critical genes (e.g., MAOA, CNDP1) showed marked downregulation.

To explore cellular heterogeneity in histidine metabolism during heart failure, we began by independently reclustering myeloid cells (Supplementary Figures S2A–C), classifying them into monocyte-macrophages (Mφ), neutrophils, and dendritic cells (DCs). Annotation accuracy was supported by subtype-specific marker gene expression, including LYVE1 in Mφ, FCER1A in DCs, and G0S2 in neutrophils (Supplementary Figure S2E). Among these, Mφ displayed the highest histidine metabolism activity (Supplementary Figure S2D), leading us to focus subsequent analysis on this population. Further subclustering of Mφ identified seven functionally distinct subtypes (Figure 3A): CCL3L3⁺, CXCL2⁺, GPNMB⁺, MYL2⁺, RNASE1⁺, FCN1⁺, and Classic macrophages, each showing clear transcriptional boundaries (Figure 3D) and defined by specific marker genes (Supplementary Figures S2F,G). We then applied CellChat to investigate communication patterns among macrophage subpopulations. The analysis revealed intensified interaction networks within “high histidine” macrophage clusters (Figures 3E–G), with CCL3L3⁺ and RNASE1⁺ subtypes serving as central signaling hubs. Differential pathway analysis further identified upregulation of MHC-II and LAIR1 signaling pathways in “high histidine” macrophages (Figure 3H). Notably, enhanced MHC-II signaling was observed across multiple subtypes within this group, especially when CCL3L3⁺ and Classic macrophages acted as signal senders (Figure 3I).

We next turned to fibroblasts to examine functional heterogeneity linked to histidine metabolism. After clustering at a stable resolution (0.1) (Figure 4A), five functionally distinct fibroblast subpopulations were identified (Figure 4B): Fib_ACTA2, Fib_FGF7, Fib_THBS4, Fib_THY1, and Fib_GPM6B. To assess the robustness of this classification, we first applied the ROGUE metric, which indicated high purity across all subtypes (ROGUE > 0.7), with Fib_GPM6B exceeding 0.9 (Figure 4C). We then performed differential gene expression analysis, which confirmed clear transcriptional distinctions among the five groups, aligning with their spatial segregation in UMAP visualization (Figure 4D). Functional characterization via Gene Ontology enrichment further highlighted specialized roles for each subtype (Figure 4E): Fib_FGF7 was associated with ribosome-related processes; Fib_ACTA2 was involved in extracellular matrix organization and cytoplasmic translation; Fib_THY1 was related to osteoblast development; Fib_GPM6B was connected to neuronal pathways; and Fib_THBS4 was enriched in insulin-like growth factor (IGF)-mediated signaling.

We then integrated these fibroblast subtypes with histidine metabolism activity profiles and analyzed their intercellular communication using CellChat. The results showed that fibroblast subpopulations with high histidine metabolism activity generally engaged in stronger cellular interactions (Figures 4F–I). Among them, the Fib_THY1 subtype emerged as the predominant signal sender within high-activity groups. Pathway-level examination identified several signaling routes—including those mediated by NCAM, TWEAK, SLIT, and ANNEXIN—that were markedly more active in high-histidine fibroblasts (Figure 4J). Focusing on interactions initiated by Fib_THY1, we observed a pronounced increase in signaling probability via the collagen-CD44 ligand–receptor pairs COL1A2-CD44 and COL1A1-CD44 between Fib_THY1 and myeloid cells in the high-histidine microenvironment (Figure 4K). Alterations in signaling patterns originating from other fibroblast subtypes are provided in Supplementary Figure S3.

To identify the most influential regulators with high precision, we implemented a more stringent, extreme-value filtering strategy for core target screening. Single-cell histidine metabolism pathway activity scores were determined (Figure 5A), enabling the classification of cells into three distinct categories: a high-activity group (top 10%, score > 90th percentile) to capture the most potent metabolic signatures, a low-activity group (absolute zero, score = 0) to serve as a clean baseline, and an intermediate group (Figure 5B). This rigorous contrast was designed to maximize the signal-to-noise ratio during the subsequent correlation and differential expression analyses. Subsequent Spearman correlation analysis between highly variable gene expression and metabolic activity scores revealed a sigmoidal relationship (Figure 5C), indicating a nonlinear association characterized by threshold behavior—gene expression initially increases gradually with rising pathway activity, then accelerates sharply beyond a critical point before plateauing under high-activity conditions. We applied a sequential filtering strategy to identify robust candidate genes: single-cell correlation analysis identified 1,671 genes positively correlated with metabolic scores (Cor > 0, P < 0.05); single-cell differential expression analysis revealed 1,969 genes significantly upregulated in high- versus low-activity groups (avg_log2FC > 0.585, pct.1 - pct.2 > 0, p_val_adj < 0.001; Figure 5E); conventional transcriptome validation using dataset GSE141910 identified 1,993 genes upregulated in HF versus controls (logFC > 0.585, p < 0.05; Figure 5D); and additional validation with dataset GSE57338 confirmed 150 differentially upregulated genes (logFC > 0.585, p < 0.05; Figure 5D). Intersection of these four gene sets yielded 23 core candidate genes (Figure 5F) that consistently demonstrated association with histidine metabolism activation in heart failure across both single-cell and bulk transcriptomic levels, establishing them as priority targets for subsequent investigation.

Building on the 23 candidate genes previously identified, we next sought to refine core targets and develop a diagnostic model. The four bulk transcriptomic datasets were initially harmonized through batch effect correction (Figures 6B,C). We then evaluated two gene sets for predictive performance: the full panel of 23 candidates, and a 9-gene subset refined by LASSO regression (using the Lambda.1se criterion) (Figure 6D). Both sets were used to train and optimize 15 machine learning algorithms—generating 117 model configurations—with performance assessed by AUC (Figure 6A), Recall, F1 Score, and Accuracy (Supplementary Figures S4A–C). Models incorporating all 23 genes exhibited strong diagnostic capability, reflected in high AUC values, reinforcing the relevance of histidine metabolism in heart failure and underscoring the collective importance of these genes. Interestingly, the 9-gene LASSO-derived subset did not surpass—and in some cases underperformed—the full gene set, particularly in non-linear classifiers such as Neural Networks and KNN. This implies that the original 23 genes form an information-rich, synergistic signature. Boruta analysis further validated this observation, confirming all 23 genes as “important” features (Figure 6H). As LASSO operates through linear constraints, it may have eliminated genes that contribute meaningfully in non-linear decision contexts. We therefore turned to machine learning-driven ranking to identify the most influential regulators among the 23 candidates. Two complementary ranking approaches were applied: at the transcriptomic level, we computed the correlation between gene expression and histidine metabolism activity scores (Figure 6F); separately, we extracted Variable Importance metrics derived from model outputs (Figure 6G). By intersecting genes ranked among the top five in both “correlation” and “model importance”, we identified two core target genes—FRZB and CRYM—following this multi-tiered validation.

To evaluate the clinical significance of the two core target genes (FRZB and CRYM), we examined their expression patterns and diagnostic performance in two independent bulk transcriptome cohorts (GSE57338 and GSE141910). In the GSE57338 dataset, expression levels of both FRZB and CRYM were significantly elevated in heart failure samples relative to controls (Figures 7A,B). Receiver operating characteristic (ROC) analysis demonstrated strong diagnostic ability, with area under the curve (AUC) values reaching 0.917 for FRZB and 0.871 for CRYM (Figure 7C). These results were consistently replicated in the GSE141910 cohort, where FRZB and CRYM again showed marked upregulation in HF specimens (Figures 7E,F) and achieved high diagnostic accuracy (AUC = 0.983 for FRZB; 0.928 for CRYM) (Figure 7G). Furthermore, we investigated the association between these genes and histidine metabolism pathway activity. Correlation analyses across both cohorts revealed two key findings: first, a strongly significant positive co-expression relationship between FRZB and CRYM (GSE57338: R = 0.54; GSE141910: R = 0.59; both P < 2.2 × 10⁻¹⁶), indicative of potential coregulation; and second, a consistent positive correlation between each gene and histidine metabolism scores (e.g., in GSE141910: FRZB R = 0.23, CRYM R = 0.34), reinforcing their roles as central drivers of histidine metabolic activation in HF. Notably, the current study focuses on transcriptome analysis at the organizational level. Further clinical validation in combination with the peripheral blood circulatory system (such as plasma or peripheral blood mononuclear cells) is required to evaluate its potential as a non - invasive biomarker.

To determine the cellular sources of FRZB and CRYM in the cardiac microenvironment, we mapped their expression patterns across three independent single-cell RNA sequencing datasets (GSE145154, GSE121893, GSE183852). The two genes displayed clearly distinct cellular expression profiles. CRYM showed highly specific localization, being predominantly and consistently expressed in cardiomyocytes in both the GSE145154 (Figures 8A,B) and GSE121893 (Figures 8C,D) datasets. This cardiomyocyte-specific expression was further supported by strong positive correlations between CRYM and canonical cardiomyocyte markers such as MYH6, MYH7, and TNNT2 (Supplementary Figure S5G). In contrast, FRZB was found to be primarily enriched within mesenchymal stromal lineages. In the GSE145154 dataset, it showed elevated expression in pericytes, fibroblasts, and smooth muscle cells (Figures 8A,B). This stromal preference was consistently observed in the GSE121893 (Figures 8C,D) and GSE183852 (Supplementary Figures S5A–D) datasets. To further confirm its stromal origin, we evaluated the relationship between FRZB expression and established stromal signatures, identifying strong positive correlations with both Fibroblasts_MCPcounter (cor = 0.42) and StromalScore_estimate (cor = 0.51) (Supplementary Figures S5E,F). In summary, the two core target genes display clear cellular lineage separation in the failing heart: CRYM is predominantly expressed by cardiomyocytes, while FRZB is mainly derived from mesenchymal stromal cells.

While our preliminary analysis identified Fib_THY1 as a key active fibroblast subpopulation, the molecular drivers connecting these stromal shifts to myocardial dysfunction remained elusive. Having identified FRZB (stromal-specific) and CRYM (cardiomyocyte-specific) as core metabolic regulators, we hypothesized that they mediate a critical stroma-myocardium communication axis that sustains HF progression. To investigate the functional specificity of cellular subpopulations expressing the core target genes in heart failure, we stratified mesenchymal cells (MSCs) and cardiomyocytes (cardiomyocytes) into four subsets according to FRZB and CRYM expression patterns: FRZB⁺ MSCs, FRZB⁻ MSCs, CRYM⁺ cardiomyocytes, and CRYM⁻ cardiomyocytes. UMAP projections (Figure 9A) and correlation analyses (Figure 9B) verified that these groups possess distinct molecular and spatial profiles. Tissue preference analysis (Figure 9C) further demonstrated that both FRZB⁺ MSCs and CRYM⁺ cardiomyocytes were preferentially localized in the HF microenvironment—in contrast to FRZB⁻ MSCs, which were more abundant in control tissues—with CRYM⁺ cardiomyocytes exhibiting the most prominent enrichment, indicating their potential role as central contributors to HF pathogenesis. We then utilized MultiNicheNet to examine communication between these enriched subpopulations, which revealed a strongly augmented bidirectional signaling network under HF conditions. In the stromal-to-myocardial direction (FRZB⁺ MSCs → CRYM⁺ cardiomyocytes), multiple ligand-receptor interactions—including FN1-ITGB1, FN1-SDC2, and TNFSF12-TNFRSF12A—were significantly elevated when FRZB⁺ MSCs served as signal sources (Figure 9D). Downstream target analysis (Supplementary Figure S6C) indicated that these signals primarily orchestrate metabolic rewiring in CRYM⁺ cardiomyocytes, affecting pathways such as nucleotide metabolism and glyoxylate/dicarboxylate metabolism. In the reciprocal myocardial-to-stromal direction (CRYM⁺ cardiomyocytes → FRZB⁺ MSCs), ligand-receptor axes such as LAMP4-ITGB1 and FGF7-NRP1 were also strengthened upon ligand emission from CRYM⁺ cardiomyocytes (Supplementary Figure S6A). Downstream enrichment (Supplementary Figure S6B) suggested that these cues from cardiomyocytes mainly provoke structural and pro-fibrotic phenotypic shifts in FRZB⁺ MSCs, including regulation of myocyte growth and actin filament-associated processes. Together, these observations support a computational model of a self-sustaining pathogenic cycle between FRZB⁺ stromal cells and CRYM⁺ cardiomyocytes in HF: stromal cells direct metabolic adaptation in cardiomyocytes, while cardiomyocytes reciprocally prime stromal cells for phenotypic transition, such as toward myofibroblast conversion. This mutually reinforcing dialogue appears to collectively propagate disease progression. Rather than a proven computational model, these inferred interactions suggest that stromal cells may contribute to metabolic adaptation.