Macrophage-enriched HMOX1 is an ischemia-reperfusion-responsive marker in the heart

Introduction

Reperfusion remains the cornerstone therapy for limiting myocardial necrosis after acute coronary occlusion (1,2). However, restoration of blood flow can also initiate myocardial ischemia-reperfusion (I/R) injury, a multifactorial process that contributes to infarct expansion, malignant arrhythmias, microvascular dysfunction, and adverse ventricular remodeling (3,4). Despite decades of investigation, clinically effective therapies that specifically target I/R injury are still lacking (5).

Acute I/R rapidly activates inflammatory and redox-sensitive signaling networks involving circulating leukocytes, resident immune cells, and stromal populations (6,7). Within this response, macrophages are notable for their context-dependent functional plasticity (8-10). Early pro-inflammatory macrophage programs may aggravate tissue injury, whereas subsequent reparative phenotypes facilitate inflammation resolution and myocardial healing (11). Defining the molecular programs that accompany these transitions may inform new therapeutic strategies (12,13).

Transcriptomic profiling offers an unbiased view of gene programs induced by I/R; however, conventional differential-expression analyses can be underpowered in small cohorts and may miss coordinated, module-level regulation. Network-based approaches such as weighted gene co-expression network analysis (WGCNA) can identify co-regulated gene modules and nominate hub genes more robustly (14). Meanwhile, feature-selection and interpretable modeling can further prioritize compact, informative gene sets and connect them to biological processes (15,16). A remaining challenge of bulk transcriptomics is that it averages signals across heterogeneous cardiac and immune cell populations, obscuring the cellular origin of candidate genes (17,18). Integrating immune deconvolution and single-cell RNA sequencing provides an avenue to assign prioritized signals to specific immune subsets and microenvironmental niches (19,20).

Here, we integrated bulk and single-cell transcriptomic datasets to prioritize I/R-associated gene programs and determine their cellular sources. By combining cross-dataset validation with complementary prioritization strategies, we identified HMOX1 as a consistently I/R-responsive gene enriched in macrophages. We further validated its upregulation in a rat myocardial I/R model using reverse transcription-quantitative polymerase chain reaction (RT-qPCR), supporting a potential role for HMOX1 in macrophage-associated responses during I/R injury. We present this article in accordance with the ARRIVE reporting checklist (available at https://jtd.amegroups.com/article/view/10.21037/jtd-2026-1-0256/rc).

Methods

Bulk transcriptomic datasets and preprocessing

Bulk myocardial transcriptomic datasets were retrieved from Gene Expression Omnibus (GEO). GSE160516 was used as the discovery cohort and comprised sham-operated mice (n=4) and I/R mice (n=12) sampled at 6, 24, and 72 h post-reperfusion, representing the early-to-subacute reperfusion window. Probe sets mapping to the same gene symbol were collapsed by averaging (limma, aveReps), and between-array normalization was performed using normalizeBetweenArrays. GSE225105 was used for external directional validation. As this dataset included sham, sham + B12, I/R, and I/R + B12 groups, only untreated sham (C1–C3) and untreated I/R samples (IR1–IR3) were analyzed to avoid treatment-related confounding. This dataset represents a single 24 h reperfusion time point with a small sample size (n=3/group); we therefore treated it as directional support rather than definitive standalone replication. Sample information for the discovery and validation bulk transcriptomic datasets is provided in Tables S1,S2. The processed expression matrix was obtained directly from GEO and already contained Gene_id and Gene_name annotations. Expression values were transformed as log₂(FPKM + 1), where FPKM is the fragments per kilobase of transcript per million mapped reads, for gene-level visualization. Group comparisons were performed using two-sided Welch’s t-test, with Benjamini-Hochberg adjustment applied for genome-wide analyses.

Identification of I/R-associated gene modules by WGCNA

WGCNA was conducted in R (version 4.5.2) using the WGCNA package (version 1.73). After filtering low-variance genes and assessing sample quality (goodSamplesGenes), we constructed an unsigned co-expression network based on pairwise Pearson correlations. A soft-thresholding power was selected to approximate scale-free topology. The adjacency matrix was transformed into a topological overlap matrix (TOM), and genes were hierarchically clustered using TOM-based dissimilarity. Modules were detected by dynamic tree cutting (minimum module size =50) and merged using an eigengene dissimilarity threshold of 0.25. Module eigengenes were correlated with I/R status, and gene significance (GS) and module membership (MM) were computed to quantify trait association and intramodular connectivity.

Differential expression analysis and candidate gene selection

Differential expression between I/R and sham groups was assessed with the limma package. Genes with an adjusted P value <0.05 and |log₂ fold change (FC)| >1 were defined as differentially expressed genes (DEGs). DEGs were visualized using volcano plots and heatmaps. Candidate genes were defined as the intersection between DEGs and genes within the key WGCNA module.

Functional enrichment analyses

Functional enrichment analyses were performed using clusterProfiler for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Disease Ontology (DO) enrichment was assessed using disease ontology semantic and enrichment (DOSE). Terms with adjusted P values <0.05 were considered significant.

Protein-protein interaction (PPI) network construction and hub gene identification

PPI information was retrieved from the search tool for the retrieval of interacting genes/proteins (STRING) database, and an interaction score >0.4 was used to construct the PPI network. Networks were visualized in Cytoscape (version 3.9.1). Hub genes were ranked using cytoHubba based on multiple topological metrics [Degree, maximal clique centrality (MCC), maximum neighborhood component (MNC), and edge percolated component (EPC)]; the top-ranked genes from each method were intersected. Molecular complex detection (MCODE) was used to identify densely connected subnetworks within the PPI graph.

Machine learning-based biomarker screening and interpretability analysis

Machine-learning-based prioritization was performed on the expression matrix of the shortlisted hub genes. Given the modest sample size (n=16), we treated modeling primarily as a feature-ranking and interpretability exercise rather than a definitive diagnostic classifier. Least absolute shrinkage and selection operator (LASSO) regression (glmnet) was used to identify features with non-zero coefficients at the optimal penalty. In parallel, a random forest model (randomForest) was fitted and predictor contributions were summarized using SHAP values. Receiver operating characteristic (ROC) curves and areas under the curve (AUCs) were calculated using pROC as descriptive summaries of within-dataset separability, not as estimates of generalizable diagnostic performance. Modeling choices (data splits, cross-validation strategy, and random seeds) were recorded.

Expression validation and immune infiltration analysis

Expression levels of selected biomarkers were compared between I/R and sham groups using two-sided Student’s t-tests, with P<0.05 considered statistically significant. As an exploratory analysis, immune-cell composition was estimated from bulk transcriptomic profiles using CIBERSORT with the LM22 signature matrix. Because LM22 was originally derived from human immune cells, its application to murine bulk data represents a cross-species approximation; results should therefore be interpreted as hypothesis-generating rather than definitive cell-fraction estimates. Correlations between biomarker expression and inferred immune-cell proportions were assessed by Pearson correlation. Effect sizes were summarized and visualized using forest plots (meta package, version 7.0-0).

Construction of an HMOX1-centered upstream regulatory network

Predicted upstream regulators of HMOX1 were integrated to construct an HMOX1-centered regulatory network. Putative microRNAs targeting HMOX1 were obtained by intersecting predictions from multiple databases, and transcription factors potentially regulating HMOX1 were collected from curated regulatory resources. The integrated network was visualized to summarize candidate upstream regulatory relationships.

Single-cell RNA sequencing analysis

Single-cell RNA sequencing data were retrieved from the GEO under accession number GSE288103. This dataset contains myocardial infarction mouse heart samples from Cd248fl/fl control mice and fibroblast-specific Cd248 conditional knockout mice. Sample information is summarized in Table S3. Cells with fewer than 300 detected genes, more than 7,500 detected genes, or a mitochondrial gene proportion exceeding 25% were excluded. Data were normalized (NormalizeData), and highly variable genes were identified (FindVariableFeatures). After scaling, principal component analysis (PCA) was performed. Batch effects were corrected using the Harmony algorithm, and uniform manifold approximation and projection (UMAP) was used for visualization. Graph-based clustering was conducted using FindNeighbors and FindClusters, marker genes were identified with FindAllMarkers, and clusters were annotated using canonical markers.

Rat myocardial ischemia-reperfusion (MI/R) model and RT-qPCR validation

Male Sprague-Dawley rats (6–8 weeks old; Animal Center of Guangxi Medical University) were randomly assigned to sham or I/R groups (n=6/group) after 3–7 days of acclimatization. A 72 h reperfusion interval was chosen to capture the subacute phase of macrophage-associated tissue remodeling, which corresponds to the latest time point sampled in the discovery transcriptomic cohort. Rats were anesthetized, intubated, and ventilated. Through a left thoracotomy (third-fourth intercostal space), the left anterior descending coronary artery (LAD) was ligated with a 5-0 suture ~1–2 mm below the right margin of the left atrial appendage to induce ischemia (confirmed by anterior left ventricle (LV) pallor and ST-segment elevation), followed by 45 min ischemia and 72 h reperfusion; sham animals underwent the same procedure without LAD ligation. At 72 h, myocardial tissues were harvested for RT-qPCR. Total RNA was extracted with TRIzol, cDNA was synthesized after genomic DNA removal, and qPCR was performed using SYBR Green Master Mix with GAPDH as an internal control. Relative HMOX1 expression was calculated using the 2^−ΔΔCt method, and primer sequences are provided in Table S4.

Genetic colocalization and molecular docking analysis

Bayesian genetic colocalization was used to explore whether variants regulating HMOX1 expression overlap with heart failure risk loci. Summary statistics were retrieved from the IEU OpenGWAS database, including HMOX1 expression quantitative trait locu (eQTL) data (eqtl-a-ENSG00000100292) and a heart failure genome-wide association study (GWAS) dataset (ebi-a-GCST009541). Analyses were restricted to the HMOX1 locus on chromosome 22 within a ±50 kb window. Posterior probabilities were computed with the coloc package, and evidence for a shared causal signal was evaluated using PP.H4. For molecular docking, the three-dimensional structure of HMOX1 was obtained from the Protein Data Bank and docked with berberine to predict potential binding conformations. The highest-scoring pose was visualized using molecular graphics software.

Ethical considerations

The animal study was reviewed and approved by the Medical Ethics Committee of The Second Affiliated Hospital of Guangxi Medical University (approval No. KY-0368, 2024). All experimental procedures were conducted in accordance with the institutional guidelines and local legislation for the care and use of laboratory animals.

Statistical analysis

Data analysis was performed in R (v4.5.2), and figures were generated with ggplot2 (v3.5.2). For comparisons involving two groups, Student’s t-test was applied, while one-way analysis of variance (ANOVA) was used when more than two groups were analyzed. A threshold of P<0.05 was considered statistically significant.

Results

Identification of I/R-associated gene modules by WGCNA

Using bulk myocardial expression data from GSE160516 (sham, n=4; I/R, n=12; spanning 6, 24, and 72 h post-reperfusion), we constructed a weighted gene co-expression network. Sample clustering confirmed the absence of obvious outliers, supporting inclusion of all samples (Figure S1). Scale-free topology analysis supported selection of a soft-thresholding power of 7 (Figure 1A). The distribution of MM values across modules is shown (Figure 1B). Dynamic tree cutting and module merging identified 25 co-expression modules (Figure 1C-1E). Among these, the midnightblue module showed the strongest positive correlation with I/R status (r=0.54, P=0.03; Figure 1F). Within this module, GS and MM were strongly correlated (r=0.999, P<1×10⁻²⁰⁰; Figure 1G), supporting the midnightblue module as the key I/R-associated program.

Figure 1 WGCNA identifies a gene module associated with myocardial ischemia-reperfusion injury. (A) Scale-free topology fit and mean connectivity across candidate soft-thresholding powers; power =7 was selected. (B) Distribution of MM values across identified modules. Sample clustering for outlier detection is shown in Figure S1. (C) Gene dendrogram with module assignment obtained by dynamic tree cutting and module merging. (D) Heatmap showing correlations among MEs. (E) Hierarchical clustering of module eigengenes after merging. (F) Module-trait relationships between MEs and I/R status. (G) Relationship between GS and MM for the key module (r=0.999, P=5e−290). GS, gene significance; I/R, ischemia-reperfusion; ME, module eigengenes; MM, module membership; WGCNA, weighted gene co-expression network analysis.

Identification of DEGs associated with myocardial I/R injury

We next compared bulk transcriptomes from sham and I/R myocardium to quantify injury-associated expression changes. PCA separated the two groups (Figure 2A), consistent with a global transcriptional response to I/R. Genome-wide gene ranking based on log₂ fold change is shown in (Figure 2B). Differential expression analysis identified DEGs using |log₂FC| >1 and adjusted P<0.05. These DEGs showed consistent within-group patterns and strong between-group separation in heatmaps (Figure 2C), and the volcano plot highlighted prominent activation of stress- and inflammation-related transcripts in I/R samples (Figure 2D).

Figure 2 Differential gene expression analysis of sham and ischemia-reperfusion myocardium. (A) PCA showing separation between sham and I/R samples. (B) Genome-wide gene ranking based on log₂ fold change between I/R and sham groups. (C) Heatmap of representative DEGs across sham (n=4) and I/R (n=12) samples (row-wise z-score scaling). (D) Volcano plot of DEGs using fold change and adjusted P value thresholds. DEGs, differentially expressed genes; I/R, ischemia-reperfusion; PCA, principal component analysis.

Functional enrichment analysis highlights immune and inflammatory involvement

Enrichment analyses of I/R-associated genes highlighted immune activation and inflammatory responses among GO biological process terms (Figure 3A), while cellular component categories were linked to membrane- and cytoskeleton-related structures (Figure 3B). A GO term network visualization summarizing enriched immune-related processes is shown in (Figure 3C). KEGG analysis identified multiple immune pathways, including cytokine-cytokine receptor interaction, chemokine signaling, and Toll-like receptor signaling (Figure 3D). DO enrichment further associated the gene set with ischemia- and inflammation-related conditions (Figure 3E).

Figure 3 Functional enrichment of I/R-associated genes. (A) GO biological process enrichment highlighting immune- and inflammation-related terms. (B) GO cellular component enrichment linking candidate genes to membrane and cytoskeletal structures. (C) GO term network visualization summarizing relationships among significantly enriched immune-related processes. (D) KEGG pathway enrichment identifying key immune and inflammatory signalling pathways. (E) DO enrichment indicating associations with ischemia- and inflammation-related conditions. DO, Disease Ontology; GO, Gene Ontology; I/R, ischemia-reperfusion; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network analysis identifies hub genes related to myocardial I/R injury

To examine protein-level connectivity, we built a STRING-derived PPI network comprising 98 nodes and 734 edges (Figure 4A). The hub-gene interaction landscape is shown in (Figure 4B). MCODE identified a densely connected subnetwork of 16 genes (Figure 4C). Hub prioritization using multiple cytoHubba metrics converged on 14 overlapping candidates (Figure 4D).

Figure 4 PPI network analysis and hub-gene identification for myocardial ischemia-reperfusion injury. (A) STRING-derived PPI network constructed from I/R-associated genes. (B) PPI network with hub genes highlighted. (C) High-density subnetwork identified by MCODE analysis. (D) Venn diagram showing overlap among the top-ranked genes identified by cytoHubba (MCC, MNC, and Degree); intersecting genes were defined as hub genes. EPC, edge percolated component; I/R, ischemia-reperfusion; MCC, maximal clique centrality; MCODE, molecular complex detection; MNC, maximum neighborhood component; PPI, protein-protein interaction; STRING, search tool for the retrieval of interacting genes/proteins.

Machine learning prioritizes HMOX1 and VCAM1 as core candidate markers

From the hub-gene panel, machine-learning-based feature selection and interpretability analyses were applied to prioritize a compact set of candidate markers. LASSO regression with cross-validation retained three genes (VCAM1, HMOX1, and Actb) with non-zero coefficients at the optimal penalty (log λ =−4.0; Figure 5A,5B). Actb was retained as a housekeeping gene whose compositional signal likely reflects overall cellular content rather than a specific biological response to I/R, and was therefore not advanced as a biological candidate. SHAP-based analysis consistently highlighted VCAM1 and HMOX1 as the most influential predictors across models (Figure 5C-5E). Only VCAM1 and HMOX1 were supported by both LASSO selection and SHAP prioritization (Figure 5F). A simple neural network built from these two genes separated sham and I/R samples within this dataset (Figure 5G); this classifier is presented as a descriptive illustration of feature separability rather than a validated diagnostic tool.

Figure 5 Machine-learning-based prioritization of core candidate markers. (A) LASSO cross-validation curve used to select the optimal penalty parameter (log λ =−4.0). (B) LASSO coefficient profiles showing genes with non-zero coefficients at the optimal λ. (C) Feature-importance ranking from SHAP analysis. (D) SHAP waterfall plot illustrating feature contributions for a representative sample. (E) SHAP dependence plots showing relationships between gene expression and SHAP values. (F) Overlap between genes selected by LASSO and those highlighted by SHAP. (G) Architecture of the ANN model built from the core biomarkers. ANN, artificial neural network; LASSO, least absolute shrinkage and selection operator; SHAP, shapley additive explanations.

HMOX1 is linked to immune infiltration patterns in myocardial I/R injury

Expression analysis in the discovery cohort (GSE160516) showed that both HMOX1 and VCAM1 were upregulated in I/R myocardium relative to sham controls (Figure 6A). To examine whether these signals were reproducible, we analyzed an independent bulk RNA-seq cohort (GSE225105) representing a single 24 h reperfusion time point. In this small external cohort (n=3/group), HMOX1 showed a directionally concordant increase in I/R samples; however, given the limited sample size, the difference did not reach statistical significance. We therefore interpret this dataset as directional support rather than definitive standalone replication. VCAM1 displayed a similar trend with greater variability. Effect sizes across the discovery dataset are summarized in Figure 6B. Exploratory CIBERSORT analysis—noting that the human-derived LM22 signature matrix represents a cross-species approximation when applied to murine data—suggested substantial remodeling of the immune microenvironment in I/R myocardium, including increased naïve B cells and M0 macrophages and decreased resting mast cells (all P<0.05). Pearson correlation analysis further indicated that HMOX1 expression was positively associated with several inferred immune populations, most prominently M0 macrophages (r=0.67), whereas associations for VCAM1 were weaker and less consistent (Figure 6C). To explore potential upstream regulation, an HMOX1-centered regulatory network integrating predicted microRNAs and transcription factors was constructed, highlighting miR-7215-5p and several transcriptional regulators as candidate modulators (Figure 6D) (r=0.67).

Figure 6 Immune infiltration and upstream regulatory analysis of HMOX1. (A) Expression of HMOX1 and VCAM1 in the discovery dataset (GSE160516) and an independent validation dataset (GSE225105). (B) Forest plot summarizing effect sizes (95% CI) for the core genes. (C) Correlations between HMOX1 expression and inferred immune-cell proportions; point size indicates the absolute value of the correlation coefficient (r), and colour reflects the P value range. (D) HMOX1-centred upstream regulatory network integrating predicted microRNAs and transcription factors. CI, confidence interval; HR, hazard ratio; NK, natural killer; RMSE, root mean square error.

Single-cell analysis localizes HMOX1 expression to macrophages

Single-cell RNA sequencing was used to define the cellular landscape of I/R myocardium and localize HMOX1 expression. After quality control and integration, 39,934 high-quality cells were retained and visualized by UMAP (Figure 7A). Based on canonical marker genes, 11 major cell types were annotated (Figure 7B), with fibroblasts, endothelial cells, and macrophages comprising the predominant compartments across samples (Figure 7C). Notably, HMOX1 expression was enriched in macrophages, supporting macrophages as the main cellular source of HMOX1 in the I/R heart.

Figure 7 Single-cell transcriptomic landscape of myocardial tissue after injury. (A) UMAP projection of 39,934 myocardial cells coloured by annotated cell type. (B) Feature plots of canonical marker genes used for cluster annotation. (C) Cell-type composition across samples. UMAP, uniform manifold approximation and projection.

Experimental validation in a rat MI/R model confirms HMOX1 upregulation

To validate the computational findings in an independent in vivo system, a rat MI/R model was established (Figure 8A). RT-qPCR confirmed higher myocardial HMOX1 mRNA levels in MI/R rats than in sham-operated controls (Figure 8B), supporting cross-species concordance with the transcriptomic analyses.

Figure 8 Rat MI/R model and RT-qPCR validation of HMOX1. (A) Representative images of LAD occlusion and reperfusion during MI/R surgery. (B) RT-qPCR quantification of HMOX1 mRNA expression (2^−ΔΔCt method; mean ± SEM). *, P<0.05 vs. sham. LAD, left anterior descending coronary artery; MI/R, myocardial ischemia-reperfusion; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; SEM, standard error of the mean.

Genetic colocalization and molecular docking provide hypothesis-generating support for HMOX1

To explore translational relevance from both genetic and structural perspectives, we performed human genetic colocalization and structure-based docking analyses. Colocalization suggested overlapping association signals between HMOX1 eQTLs and heart failure GWAS summary statistics across the HMOX1 locus on chromosome 22 (Figure 9A), consistent with a shared regulatory component. Docking indicated that berberine could occupy a defined pocket within HMOX1 and adopt a plausible binding conformation (Figure 9B,9C). These analyses provide hypothesis-generating support for translational relevance but do not establish causality or pharmacological efficacy.

Figure 9 Translational analyses of HMOX1. (A) Regional association and Bayesian colocalization plot for the HMOX1 locus comparing eQTL and heart failure GWAS signals. (B) Predicted docking pose of berberine in the HMOX1 structure. (C) Close-up view of the predicted binding pocket. eQTL, expression quantitative trait locus; GWAS, genome-wide association study.

Discussion

MI/R injury elicits a burst of innate immune activity that can amplify cardiomyocyte loss and shape the trajectory of tissue remodeling (21,22). Here, we adopted a multi-tiered analytical strategy—combining bulk co-expression networks, differential expression, protein-interaction topology, model-based feature ranking, immune profiling, single-cell mapping, in vivo validation, and exploratory translational analyses—to nominate I/R-associated genes with defined cellular context (23). Because the discovery cohort spanned 6, 24, and 72 h post-reperfusion, our findings reflect the early-to-subacute reperfusion window rather than a single-time-point snapshot. Across these analytical layers, HMOX1 and VCAM1 emerged as informative candidates, with the most consistent evidence supporting HMOX1 as a macrophage-enriched I/R-responsive gene (24).

At the network level, WGCNA highlighted the midnightblue module as tightly linked to the I/R phenotype. The strong concordance between GS and MM suggested a coherent transcriptional program in which highly trait-associated genes also occupy central network positions. Intersecting this module with differential expression results helped focus downstream analyses on genes that were both network-coherent and robustly altered after I/R. Enrichment analyses consistently pointed to innate immune activation, leukocyte chemotaxis, and Toll-like receptor-related signaling, reinforcing inflammation as a core transcriptional signature of reperfusion injury (25,26).

Within this candidate pool, PPI topology and complementary centrality metrics identified a compact set of hub genes. Feature-selection models added a second, orthogonal layer of ranking: LASSO regression retained three genes at the optimal penalty (VCAM1, HMOX1, and Actb), while SHAP analysis independently highlighted the same top two—HMOX1 and VCAM1—as the most discriminatory features. Actb was interpreted as a housekeeping-related signal reflecting overall cellular content rather than a specific I/R biological response, and was not pursued further. VCAM1 encodes an endothelial adhesion molecule central to leukocyte recruitment (27,28), and its selection aligns with the inflammatory character of the gene set. HMOX1, the inducible heme-degrading enzyme, occupies the interface of oxidative stress and inflammatory control through its generation of biliverdin, carbon monoxide, and ferrous iron. We ultimately focused on HMOX1 for downstream analyses because its correlation with macrophage-related signatures in bulk data was markedly stronger than that of VCAM1, and because single-cell profiling confirmed macrophage-predominant expression (29-32).

Cellular context was critical for interpretation. In bulk profiles, HMOX1 tracked inferred immune-cell shifts and showed the strongest association with macrophage-related signals (33). Single-cell mapping further localized HMOX1 expression predominantly to macrophages within injured myocardium. This pattern is consistent with the time-dependent roles of macrophages in I/R injury, where early inflammatory programs can exacerbate damage while later reparative programs support healing (34). The upstream regulatory network provides testable hypotheses for transcriptional and post-transcriptional control of HMOX1 in this setting (35,36).

Cross-species validation in an independent rat MI/R model confirmed elevated HMOX1 transcript levels at 72 h post-reperfusion, a time point chosen to capture subacute macrophage-associated tissue remodeling and corresponding to the latest phase sampled in the discovery cohort (37). An external murine bulk RNA-seq dataset (GSE225105; 24 h reperfusion, n=3/group) showed a directionally concordant HMOX1 increase that did not reach statistical significance; given the small sample size, this cohort provides trend-level support rather than definitive replication. The human colocalization and docking analyses were included as hypothesis-generating exercises that may guide future translational work; they do not establish causality or pharmacological efficacy (38).

Several limitations warrant consideration. First, the discovery cohort pooled samples across 6, 24, and 72 h post-reperfusion; while this breadth may enhance the detection of genes robustly regulated across the early-to-subacute window, it also introduces time-point heterogeneity that could obscure phase-specific signals. Second, the discovery dataset contained only four sham samples against 12 I/R samples, which may increase variance in network inference and feature ranking; we mitigated this through external-dataset directional validation and cross-species RT-qPCR. Third, bulk transcriptomic profiles represent cellular mixtures, and the CIBERSORT/LM22-based deconvolution applied here uses a human-derived signature matrix on murine data—a cross-species approximation whose results should be interpreted cautiously (39). Although single-cell mapping confirmed macrophage-enriched HMOX1 expression, it remains unclear whether the observed HMOX1 upregulation reflects increased macrophage abundance, elevated per-cell transcription, or both; direct protein-level quantification and flow-cytometric validation in matched specimens will be needed to resolve this distinction. Fourth, genetic colocalization and molecular docking carry inherent assumptions and should be treated as hypothesis-generating rather than confirmatory; independent functional experiments are essential to establish mechanistic causality (40,41).

In summary, by layering co-expression analysis, feature-selection modeling, single-cell mapping, and cross-species experimental validation, this study identifies HMOX1 as a macrophage-enriched, reproducibly I/R-responsive marker in the heart. The principal contribution lies not in HMOX1 as a novel gene—its cytoprotective functions are well established—but in the systematic, multi-platform demonstration that its I/R-responsive expression maps specifically to the macrophage compartment. These findings provide a defined starting point for mechanistic studies and for evaluating whether targeted modulation of HMOX1-related pathways can attenuate macrophage-driven reperfusion damage.

Conclusions

Integrated bulk and single-cell transcriptomic analyses, reinforced by cross-species RT-qPCR validation, consistently prioritize HMOX1 as a macrophage-enriched, I/R-responsive marker. These findings establish a cell-type-resolved candidate and motivate future functional studies to determine whether HMOX1 actively regulates macrophage phenotypic switching during reperfusion injury and whether its modulation confers measurable cardioprotection.

Acknowledgments

The authors thank Citexs for English language editing and Stork for linguistic assistance and literature monitoring. The authors also acknowledge Home for Researchers for general academic support during manuscript preparation, as well as members of GZLab for valuable advice and assistance with bioinformatics analyses.

Footnote

Reporting Checklist: The authors have completed the ARRIVE reporting checklist. Available at https://jtd.amegroups.com/article/view/10.21037/jtd-2026-1-0256/rc

Data Sharing Statement: Available at https://jtd.amegroups.com/article/view/10.21037/jtd-2026-1-0256/dss

Peer Review File: Available at https://jtd.amegroups.com/article/view/10.21037/jtd-2026-1-0256/prf

Funding: This work was supported by the Guangxi Natural Science Foundation (No. 2024GXNSFAA010463), the Guangxi Medical and Health Key Discipline Construction Project, and the 2025 Guangxi Higher Education Institutions Young and Middle-aged Teachers’ Research Basic Ability Enhancement Project (No. 2025KY0122).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jtd.amegroups.com/article/view/10.21037/jtd-2026-1-0256/coif). X.B. reports that this study has received funding support from 2025 Guangxi Higher Education Institutions Young and Middle-aged Teachers’ Research Basic Ability Enhancement Project (grant No. 2025KY0122). L.X. reports that this study has received funding support from Guangxi Natural Science Foundation (grant No. 2024GXNSFAA010463) and Guangxi Medical and Health Key Discipline Construction Project. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The animal study was reviewed and approved by the Medical Ethics Committee of The Second Affiliated Hospital of Guangxi Medical University (approval No. KY-0368, 2024). All experimental procedures were conducted in accordance with the institutional guidelines and local legislation for the care and use of laboratory animals.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Byrne RA, Rossello X, Coughlan JJ, et al. 2023 ESC Guidelines for the management of acute coronary syndromes. Eur Heart J 2023;44:3720-826. [Crossref] [PubMed]
2. Rao SV, O’Donoghue ML, Ruel M, et al. 2025 ACC/AHA/ACEP/NAEMSP/SCAI Guideline for the Management of Patients With Acute Coronary Syndromes: A Report of the American College of Cardiology/American Heart Association Joint Committee on Clinical Practice Guidelines. Circulation 2025;151:e771-862. [Crossref] [PubMed]
3. Yellon DM, Hausenloy DJ. Myocardial reperfusion injury. N Engl J Med 2007;357:1121-35. [Crossref] [PubMed]
4. Hausenloy DJ, Yellon DM. Myocardial ischemia-reperfusion injury: a neglected therapeutic target. J Clin Invest 2013;123:92-100. [Crossref] [PubMed]
5. Heusch G. Myocardial ischaemia-reperfusion injury and cardioprotection in perspective. Nat Rev Cardiol 2020;17:773-89. [Crossref] [PubMed]
6. Frangogiannis NG. The inflammatory response in myocardial injury, repair, and remodelling. Nat Rev Cardiol 2014;11:255-65. [Crossref] [PubMed]
7. Swirski FK, Nahrendorf M. Leukocyte behavior in atherosclerosis, myocardial infarction, and heart failure. Science 2013;339:161-6. [Crossref] [PubMed]
8. Epelman S, Lavine KJ, Randolph GJ. Origin and functions of tissue macrophages. Immunity 2014;41:21-35. [Crossref] [PubMed]
9. Murray PJ, Wynn TA. Protective and pathogenic functions of macrophage subsets. Nat Rev Immunol 2011;11:723-37. [Crossref] [PubMed]
10. Bajpai G, Schneider C, Wong N, et al. The human heart contains distinct macrophage subsets with divergent origins and functions. Nat Med 2018;24:1234-45. [Crossref] [PubMed]
11. Wynn TA, Vannella KM. Macrophages in Tissue Repair, Regeneration, and Fibrosis. Immunity 2016;44:450-62. [Crossref] [PubMed]
12. Nahrendorf M, Swirski FK. Monocyte and macrophage heterogeneity in the heart. Circ Res 2013;112:1624-33. [Crossref] [PubMed]
13. Nahrendorf M, Swirski FK, Aikawa E, et al. The healing myocardium sequentially mobilizes two monocyte subsets with divergent and complementary functions. J Exp Med 2007;204:3037-47. [Crossref] [PubMed]
14. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
15. Tibshirani R. Regression Shrinkage and Selection Via the Lasso. Journal of the Royal Statistical Society: Series B (Methodological) 1996;58:267-88.
16. Lundberg SM, Lee S-I. A unified approach to interpreting model predictions. Proceedings of the 31st International Conference on Neural Information Processing Systems; Long Beach, California, USA: Curran Associates Inc.; 2017:4768-77.
17. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
18. Stuart T, Butler A, Hoffman P, et al. Comprehensive Integration of Single-Cell Data. Cell 2019;177:1888-1902.e21. [Crossref] [PubMed]
19. Litviňuková M, Talavera-López C, Maatz H, et al. Cells of the adult human heart. Nature 2020;588:466-72. [Crossref] [PubMed]
20. Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 2019;16:1289-96. [Crossref] [PubMed]
21. Timmers L, Pasterkamp G, de Hoog VC, et al. The innate immune response in reperfused myocardium. Cardiovasc Res 2012;94:276-83. [Crossref] [PubMed]
22. Arslan F, de Kleijn DP, Pasterkamp G. Innate immune signaling in cardiac ischemia. Nat Rev Cardiol 2011;8:292-300. [Crossref] [PubMed]
23. Prabhu SD, Frangogiannis NG. The Biological Basis for Cardiac Repair After Myocardial Infarction: From Inflammation to Fibrosis. Circ Res 2016;119:91-112. [Crossref] [PubMed]
24. Weinberger T, Denise M, Joppich M, et al. Resident and recruited macrophages differentially contribute to cardiac healing after myocardial ischemia. Elife 2024;12:RP89377. [Crossref] [PubMed]
25. Arslan F, Smeets MB, O’Neill LA, et al. Myocardial ischemia/reperfusion injury is mediated by leukocytic toll-like receptor-2 and reduced by systemic administration of a novel anti-toll-like receptor-2 antibody. Circulation 2010;121:80-90. [Crossref] [PubMed]
26. Shimamoto A, Chong AJ, Yada M, et al. Inhibition of Toll-like receptor 4 with eritoran attenuates myocardial ischemia-reperfusion injury. Circulation 2006;114:I270-4. [Crossref] [PubMed]
27. Ley K, Laudanna C, Cybulsky MI, et al. Getting to the site of inflammation: the leukocyte adhesion cascade updated. Nat Rev Immunol 2007;7:678-89. [Crossref] [PubMed]
28. Cybulsky MI, Gimbrone MA Jr. Endothelial expression of a mononuclear leukocyte adhesion molecule during atherogenesis. Science 1991;251:788-91. [Crossref] [PubMed]
29. Yet SF, Melo LG, Layne MD, et al. Heme oxygenase 1 in regulation of inflammation and oxidative damage. Methods Enzymol 2002;353:163-76. [Crossref] [PubMed]
30. Otterbein LE, Bach FH, Alam J, et al. Carbon monoxide has anti-inflammatory effects involving the mitogen-activated protein kinase pathway. Nat Med 2000;6:422-8. [Crossref] [PubMed]
31. Gozzelino R, Jeney V, Soares MP. Mechanisms of cell protection by heme oxygenase-1. Annu Rev Pharmacol Toxicol 2010;50:323-54. [Crossref] [PubMed]
32. Ryter SW. Heme Oxygenase-1: An Anti-Inflammatory Effector in Cardiovascular, Lung, and Related Metabolic Disorders. Antioxidants (Basel) 2022;11:555. [Crossref] [PubMed]
33. Ryter SW, Alam J, Choi AM. Heme oxygenase-1/carbon monoxide: from basic science to therapeutic applications. Physiol Rev 2006;86:583-650. [Crossref] [PubMed]
34. Sun J, Hoshino H, Takaku K, et al. Hemoprotein Bach1 regulates enhancer availability of heme oxygenase-1 gene. EMBO J 2002;21:5216-24. [Crossref] [PubMed]
35. Kensler TW, Wakabayashi N, Biswal S. Cell survival responses to environmental stresses via the Keap1-Nrf2-ARE pathway. Annu Rev Pharmacol Toxicol 2007;47:89-116. [Crossref] [PubMed]
36. Yamamoto M, Kensler TW, Motohashi H. The KEAP1-NRF2 System: a Thiol-Based Sensor-Effector Apparatus for Maintaining Redox Homeostasis. Physiol Rev 2018;98:1169-203. [Crossref] [PubMed]
37. Duckers HJ, Boehm M, True AL, et al. Heme oxygenase-1 protects against vascular constriction and proliferation. Nat Med 2001;7:693-8. [Crossref] [PubMed]
38. Giambartolomei C, Vukcevic D, Schadt EE, et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet 2014;10:e1004383. [Crossref] [PubMed]
39. Avila Cobos F, Alquicira-Hernandez J, Powell JE, et al. Benchmarking of cell type deconvolution pipelines for transcriptomics data. Nat Commun 2020;11:5650. [Crossref] [PubMed]
40. Garmire LX, Li Y, Huang Q, et al. Challenges and perspectives in computational deconvolution of genomics data. Nat Methods 2024;21:391-400. [Crossref] [PubMed]
41. Huang Z, Han Z, Ye B, et al. Berberine alleviates cardiac ischemia/reperfusion injury by inhibiting excessive autophagy in cardiomyocytes. Eur J Pharmacol 2015;762:1-10. [Crossref] [PubMed]