Molecular Mechanisms of DBNL in Heart Failure: From Macrophage Immunometabolism to Therapeutic Implications

Summary

We present a protocol to integrate multiomics Mendelian randomization, single-cell transcriptomics, and structure-based in silico screening to systematically delineate the Drebrin Like (DBNL)-macrophage axis in heart failure. We suggest DBNL as a druggable immunometabolic target.

Abstract

Heart failure (HF) remains a major global health challenge, with limited effective treatments targeting its core pathophysiological mechanisms. In this study, an integrated multiomics approach combining Mendelian randomization (MR) and single-cell RNA sequencing (scRNA-seq) was used to identify potential biomarkers and therapeutic targets for HF. We utilized data from genome-wide association studies (GWASs), expression quantitative trait loci (eQTLs), methylation quantitative trait loci (mQTLs), and protein quantitative trait loci (pQTLs) to investigate the genetic mechanisms of HF. Single-cell RNA sequencing was performed to analyse gene expression in cardiac macrophages, and pseudotime analysis was used to study the dynamic regulation of DBNL during macrophage differentiation. Molecular docking and dynamics simulations identified pirinixic acid (WY-14643; PubChem CID: 4594; synonyms: 50892-23-4; WY-14643) as a potential regulator of DBNL. Multiomics analysis revealed that DBNL (Drebrin-like) is a key gene associated with HF risk. Single-cell RNA sequencing revealed that DBNL is expressed mainly in cardiac macrophages and is upregulated under pathological conditions. The expression of DBNL by macrophages is associated with immune metabolism and profibrotic pathways, particularly the IL-6/JAK/STAT3, PI3K/AKT/mTOR, and TGF-β signalling pathways. Pseudotime analysis indicated that DBNL has a dynamic regulatory effect on macrophage differentiation, especially in chronic inflammation. Molecular docking and dynamic simulations have shown that pyridine acid has a potential role in regulating DBNL. This study elucidates the role of DBNL in the progression of heart failure and suggests its potential as a therapeutic target. Importantly, the therapeutic relevance is inferred from computational predictions. These findings offer novel insights into immune metabolism and macrophage-mediated intercellular communication, establishing a theoretical basis for future applications in the optimization of DBNL-targeted drugs and functional research. Nonetheless, additional experimental validation and preclinical studies are needed to substantiate the clinical efficacy of DBNL as a therapeutic target.

Introduction

Heart failure (HF) represents a significant and growing global health challenge that affects millions of individuals worldwide and places considerable economic strain on healthcare systems¹. Consequently, there is a pressing need for efficacious interventions in heart failure (HF) management.

The complex pathophysiology of heart failure originates from impaired myocardial contractility of the left ventricle, leading to a cascade of systemic complications, such as pulmonary and peripheral congestion². However, directly addressing the fundamental pathophysiological mechanisms that underlie the progression of heart failure, particularly adverse myocardial remodelling, remains a critical and unmet challenge³. In recent decades, the domain of biomedical research has undergone a substantial paradigm shift, primarily facilitated by the incorporation of multiomics technologies and advanced bioinformatics methodologies. These tools offer exceptional opportunities for identifying novel disease biomarkers and therapeutic targets. Mendelian randomization (MR), a rigorous genetic epidemiological approach, uses single-nucleotide polymorphisms (SNPs) as instrumental variables to infer causal relationships between genetic exposures and clinical outcomes⁴. By adopting the methodological framework of randomized controlled trials, Mendelian randomization (MR) effectively addresses the problems of confounding and reverse causality to provide robust evidence for the identification and repurposing of drug targets. In parallel, single-cell RNA sequencing (scRNA-seq) reveals tissue heterogeneity with unparalleled resolution,offering profound insights into cell type-specific gene expression and intercellular communication networks within intricate organs such as the heart⁵.The multiomics approach, which synthesizes genetic, transcriptomic, and proteomic data, offers substantial potential for elucidating disease mechanisms and identifying therapeutic pathways.

Despite these technological advances, significant knowledge gaps persist regarding the specific cellular and molecular mechanisms that drive cardiac inflammation and remodelling in heart failure. Accumulating evidence suggests that immune cells, particularly macrophages, are crucial for orchestrating inflammatory responses and subsequent tissue remodelling during myocardial injury and chronic heart failure processes, which are in turn crucial to disease progression⁶. Previous research has concentrated predominantly on broad inflammatory mediators or cell surface markers, resulting in a significant gap in the exploration of intracellular regulatory targets within these crucial immune cell populations. Additionally, although traditional genome-wide association studies (GWASs) have identified numerous genetic associations related to heart failure, they frequently encounter challenges in establishing definitive causal relationships, thereby limiting their applicability for direct therapeutic translation⁷.Consequently, the implementation of a comprehensive multiomics Mendelian randomization (MR) framework, which incorporates genetic variants that influence DNA methylation (mQTLs), gene expression (eQTLs), and protein abundance (pQTLs), is essential for achieving more robust causal inferences across multiple biological regulatory layers. This approach addresses the inherent limitations associated with single-layer QTL studies. Notably, the Drebrin-like (DBNL) protein has been the subject of extensive research in the context of neuronal development and cytoskeletal dynamics⁸; however, its functional significance and therapeutic potential in cardiovascular diseases, particularly heart failure, remain largely unexplored. This significant gap in our understanding impedes the consideration of DBNL as a novel and potentially viable target for drug development in heart failure treatment.

To address these critical knowledge gaps, this study aims to robustly identify key genes causally associated with heart failure risk through an integrated multiomics MR framework that combines eQTL, mQTL, and pQTL data. This study presents a novel approach involving the integration of multiomics Mendelian randomization analyses -- including expression quantitative trait loci (eQTLs), methylation quantitative trait loci (mQTLs), and protein quantitative trait loci (pQTLs) -- with high-resolution single-cell transcriptomic data (GSE161470) and bulk transcriptomic data (GSE161472) to systematically identify the DBNL gene as a critical causal risk factor for heart failure. These findings demonstrate that DBNL is predominantly and specifically expressed in cardiac macrophages and that its expression changes dynamically throughout their differentiation process. Moreover, the results of this study elucidate the pivotal role of DBNL in modulating the cardiac immune microenvironment and orchestrating heart failure-associated immunometabolic pathways, such as the IL-6/JAK/STAT3 signalling cascade, as well as profibrotic signalling pathways. Additionally, by employing computational biology techniques, this research predicts that pyridyl acid (WY-14643) is a high-affinity ligand for DBNL, thereby providing a comprehensive, multidimensional theoretical framework that supports the designation of DBNL as a novel therapeutic target for heart failure and facilitates the development of targeted pharmacological interventions. Using single-cell RNA sequencing, we characterized the cell-specific expression and functional pathways of the key gene DBNL identified within the cardiac microenvironment, and we investigated the role of DBNL in immunometabolic and profibrotic signalling pathways, as well as in cell-cell communication among cardiac cell types. Potential small-molecule ligands for DBNL were identified through computational molecular docking and molecular dynamics simulations. We hypothesize that DBNL, which is expressed primarily in cardiac macrophages, exerts a causal effect on heart failure progression by regulating immunometabolic pathways and intercellular crosstalk. In summary, by integrating advanced multiomics Mendelian randomization with high-resolution single-cell analysis, we have systematically identified and characterized novel therapeutic targets for heart failure and explored potential small-molecule modulator targets.

Protocol

Data Download
Gene Expression Data
The single-cell RNA sequencing (scRNA-seq) data utilized in the present study were sourced from the Gene Expression Omnibus (GEO) repository maintained by the National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/geo/), specifically from the dataset with the accession number⁹ GSE161470 (human cardiac tissue comprising four control specimens and one pathological specimen). This dataset was originally published by Zhang et al. in 2022¹⁰. The primary objective of the original investigation was to examine the cellular heterogeneity and molecular regulatory mechanisms in human cardiac tissue under conditions of heart failure. For the current analysis, five samples from this dataset, each comprising comprehensive single-cell expression profiles derived from human heart tissue, were selected. The other dataset utilized in this study was likewise obtained from the NCBI GEO public repository, specifically the Series Matrix File corresponding to accession number GSE161472, accompanied by the annotation file GPL11154. The expression profile comprises a total of 84 samples, including 37 control samples and 47 disease samples. This research encompasses an integrative multiomics analysis, with all investigations conducted using publicly accessible data.

eQTL data
The eQTL data, obtained from the eQTLGen consortium, focus on elucidating the genetic framework of gene expression in blood and the genetic factors influencing complex traits³. The consortium is presently engaged in the second phase of its extensive project, undertaking meta-analyses of genome-wide data pertaining to blood gene expression.

Exposure data - mQTLs
The mQTL data were obtained from a published meta-analysis of the European (EUR) cohort, which examines whole-blood DNA methylation within the genetic framework of 3,701 samples from populations of European ancestry¹¹. The dataset included information on 426,636 mQTL traits.

Exposure data-pQTL
Plasma pQTL data were acquired from the deCODE database (https://www.decode.com/summarydata/)⁴. This study employed the 2021 data release of the deCODE pQTL dataset, encompassing a genome-wide association study (GWAS) of plasma protein levels measured using 4,907 aptamers in a cohort of 35,559 individuals of European descent.

Outcome data
Summary statistics for heart failure were obtained from a large-scale genome-wide association study (GWAS) primarily involving participants of European ancestry, accessed through the European Bioinformatics Institute (EBI) database (GCST90162626). The heart failure dataset included 115,150 cases and 1,550,331 controls. The GWAS Catalogue, which encompasses publications, leading associations, and detailed summary statistics, currently offers data mapped to the Genome Assembly and dbSNP Build.

Mendelian randomization analysis of mQTLs, eQTLs, and pQTLs
To systematically investigate the potential causal relationships between gene expression, protein abundance, DNA methylation levels, and the risk of heart failure, Mendelian randomization (MR) analyses were performed utilizing expression quantitative trait loci (eQTLs), protein quantitative trait loci (pQTLs), and methylation quantitative trait loci (mQTLs) datasets. During the exposure data preprocessing phase, single-nucleotide polymorphisms (SNPs) associated with each exposure variable (gene, protein, or methylation site) were extracted from the respective databases at a genome-wide significance threshold of P < 1 × 10⁻⁵ to serve as initial candidate instrumental variables (IVs). Subsequently, linkage disequilibrium (LD) clumping was conducted for each exposure factor's IVs using a window size of 10,000 kilobases (kb) and an LD R² threshold of 0.001 to ensure independence among instruments. These selected IVs were then harmonized with summary statistics from a heart failure genome-wide association study (GWAS; ID: GCST90162626) by employing the read_outcome_data function, retaining only those SNPs exhibiting an association P value below 5×10⁻⁵ in the outcome dataset. To mitigate weak instrument bias, the F statistic for each IV was calculated as F = (β_exposure/SE_exposure)², and only instruments with F > 10 were included in subsequent analyses. For causal effect estimation, allele alignment between the exposure and outcome datasets was performed using the harmonize_data function from the TwoSampleMR package. MR analyses were then conducted employing four complementary statistical approaches: (1) the inverse variance weighted (IVW) method, which provides meta-analyses of Wald ratio estimates across SNPs; (2) MR-Egger regression, which accounts for directional pleiotropy by incorporating an intercept term under the Instrument Strength Independent of Direct Effect (InSIDE) assumption; (3) the weighted median method, which yields consistent causal estimates even if up to 50% of the instruments are invalid; and (4) the weighted mode method, which identifies the most frequent causal effect estimate cluster, offering enhanced statistical power and reduced type I error relative to MR-Egger. In instances where only a single instrumental variable was available, the Wald ratio method was applied exclusively. To assess the robustness of the findings, comprehensive sensitivity analyses were performed, including heterogeneity testing via the mr_heterogeneity function, pleiotropy assessment using mr_pleiotropy_test, and leave-one-out analyses implemented through the mr_leaveoneout function, which iteratively excludes each SNP to determine the influence of individual variants on the overall results. Significant associations were visualized using graphical tools such as mr_scatter_plot and mr_forest_plot. This analytical pipeline was uniformly applied across the eQTL, pQTL, and mQTL datasets to maintain methodological consistency throughout the study.

Colocalization analysis
A colocalization analysis was conducted through the Coloc method, eQTL summary data, and a GWAS of heart failure⁵. The index single-nucleotide polymorphism (SNP) was utilized to compute the posterior probability within a 100 kb clumping window. In the colocalization (coloc) analysis, Hypothesis H3 denotes the posterior probability that the two traits, namely, gene expression and heart failure, are correlated but possess distinct causal variants. Conversely, Hypothesis H4 indicates the posterior probability that the association between the two traits is attributable to a single, shared causal variant. A threshold of SNP.PP.H4 greater than 0.90 was used to determine colocalization.

Immune infiltration
The CIBERSORT method is a widely adopted technique for assessing immune cell types within the microenvironment¹². By utilizing support vector regression principles, deconvolution analysis of the expression matrix of immune cell subtypes can be performed. By incorporating 547 biomarkers, CIBERSORT can differentiate 22 human immune cell phenotypes, including T cells, B cells, plasma cells, and various myeloid cell subpopulations. Utilizing the GSE161472 dataset, an analysis was conducted employing the CIBERSORT algorithm in conjunction with its integrated LM22 signature matrix, which characterizes the gene expression profiles of 22 distinct human immune cell types. The infiltration levels of these 22 immune cell populations were quantified for each individual sample. Subsequently, the cor.test function was applied to assess the correlations between the expression of key genes and the corresponding immune cell infiltration levels.

Single-cell RNA sequencing data processing and quality control
The single-cell expression profile data were processed using the Seurat (V4.3.0) package in the R (V4.3.0) environment⁶.This study utilized a conventional workflow for the analysis of single-cell RNA sequencing data. Initially, the expression profiles were imported utilizing the Seurat package. Cells were filtered on the basis of several quality metrics, including each cell's total UMI count, number of expressed genes, proportion of mitochondrial reads, and proportion of ribosomal reads. Outliers were identified as values deviating from the median by more than three median absolute deviations (MADs). The specific filtering thresholds applied were as follows: nFeature_RNA ≥ 200, percent.mt ≤ 2.15718, nFeature_RNA ≤ 2840.941, and nCount_RNA ≤ 5194.27. Typically, cells exhibiting excessively high total UMI counts and numbers of expressed genes were classified as doublets, whereas cells with elevated percentages of mitochondrial or ribosomal reads were considered to be of low quality, potentially undergoing apoptosis or fragmentation. Following these filtering steps, DoubletFinder (version 2.0.4) was employed to identify and remove doublets from each sample individually, completing the cell quality control process. Initially, data normalization was conducted using the normalizeData function. The cell cycle status was subsequently evaluated via the CellCycleScoring function, and highly variable genes were identified through the FindVariableFeatures method. The dataset was then scaled using ScaleData to standardize the data and to mitigate the influence of mitochondrial genes, ribosomal genes, and cell cycle effects on downstream analyses. Linear dimensionality reduction was performed using principal component analysis (PCA) through the RunPCA function, with significant principal components selected for further analysis. To address batch effects across different samples, the Harmony algorithm (version 1.1.0) was employed. This approach iteratively clusters similar cells from distinct batches within PCA space while maintaining batch diversity within clusters. Given the relatively mild batch effects observed in the dataset, default parameters (θ = 2) were applied. Nonlinear dimensionality reduction was subsequently carried out using RunUMAP, followed by the construction of a cell neighbourhood graph via FindNeighbors and cell clustering through FindClusters. For cell type annotation, a hierarchical annotation framework was implemented: primary manual annotation was based on characteristic gene expression patterns informed by the CellMarker database and pertinent literature; this was complemented by automated annotation results obtained from the SingleR software as a reference. To further increase the accuracy and comprehensiveness of cell type identification, multiple authoritative databases, including the Human Primary Cell Atlas (HPCA), BlueprintEncode, MonacoImmune, DatabaseImmuneCell, and NovershternHaematopoietic, were consulted. Cell annotation was conducted by querying the CellMarker (http://117.50.127.228/CellMarker/CellMarkerBrowse.jsp) database and reviewing the literature, aided by automated annotation support provided by SingleR (V2.4.0) software¹³. It aims to identify the cell types present in the corresponding tissue and their associated marker genes¹⁴.

Analysis of ligand - receptor interactions
In this study, CellCall (version 1.0.7) was used to perform a comprehensive analysis of intercellular communication networks¹⁵. By utilizing cell type annotations derived from Seurat alongside the raw count matrix, a normalized analysis object was constructed with parameters configured for the human genome. The TransCommuProfile function was applied to quantify the strength of cell-cell interactions through a weighted algorithm, implementing a significance threshold of a p value < 0.05 to identify reliable ligand-receptor pairs. The significant interaction pairs were subsequently subjected to KEGG pathway enrichment analysis via the getHyperPathway function, and the relationships between cell types and pathways were illustrated using bubble plots. The overall communication network was ultimately visualized through a circular plot, wherein eight distinct colours denoted different cell types. Interaction strength and directionality were represented by arrow features, providing a detailed characterization of intercellular signalling dynamics.

Pseudotime analysis
To investigate the dynamic transcriptional regulation of macrophages throughout the progression of heart failure, this study utilized the Monocle algorithm to conduct pseudotime analysis on macrophage subpopulations. The gene expression matrix corresponding to the target cell subpopulations was extracted to construct single-cell trajectory analysis objects, with highly variable genes selected as ordering features. By employing the DDRTree dimensionality reduction technique, cells were mapped onto a two-dimensional space to reconstruct the differentiation trajectory. Visualization analyses were performed to determine the distribution of cells along the pseudotime axis and to identify genes whose expression changed significantly over pseudotime. Subsequent analyses concentrated on the key gene DBNL and characterized its expression dynamics along the cell trajectory, elucidating the transcriptional reprogramming mechanisms of macrophages during heart failure progression¹⁶.

Gene set enrichment analysis (GSEA)
In this study, a GSEA approach was used to elucidate the regulatory mechanisms associated with key genes involved in heart failure. Utilizing previously identified key genes, samples were stratified into high and low-expression cohorts on the basis of the median expression value. Differential expression analysis was conducted employing the limma package, generating a ranked gene list according to log fold change (logFC). Subsequent KEGG pathway enrichment analysis was performed using the clusterProfiler tool, with gene sets sourced from the MsigDB database serving as the reference background. The GSEA algorithm was then applied to identify signalling pathways that were significantly enriched between the two expression groups, and an adjusted p-value threshold of less than 0.05 was used to determine statistical significance. To illustrate the regulatory functions of the core genes within critical pathways, various visualization techniques, including multipathway GSEA plots and circular network diagrams, were used.

Gene set variation analysis (GSVA)
GSVA is a nonparametric, unsupervised approach utilized for assessing gene set enrichment within transcriptomic data. This method transforms gene-level variations into pathway-level variations by calculating composite scores for specific gene sets, facilitating the evaluation of biological functional changes across various samples. In the present study, gene sets were sourced from the Molecular Signatures Database. The GSVA algorithm was employed to compute composite scores for each gene set, enabling the evaluation of potential biological functional alterations across different samples. The results of the GSVA enrichment analysis are provided in the supplementary material (Supplementary Table 1).

CTD drug prediction
The target gene (DBNL) was input into the search field of the Comparative Toxicogenomics Database (CTD), the disease category "cardiovascular disease" was selected, and the query was executed to retrieve drug prediction data associated with the condition "heart failure". The obtained prediction results were subsequently imported into Cytoscape software to facilitate data visualization and enable the construction of a gene-chemical interaction network map.

Molecular docking methods
Owing to the unresolved three-dimensional crystal structure of the human DBNL protein (UniProt ID: Q9UJU6), this study predicted the three-dimensional structure of DBNL on the basis of AlphaFold3¹⁷. Pirinixic acid (WY-14643) is available for download from the PubChem database (PubChem CID: 5694). Afterwards, the protein structure was preprocessed using MGLTools software (version 1.5.7)¹⁸, including steps such as the addition of hydrogen atoms. At the same time, proteins and small molecules were converted into the PDBQT format required for docking. AutoDock Vina software (version 1.1.2)¹⁹ was used for global molecular docking (exogenicity=16, num_modes=30) to explore potential binding modes. Upon completion of the docking calculations, the complex conformation with the highest affinity, indicated by the lowest binding free energy, should be selected as the initial structure for subsequent molecular dynamics simulations.

Molecular dynamics simulation method
To systematically investigate the binding stability and interaction mechanisms between candidate compounds and proteins, conventional molecular dynamics (MD) simulations were conducted utilizing the GROMACS software package (version 2024.03)²⁰. The protein parameters were generated using the Amber14SB force field²¹, the water molecule model was generated using the TIP3P model²², and the ligand topology parameters were generated utilizing the Antechamber Python Parser Interface (ACPYPE) tool, which is based on the General Amber Force Field (GAFF). The ligand-protein complex system was subsequently situated within a periodic boundary octahedral box filled with TIP3P water molecules. Sodium (Na⁺) and chloride (Cl⁻) ions were introduced to achieve a concentration of 0.15 mol/L and to neutralize the overall charge of the system. Upon completion of the system construction, the initial step involved minimizing the energy using the steepest descent method for 50,000 steps, with the aim of eliminating any potentially unreasonable conformations within the structure. Two stages of system equilibration were subsequently performed: a 100 ps NVT (constant particle number, volume, and temperature) simulation followed by a 100 ps NPT (constant particle number, pressure, and temperature) simulation. During these simulations, positional constraints were applied to the heavy atoms of the protein backbone to preserve the structural integrity of the protein. The temperature was maintained at 300 K using the V-rescale thermostat, and the pressure was regulated at 1 bar by employing the Parrinello-Rahman barostat for pressure coupling.Upon completion of the equilibration phase, a production phase simulation spanning 100 nanoseconds was conducted, during which all positional constraints are eliminated. The trajectory was integrated using a 2-femtosecond time step, and the particle-mesh Ewald (PME) method was employed to precisely manage long-range electrostatic interactions. The trajectory was saved every 10 ps, and a total of 10000 frames were output for subsequent analysis. In addition, stable trajectories in the 90-100 ns interval were extracted from the simulation, and the binding free energy of ligand protein complexes was calculated using the gmx MMPBSA tool²³.

Statistical analysis
The validity of this Mendelian randomization (MR) analysis is contingent upon three fundamental assumptions. (1) Relevance: The instrumental variables (IVs) must exhibit a strong association with the exposure. (2) Independence: The IVs must be independent of any confounders that affect both the exposure and the outcome. (3) Exclusion Restriction: The IVs should influence the outcome exclusively through their impact on the exposure. A breach of this assumption, where an IV affects the outcome through pathways not involving the exposure, is referred to as horizontal pleiotropy. All the statistical analyses were performed using R version 4.3.0, with two-sided tests, and a p value of less than 0.05 was generally considered to indicate statistical significance.

Results

Identification of DBNL as a gene associated with heart failure risk through Mendelian randomization and colocalization analysis
The outcome identifier for heart failure-related samples, GCST90162626, was acquired from the summary statistics dataset. Causal relationships among 1,247 pairs of expression quantitative trait loci (eQTLs) and outcomes were determined utilizing the extract_instruments and extract_outcome_data functions (Supplementary Table 2). Subsequent Mendelian randomiza...

Discussion

Although inflammation constitutes a fundamental aspect of the pathophysiology of heart failure (HF), previous investigations have predominantly concentrated on generalized inflammatory mediators and cell surface markers^17,18,24. This focus has constrained the elucidation of the specific cellular mechanisms involved. To overcome these limitations, we utilized an integrative approach combining multiomics Mendelian randomization (M...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
CellMarker Database	Harbin Medical University (CellMarker 2.0)	http://bio-bigdata.hrbmu.edu.cn/CellMarker/	CellMarker Database; For querying tissue/cell type marker genes during single-cell annotation
Comparative Toxicogenomics Database (CTD)	CTD (Duke University et al.)	http://ctdbase.org/	CTD Database; For retrieving interactions between chemicals, genes, phenotypes, and diseases, used for drug prediction
deCODE Plasma pQTL Dataset (2021, 4,907 aptamers, 35,559 European samples)	deCODE genetics	https://www.decode.com/summarydata/	2021 pQTL data release; Used as pQTL exposure data for MR analysis
eQTLGen Whole Blood eQTL Data	eQTLGen Consortium	https://www.eqtlgen.org/	eQTLGen phase I/II (specific version as per usage); Used as eQTL exposure data for MR and colocalization analysis
GEO Single-cell RNA-seq Dataset GSE161470 (Human Cardiac Tissue: 4 control + 1 pathological sample)	NCBI Gene Expression Omnibus (GEO)	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE161470	GSE161470; To obtain single-cell RNA sequencing data of human heart tissue for downstream analysis
GWAS Catalog	NHGRI–EBI	https://www.ebi.ac.uk/gwas/	GWAS Catalog; To retrieve HF-related GWAS information, annotations, and auxiliary analysis
Heart Failure GWAS Summary Statistics (GCST90162626)	NHGRI–EBI GWAS Catalog	https://www.ebi.ac.uk/gwas/studies/GCST90162626	Study accession: GCST90162626; Used as outcome data for MR and colocalization analysis
Molecular Signatures Database (MSigDB) v7.0	Broad Institute	https://www.gsea-msigdb.org/gsea/msigdb	MSigDB v7.0; For GSEA / GSVA background gene sets and pathway annotations
Pirinixic acid (WY-14643) Structure Data (PubChem)	NCBI PubChem	https://pubchem.ncbi.nlm.nih.gov/compound/5694	PubChem CID: 5694; Used to obtain ligand 3D structure for molecular docking and dynamics simulation
UniProt Protein Database (DBNL, UniProt ID: Q9UJU6)	UniProt Consortium	https://www.uniprot.org/uniprot/Q9UJU6	UniProt: Q9UJU6; Used to obtain the amino acid sequence of DBNL for structure prediction and docking
Whole Blood mQTL Meta-Analysis (Ref 11, 3,701 European samples, 426,636 mQTL traits)	Authors / Consortium (according to your reference 11)	Pending author to add original database URL / DOI	Reference DOI or Dataset ID; Used as mQTL exposure data for MR analysis