Inflammatory Biomarkers in Coronary Artery Disease: Insights From Mendelian Randomization and Transcriptomics

Introduction

Coronary artery disease (CAD), caused by atherosclerosis (AS), is a leading cause of mortality worldwide and presents a huge global economic burden.^1,2 Currently, noninvasive coronary evaluation using computed tomography (CT) and invasive coronary evaluation using angiography are the principal means of CAD diagnosis.³ Although percutaneous coronary intervention (PCI) and coronary artery bypass surgery (CABG) have undergone major advances in the last few decades, therapeutic agents, including antiplatelet agents, lipid-lowering medications, and β-blockers³ are still the baseline treatment. Although there are various prediction methods, such as risk factor evaluation and genetic variation,⁴ the prevalence of CAD continues to increase annually.⁵ Therefore, identifying novel biomarkers associated with CAD and developing therapeutic strategies may be useful in reducing the incidence and mortality of CAD.

Inflammatory response initiated by the damage to the arterial walls plays a key role in the pathophysiology of CAD,⁶ which results in an increased permeability to low-density lipoprotein (LDL), the formation of oxidized LDL and “foam cells”, accompanied by the secretion of inflammatory biomarkers such as IL-1β, IL-6 and tumor necrosis factor (TNF)-α,^7,8 with the activation of both pro-atherogenic and anti-atherogenic signaling pathways.⁹ Although few drugs that specifically block inflammatory cytokine pathways can reduce the risk of cardiovascular disease, the risk of infection has increased.⁷ Recently, immunomodulatory therapies targeting the immune system to reduce inflammation, such as cyclosporine and colchicine, have emerged. The blockage of other potential targets, such as the IL-6 pathway, may be beneficial in CAD.¹⁰ Elucidating the mechanisms of inflammation and immunity in CAD could provide novel insights into and intervention targets for CAD prevention and treatment.

Genetic risk factors substantially contribute to the development of CAD.¹¹ Genetic variants have been proposed as potential instrumental variables (IVs) to simulate the effects of modifiable environmental exposures on disease susceptibility, referred to as Mendelian randomization (MR).¹² MR, which resembles RCT, can be used as a potent tool to derive evidence for a direct causal relationship, overcoming the limitations of observational epidemiology prone to reverse causation. The advantages of MR include its capacity to establish causal links between exposures and outcomes, a task often complicated by confounding elements and biases in traditional observational studies. Additionally, by utilizing genetic variability as an instrumental variable, MR analysis reduces the influence of confounding variables that could distort the relationship between exposures and outcomes. Lastly, the insights gained from MR analysis can inform the development of targeted interventions or treatment strategies tailored to individual genetic profiles, thereby advancing the field of personalized healthcare.¹³ Over the past decade, epidemiological research has revealed numerous inflammatory biomarkers that are associated with an increased risk of CAD. However, the causality remains largely undefined. To date, only a few MR studies have provided convincing evidence that serum CRP level¹⁴ or fibrinogen¹⁵ has no causative role in CAD, but the IL-6R SNP (rs7529229 or rs2228145) is associated with an increased risk of CAD.¹⁶ Here, transcriptomics combined with MR was performed to provide a powerful framework to ascertain the causative role of inflammatory biomarkers in CAD. Differentially expressed genes (DEGs) were screened between the CAD and control groups in the dataset downloaded from the Gene Expression Omnibus (GEO) database. Differentially expressed inflammation-related genes (DE-IRGs) were identified by overlapping DEGs, key module genes, and inflammation-related genes (IRGs). Machine learning and MR analyses were performed to investigate the causal relationship between the DE-IRGs and CAD. Ultimately, 12 feature genes were acquired via four machine learning algorithms, and two key genes, hypoxia-inducible factor 1 alpha (HIF-1A) (odds ratio (OR) = 1.031, P = 0.024) and tumor necrosis factor-alpha-induced protein 3 (TNFAIP3, also known as A20) (OR = 1.104, P = 0.007), showed a causal association with CAD in the MR analysis. Activated NK cells and TNFAIP3 levels showed a moderately positive correlation (r = 0.52). Two drugs (USTEKINUMAB and METHOTREXATE) of TNFAIP3 and 136 drugs (BOLDINE and FLUPIRTINE, etc) of HIF1A were obtained from the DGIDB database. Finally, TNFAIP3 and HIF1A were significantly related to ribonucleic acid splicing, the T-cell receptor signaling pathway, and the NOD-like receptor signaling pathway. These results verified the causal relationship between HIF1A, TNFAIP3 and CAD, providing novel insights into drug R&D and potential therapeutic targets for CAD from the perspective of inflammation and immunity.

Materials and Methods

Data Source

Transcriptome data were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). There were 93 samples from CAD patients and 48 control peripheral blood samples from the GSE113079 dataset.¹⁷ The GSE12288 dataset includes 110 peripheral blood samples from patients with CAD and 112 control peripheral blood samples.¹⁸ A total of 935 IRGs were obtained through the search for “HALLMARK_INFLAMMATORY_ RESPONSE” and “GOBP_INFLAMMATORY_RESPONSE” in the Molecular Signatures Database (MSigDB) (https://www.pathwaycommons.org/).¹⁹

Differentially Expressed IRGs (DE-IRGs) Screening

Limma (version: 3.50.1) in the R package was used to screen DEGs between CAD and control samples in the GSE113079 dataset.²⁰ P adjust < 0.05 and |log2(fold change)| > 0.5 were regarded as filterable criteria. The ggplot2 (version 3.4.1)²¹ and ComplexHeatmap R packages (version 2.14.0)²² were used to draw volcano plots and heatmaps of the DEGs. Subsequently, to search for key module genes related to phenotypic traits (CAD and control groups), weighted gene co-expression network analysis (WGCNA) was performed using the GSE113079 dataset, based on the WGCNA R package (version: 1.71).^23,24 Specifically, key modular genes are genes that play important roles in specific biological processes or gene networks. These genes are usually involved in the core of gene regulation, signaling, or disease development, and serve as core nodes in gene co-expression networks. They have an important impact in maintaining cellular functions and systematic biological states, and are key factors in understanding biological mechanisms. All the samples were clustered to eliminate outliers. Based on the optimal power β and the scale-free evaluation factor, the samples in the GSE113079 dataset were divided into several modules using a Dynamic Programming Tree Cutting Algorithm (DTPTA), an algorithm commonly used for solving specific optimization problems, especially when searching for an optimal solution in a tree structure. The algorithm improves computational efficiency by decomposing a complex problem into smaller subproblems, solving each subproblem step-by-step, and avoiding repetitive computations by storing the results of the subproblems.Key module genes were obtained by computing the correlation coefficient matrix between the module eigenvectors and phenotypic traits. CAD-related genes were identified between DEGs and key module genes using the Venn diagram package in R (version 1.7.1). Finally, DE-IRGs were screened by intersecting the CAD-related genes and IRGs.

Functional Analyses of Biomarkers

To explore the biological pathways of the DEGs, a single-gene set enrichment analysis (GSEA) was performed based on the |logFC| sequences of all genes in the GSE113079 dataset. The reference gene sets were the KEGG and GO gene sets from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/). The screening criteria were adjusted to P < 0.05. The top five most significant functions enriched in Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were visualized using the R package enrichment plot (version: 1.18.0). To further characterize the biological functions and signaling pathways in which DE-IRGs were involved, the R package clusterProfiler (version: 4.2.2)²⁵ was used to perform GO and KEGG analyses. The screening criteria were adjusted to P < 0.05.

Establishment of Protein-Protein Interaction (PPI) Network and Analyses of Machine Learning

DE-IRGs were used to establish a PPI network (confidence > 0.4) using STRING Database (https://string-db.org/). Then, based on the cytoHubba plugin in Cytoscape software, the MCC algorithm was applied to identify the top 20 genes. Simultaneously, four machine learning methods were used to obtain feature genes. The glmnet R package (version 4.1-4)²⁵ was used to perform least absolute shrinkage and selection operator (LASSO) regression analysis. The caret R package (version: 6.0-93) (https://CRAN.R-project.org/package=caret) was used to execute the support vector machine‐recursive feature elimination (SVM-RFE) algorithm. The Boruta algorithm was applied using the Boruta R package (version 8.0.0 (https://CRAN). R-project. org/package=caret). The XGBoost algorithm was executed using the XGBoost R package (version 1.7.3.1) (https://CRAN.R-project.org/package=xgboost). Overlapping feature genes were acquired using LASSO, SVM-REF, Boruta and XGBoost. The GoSemSim (version: 2.24.0)²⁶ and corrplot R packages (version: 0.92) (https://github.com/taiyun/corrplot) were used to estimate the functional similarity and correlation of feature genes, respectively. The correlation coefficient (r-value) and P-value were computed for the feature genes based on Spearman correlation analysis.

MR Analyses

In order to investigate the causal relationship between CAD and characterized genes, two-sample Mendelian randomization was performed. The genome-wide association study (GWAS) data of “coronary artery disease” as the outcome and feature genes as exposure factors were obtained from the Integrative Epidemiology Unit (IEU) OpenGWAS database (https://gwas.mrcieu.ac.uk/) in the forward MR analysis. The feature genes were considered outcomes, and “coronary artery disease” acted as the exposure factor in the reverse MR analysis. The exposure factors were read and IVs were screened using the TwoSampleMR R package (version: 0.5.6).²⁷

In Mendelian randomization (MR), SNPs were used as instrumental variables to estimate the causal effect of an exposure on an outcome. The purpose of using SNPs in MR are to provide a way to infer causality between an exposure and an outcome that are less prone to confounding and reverse causation compared to observational studies. By leveraging the genetic variation that naturally occurred in populations, MR aimed to mimic the random assignment of treatments in a randomized controlled trial, thus providing a stronger basis for establishing a causal relationship. SNPs of the exposure factors were chosen as candidate IVs for further MR analysis. The SNPs were singled out with a genome-wide significance level of P < 5×10⁻⁸ and IVs with linkage disequilibrium were excluded with the criteria of R² = 0.001 and kb = 10 in the forward MR analysis. SNPs with a genome-wide significance level of P < 5×10⁻⁶ and IVs of linkage disequilibrium were excluded with the criteria of R² = 0.001 and kb = 10000) in the reverse MR analysis.

The MR study met the following hypotheses: (i) IVs are closely related to the exposure factors. (ii) IVs are not associated with confounding factors. (iii) IVs influence the outcome only through exposure factors. Five algorithms–inverse variance weighted (IVW),²⁸ MR Egger,²⁹ weighted median (WM),³⁰ simple mode,²⁷ and weighted mode³¹ were applied to the bidirectional MR analysis. Sensitivity analyses (Cochran’s Q statistical test, MR-Egger test, and leave-one-out analysis) were performed to assess the reliability of MR analysis. The exposure factors that significantly affected CAD were identified as the key genes.

Building of Alignment Diagram Model Based on Key Genes

To accurately predict the occurrence of CAD, a logistic model was built in the GSE113079 dataset based on the expression levels of the key genes. The rms R package (version 6.5-0) was used to construct an alignment diagram of the key genes in the GSE113079 dataset. The predictive power of the model was determined using calibration and decision curve analysis (DCA).

Correlation Analysis of Key Genes and Interleukin 6 (IL6)

IL6 was associated with the cause and mortality of CAD,³² and the correlation between key genes and IL6 was accessed using the corrplot package. PPI networks interacting with HIF1A, TNFAIP3, and IL6 were obtained from the STRING database. In addition to the protein interaction network, the biological pathways of key genes and IL6 are understand. The samples in the GSE113079 dataset were divided into high and low expression groups based on key genes, and single-gene GSEA was performed to investigate the biological pathways related to key genes and IL6.

Analysis of the Immune-Infiltrating Cells and Prediction of Small-Molecule Drug

CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of known RNA Transcripts) was an algorithm designed to estimate the relative abundance of cell types in complex tissues by parsing mixed gene expression data. Its basic principle involved using linear support vector regression (SVR) to decompose the mixed expression data and infer the contribution of individual cell types. Typically, CIBERSORT utilized the LM22 gene set, which comprised a set of gene expression signature matrices that define 22 immune cell types. A linear support vector regression model was trained using the SIGNATURE gene matrix as the feature matrix and the mixed gene expression matrix as the response variable. The relative abundance of each cell type in each sample was predicted by the regression model. Signature gene matrices were predefined to include the genes that characterized each cell type and their expression levels in that cell type. For example, the LM22 gene set defined 22 immune cell types. The CIBERSORT algorithm was used to determine the proportion of the 22 immune-infiltrating cells in the GSE113079 dataset. In addition, the abundance of immune-infiltrating cells in CAD and control samples was calculated using the Wilcoxon test. The R package psych (version: 2.2.9) (https://CRAN.R-project.org/package=psych) was used to analyze the Spearman correlation for both differential immune-infiltrating cells and key genes.

Small molecule drug prediction aimed to identify and prioritize potential therapeutics by evaluating their likelihood to interact effectively with specific biological targets, thereby achieving the desired pharmacological effects. This approach expedited the initial phases of drug development by minimizing the requirement for extensive lab trials and enabling researchers to concentrate on the most promising candidates. This study, the DGidb database (https://dgidb.org/) was used to predict small-molecule drugs corresponding to key genes and IL6. Hub gene and small molecule drug correspondences were imported into the Cytoscape software to visualize their relationships.

Validation of Key Genes Expression

The Wilcoxon test was applied to compute key genes and IL6 expression between CAD and control samples in the GSE113079 and GSE12288 datasets. In addition, whole blood samples were obtained from five CAD and five control individuals. Total RNA was extracted using the TRIzol reagent (Ambion, China). The mRNA was reverse-transcribed using a SureScript-First-strand-cDNA-synthesis-kit (Servicebio, China) and cDNA was collected. Real-time fluorescence quantitative PCR (RT-qPCR) was performed using the SYBR Green qPCR Master Mix. The relative gene expression was calculated using the 2^−ΔΔCt method. Primer sequences are shown (Table S1).

Statistical Analysis

All analyses were performed using the R package. The Wilcoxon test was used for group comparisons, except for the special instructions. *P < 0.05, considered statistically significant.

Results

Identifying 92 DE-IRGs in CAD and Control Samples

First, 3722 DEGs were identified in GSE113079. Among of which, 1772 genes were upregulated and 1950 genes were downregulated in CAD group than that in control group. The top 10 upregulated and downregulated genes are marked in the volcano plot, and the expression of the marked genes is shown in the heatmap (Figure 1A and B). In addition, based on the optimal power β of nine and scale-free evaluation factor R2 = 0.85, a scale-free network was built. According to the correlation coefficient between genes, genes with similar expression patterns are classified into the same gene module, and different gene modules are distinguished by different colors. These genes were divided into 13 co-expression modules by WGCNA analysis (Figure 1C–F). The yellow module with the highest positive correlation with CAD (r = 0.79, P < 0.0001) and brown module with the highest negative correlation (r =−0.77, P < 0.0001) were selected as key modules (Figure 1G). A total of 1839 common DEGs (1108 upregulated genes and 731 downregulated genes) were identified among the DEGs and key module genes (Figure 1H and I). Finally, 92 DE-IRGs were identified between common DEGs and IRGs (Figure 1J).

Figure 1 DE-IRGs identified in CAD. (A) Volcano plot of the distribution of differentially expressed genes between CAD and Control. (B) Heatmap of differentially expressed genes between CAD and Control. (C) Sample-level clustering with pre-sample CAD and Control removed. (D) Sample-level clustering after introduction of sample traits. (E) soft threshold filtering. (F) WGCNA analysis to identify co-expression modules. (G) Heatmap of correlation between modules and phenotypes. (H) Wayne plots of up-regulated DEGs and yellow module genes. (I) Wayne diagram of DEGs and brown module genes. (J) Wayne diagram of intersection of CAD-related and inflammation-related genes.

A Total of 12 Feature Genes of CAD Were Acquired

To identify the feature genes, a PPI network was constructed based on 92 DE-IRGs (Figure 2A). The top 20 genes (IL18R1, CX3CL1, ICAM1, and others) were screened using the MCC algorithm (Figure 2B and Table S2). To gain further feature genes, four machine learning methods (LASSO, SVM-RFE, Boruta, and XGBoost) were applied based on the top 20 genes. Twelve common feature genes, including NFKB1, TNFAIP3, and HIF1A, were acquired using four algorithms (Figures 2C and S1). The CCL4 and TNFAIP3 were the highest positive correlation (r = 0.82, P < 0.001), whereas CCL24 and NFKB1 showed the highest negative correlation (r = −0.69, P < 0.001) (Figure 2E and F).

Figure 2 Feature genes of CAD acquired by machine network learning. (A) PPI network was constructed based on the 92 DE-IRGs. (B) Gene PPI interaction network of top20. (C) 12 feature genes were acquired by four methods of machine learning. (D) Functional similarity estimated of these feature genes. (E) Heatmap of correlations between characterised genes. (F) Scatterplot of correlation between character genes.

HIF1A and TNFAIP3 as Risk Factors for CAD

The causal relationship between the feature genes and CAD was explored using MR analysis. We initially screened the 12 characterized genes for SNPs. Subsequently, we found that SNPs were present in the CCR6, NFKB1, XCR1, HIF1A, and TNFAIP3 genes (Tables S3–S7). In MR analysis, two feature genes (HIF1A and TNFAIP3) showed a causal association with CAD, and we chose them for follow-up analysis. The results of univariate MR analysis indicated that HIF1A (OR_IVW = 1.031, P_IVW = 0.024; OR_WM = 1.041, P_WM = 0.039) and TNFAIP3 (OR_IVW = 1.104, P_IVW = 0.007; OR_WM = 1.141, P_WM = 0.006) were risk factors for CAD (Tables 1 and 2). The scatter plot and forest plot of the IVW analysis also confirmed the above result (Figure 3A–D). In addition, CAD was not related to HIF1A or TNFAIP3 in the reverse MR analysis (Tables S8–S10 and Figure S2).

Table 1 The Association of HIF1A With CAD by Univariate MR Analysis

Table 2 The Association of TNFAIP3 With CAD by Univariate MR Analysis

Figure 3 SNPs related to risk factors for CAD. (A) Scatterplot of forward MR analysing the effect of instrumental variables (SNPs) on exposure factors (HIF1A) and the effect of instrumental variables (SNPs) on outcome (Coronary artery disease). (B) Forest plot of effect sizes of exposure factors on outcome variables analysed by IVW for positive MR. (C) Scatterplot of the effect of TNFAIP3-Coronary artery disease positive MR analysis instrumental variable (SNP) on exposure factor (TNFAIP3). (D) Positive MR analysis of TNFAIP3-Coronary artery disease Exposure factors on junction Forest plot of effects of local variables (E) Funnel plot of instrumental variables for positive MR. (F) Scatterplot of instrumental variables for positive MR of TNFAIP3-coronary artery disease.

The MR Analysis Results Were Reliable

Sensitivity analysis was performed to assess the reliability of the MR analysis results. There was no heterogeneity and pleiotropy (P > 0.05) among the samples in the Cochran’s Q statistical test (Tables 3 and 4) and the MR-Egger test (Tables 5 and 6). The selection of SNPs met the requirements of exclusivity and independence in the leave-one-out analyses (Figure 4A–D). In summary, all MR analyses revealed that HIF1A and TNFAIP3 were risk factors for CAD.

Table 3 Sensitivity Analysis of MR Analysis About the Association of HIF1A With CAD

Table 4 Sensitivity Analysis of MR Analysis About the Association of TNFAIP3 With CAD

Table 5 Pleiotropy Assessment of MR Analysis About the Association of HIF1A With CAD

Table 6 Pleiotropy Assessment of MR Analysis About the Association of TNFAIP3 With CAD

Figure 4 The exclusivity and Independence assumptions of the selection of SNPs in the leave-one-out analyses. (A and C) Forward and Backward MR’s Reject-by-Removal Test Forest Plots. (B and D) Forward and reverse culling test forests in TNFAIP3-Coronary artery disease.

Analyses of HIF1A and TNFAIP3 for CAD Prediction

A logistic regression model was generated based on the expression of HIF1A and TNFAIP3. The logistic model predicted the occurrence of CAD with 81% (Figure 5A–C). The ROC curves of HIF1A (area under the curve (AUC)=0.894) and TNFAIP3 (AUC=0.816) showed that they could accurately predict CAD alone (Figure 5D). Subsequently, an Alignment Diagram was constructed, which showed better consistency for CAD prediction using the calibration curve (P = 0.5) (Figure 5E and F). In addition, DCA indicated that the alignment diagram was effective in clinical practice (Figure 5G).

Figure 5 The predictive value of HIF1A and TNFAIP3 for CAD. (A) Confusion matrix plot for the training set. (B) ROC plot for the training set. (C) ROC plot for the validation set. (D) Individual ROC curves for the hub gene plotted in the training set (GSE113079). (E) hub Gene Construction Columnar Diagram Model. (F) Calibration curves for line diagram models. (G) Decision curve for DCA.

Positive Correlation Between IL6 and TNFAIP3

The IL6 and TNFAIP3 levels were positively correlated (r = 0.5, P < 0.001) (Figure 6A and B). Furthermore, a consistent trend was observed in the GSE113079 and GSE12288 datasets. TNFAIP3, HIF1A, and IL6 were expressed at lower levels in patients with CAD than in controls (Figure 6C–G). There were 11 genes (ARNT, ARNT2, etc) related to HIF1A, ten genes (IKBKG, TNFRSF1A, etc) related to TNFAIP3, and ten genes (IL10, TNFRSF1A, etc) related to IL6 in the PPI network (Figure 7A–C). TNFAIP3, HIF1A, and IL6 were significantly related to the NOD-like receptor signaling pathway (Figure 7D–F).

Figure 6 Correlation of IL6 with HIF1A and TNFAIP3. (A and B) Scatterplot of correlation between IL6 and hub genes. (C–E) The expression of TNFAIP3, HIF1A and IL6 in CAD patients and control individuals, with red representing the CAD group and blue representing the Control group. *p < 0.05; **, p < 0.01; **** p < 0.0001 compared with the control group. (F and G) The results of RT-qPCR showed a lower expression of TNFAIP3 and HIF1A in CAD group.

Figure 7 Individual gene PPI analysis and hub gene GSEA enrichment analysis. (A–C) PPI networks indicated genes related to TNFAIP3, HIF1A and IL6. (D–F) Single-gene GSEA analysis investigated the biological pathways related to TNFAIP3, HIF1A and IL6.

Positive Correlation Between Activated NK Cells and TNFAIP3

The abundance of 22 immune-infiltrating cells in all samples was determined. There were higher cell compositions of M0 Macrophages, Monocytes, naive CD4⁺ T cells, and regulatory T cells in the CAD samples than in the control samples (Figure 8A and B). TNFAIP3 and NK cells had a moderate positive correlation (r = 0.52, P < 0.001), and HIF1A (r = −0.40, P < 0.001) and TNFAIP3 (r = −0.47, P < 0.001) were negatively correlated with naive CD4⁺ T cells (Figure 8C). These results may enhance our understanding of the role of the CAD inflammation-associated hub gene in TME formation. Finally, in the gene-drug correspondence network, 79 Nodes (3 genes, 76 drugs) and 76 Edges were included, of which 1 drug was predicted for TNFAIP3, 48 drugs were predicted for IL6, and 144 drugs were predicted for HIF1A. All drugs were approved by FDA (Figure 8D).

Figure 8 Analysis of immune-infiltrating cells and prediction of small-molecule drug. (A) The abundance of 22 immune-infiltrating cells and different infiltration of immune cells between CAD samples and control samples. (B) A box plot of differential immune infiltrating cell abundance was constructed based on CAD and Control. Red represents the CAD group, and blue represents the Control group. ** p < 0.01; **** p < 0.0001 compared with the control group. (C) Lollipop plot of correlation between differentially immune infiltrating cells and hub gene, IL6. (D) Small molecule drugs of TNFAIP3, HIF1A and IL6 were forecasted from DGidb database.

Discussion

The basic pathogenesis of CAD and coronary atherosclerosis is a complex, persistent, and progressive inflammatory process involving interleukins (ILs), such as IL-6.³³ The process was initiated with the retention and oxidation of low-density lipoprotein (LDL) inside the intima, followed by the dysfunction of endothelial cells, the formation of the “foam cells”, the proliferation and migration of vascular smooth muscle cells (VSMCs), and ultimately the formation of atherosclerotic plaques.³⁴ When LDL cholesterol is controlled by statins and other lipid-lowering drugs, the risk of major adverse cardiovascular events (MACE) is believed to be inflammatory in nature.³⁵ Therefore, there is a need to develop anti-inflammatory intervention strategies that target specific inflammatory mediators to reduce the risk of cardiovascular disease.³⁶ In the present study, we investigated the genes that contribute to the inflammatory response to CAD. By performing a bidirectional Mendelian randomization analysis, two key genes, HIF1A and TNFAIP3, were confirmed to be risk factors for CAD. However, CAD was not related to HIF1A and TNFAIP3, providing a potential reference for the treatment and prevention of CAD.

In the MR analysis, we selected only six feature genes from the 12 because the SNPs related to the other six genes were too few to meet the analysis requirements. Among the six feature genes, P-values of HIF1A, NFKB1, CCR6, TNFAIP3, and XCR1 were lower than 0.05, indicating that they were significantly associated with CAD. However, sensitivity tests for NFKB1, CCR6, and XCR1 indicated confounding factors; therefore, their relationship with CAD is unreliable. Therefore, we chose HIF1A and TNFAIP3 for further analysis.

HIF-1A is an oxygen-sensitive transcription factor that mediates adaptive metabolic responses to hypoxia.³⁷ Previous studies have demonstrated that HIF-1A promotes mitochondrial damage, neuronal apoptosis, and expression of inflammatory factors.³⁸ Individuals with a high risk of atherosclerosis show significantly high level of HIF-1A expression.³⁹ However, the causal relationship between HIF-1A and CAD remains unclear. Here, we report, for the first time, that HIF1A is a risk factor for CAD and can predict CAD occurrence with an accuracy of 89.4%.

TNFAIP3 is a ubiquitin-modifying protein known to preserve immune homeostasis and prevent autoimmune diseases by negatively regulating NF-κB signaling.⁴⁰ In this study, the mRNA expression of TNFAIP3 was lower in patients with CAD than in control individuals. Consistent with our results, overexpression of TNFAIP3 attenuates atherosclerosis progression in both hyperlipidemic mouse models and humans.^41,42 Indeed, TNFAIP3 exerts numerous anti-atherogenic functions by reducing endothelial cell dysfunction, limiting immune cell infiltration, and inhibiting smooth muscle cell migration and proliferation.⁴³ Genome-wide association studies have demonstrated that SNPs at the TNFAIP3 locus are risk factors for coronary artery disease,⁴⁴ and mutations in TNFAIP3 are associated with a variety of inflammatory diseases.⁴⁵ Consistent with these reports, it was also found that TNFAIP3 had a causal association with CAD in our study and could predict CAD occurrence with an accuracy of 81.6%.

Much evidence has indicated an inseparable connection between the immune microenvironment and CAD.^31,46 Our analysis showed that the relative percentages of M0 Macrophages, Monocytes, CD4+ naive T cells, and regulatory T cells were higher in the CAD group, whereas the proportions of activated CD4+ memory T cells, activated dendritic cells, activated NK cells, and CD8+ T cells were lower than those in the normal group. Damaged endothelial cells release monocyte chemotactic protein-1, which recruits monocytes to lesions and then differentiates into macrophages, leading to inflammation and atherosclerotic plaque development.⁴⁷ In addition, HIF1A and TNFAIP3 levels negatively correlated with naive CD4⁺ T cells. GO and KEGG analyses of these targeted genes revealed that they were enriched in the T-cell receptor signaling pathway, suggesting that they may be involved in CAD via T cells. A functional and homeostatic defect in regulatory T cells is associated with a subclinical pro-inflammatory and atherogenic state, and decreased CD4+ naive T cell counts are independently associated with cardiovascular events.^48,49 We believe that there is a need for immune restoration strategies for T cells in CAD patients.

As the core of the inflammatory cascade, IL-6 plays a pivotal role in the initiation and progression of atherosclerosis.^50,51 In recent years, emerging anti-inflammatory approaches targeting the IL-6 pathway have been shown to significantly reduce the incidence of cardiovascular events.^10,52 In our study, TNFAIP3 and HIF1A were positively correlated with the gene expression IL-6 according to the correlation analysis. In addition, they were enriched in the NOD-like and Toll-like receptor signaling pathways, both of which are related to the innate immune response. NODs and TLRs recognize pathogen-associated molecular patterns (PAMPs), followed by the activation of NF-κB and MAPK signaling mediated by myeloid differentiation factor 88 (MyD88) or non-MyD88, leading to the secretion of pro-inflammatory cytokines and eventually the occurrence and development of atherosclerosis.^53,54

Finally, we analyzed the key genes for gene-targeted drugs and two small-molecule drugs (USTEKINUMAB and METHOTREXATE) of TNFAIP3, 136 drugs (BOLDINE, FLUPIRTINE, etc) of HIF1A, and 25 drugs (INSULIN, COR-001, etc) of IL6 predicted from the DGidb database. Among these, METHOTREXATE, a traditional immunomodulatory and anti-inflammatory drug, has been proved to have beneficial effects on atherosclerosis and cardiovascular clinical endpoints in some experimental and clinical studies.^55,56 However, more randomized controlled studies are needed to confirm the therapeutic potential of METHOTREXATE in cardiovascular prevention.⁵⁷ USTEKINUMAB and COR-001, monoclonal antibodies against anti-IL-12, IL-23, and IL-6, have been used in skin and kidney-related diseases,^58,59 and their application in cardiovascular diseases based on evidence-based medicine needs to be confirmed. BOLDINE is a plant-derived bioactive compound with beneficial effects on human health. For example, boldine protects endothelial cells from hypertension by increasing vascular NO production and reducing ROS overproduction.⁶⁰ However, its application to CAD has not yet been reported. Further studies are required to elucidate the role of these gene targeting drugs in CAD.

Strengths and Limitations

The strength of this study is the causal relationship between inflammatory markers and coronary heart disease. Despite these meaningful findings, this study has several limitations. First, it was based on publicly available data, and additional external cohorts are required to validate our findings. In future studies, we will conduct experiments using cellular models or animal models, etc., to further confirm our findings and to delve into the specific mechanism of action of the relevant genes in the development of CAD. Second, gene expression may not be directly equivalent to protein expression, and continued attention should be paid to the roles of these potential key genes in regulating the inflammatory immune network in coronary heart disease.

Conclusions

In this study, we identified two key inflammation-related genes (HIF1A and TNFAIP3) involved in CAD, and proved that these two hub genes could predict the occurrence of CAD with high accuracy. The study presents novel and impactful findings that contribute to the understanding of CAD’s genetic and inflammatory mechanisms. By integrating Mendelian randomization, transcriptomics, and machine learning, it advances the field of CAD inflammation-immune research, particularly in identifying novel biomarkers for early detection and therapeutic targeting of CAD.