Novel Diagnostic Biomarkers Related to Necroptosis and Immune Infiltration in Coronary Heart Disease

Introduction

Coronary heart disease (CHD) remains a major public health threat globally, serving as a major cause of death and contributing to rising healthcare costs.¹ The main pathological basis of CHD is atherosclerosis (AS), the buildup of subendothelial lipid lesions, and the accompanying inflammatory response.^2,3 Atherosclerotic plaques usually accumulate silently, maintaining a subclinical atherosclerotic status for decades; during this time, the condition may remain undetected in patients due to a lack of typical symptoms.⁴ However, acute coronary syndrome (ACS), which occurs when a plaque becomes eroded or ruptures, can rapidly become fatal. Despite significant advances in CHD management, there have been few improvements in the clinical situation and prognosis of CHD. The DISCHARGE study showed that there is no effective examination that could improve the prognosis of CHD.⁵ Currently, no specific and sensitive biological markers are available for the diagnosis of CHD until ACS develops; therefore, there is an urgent need for validated peripheral blood markers.

AS progression is marked by diverse forms of cell death. Among these, necrosis is a major mode of cell death in the microstructure of human atherosclerotic plaques.^6,7 Necrosis is usually thought to occur unexpectedly and passively. However, this view has been challenged by the identification of necroptosis, which is initiated by a meticulously regulated signaling cascade. Necroptosis is induced through a strictly regulated pathway mainly in which cell death ligands (eg tumor necrosis factor, Toll-like receptor ligands, etc) bind to receptors, activate receptor-interacting protein kinase 1 (RIPK1), recruit receptor-interacting protein kinase 3 (RIPK3) and phosphorylated mixed lineage kinase-like (MLKL) after forming the RIPK1/RIPK3 complex, and oligomerize phosphorylated MLKL to form necrotic vesicles that mediate cell death.^8–11 Emerging research has indicated that necroptosis plays a role in cardiac pathophysiology, including conditions such as AS, ischemic reperfusion injury, and myocardial infarction.^12–14 Karunakaran et al found that RIPK3, MLKL, and p-MLKL levels were significantly increased in human AS plaques and that administering programmed necrosis inhibitor, Nec-1, in animal experiments yielded a significant decrease in atherosclerotic plaques, as well as an increase in plaque stability.¹⁵ Therefore, necroptosis plays a major role in CHD.

The evolution of AS is essentially a battle between anti-inflammatory and pro-inflammatory processes dominated by immune cells, and AS progression is considered to be correlated with immune cell infiltration in the vascular wall, including the location, quantity, and function of immune cells.¹⁶ Necroptosis is recognized as a trigger of inflammation¹⁷ and can induce a proinflammatory state.¹⁵ Therefore, our study sought to ascertain whether necroptosis genes can serve as biological markers of CHD and whether they play a role in CHD through immune inflammation or infiltration.

Herein, we first established the landscape of necroptosis genes in CHD and identified the hub genes using machine learning algorithms. A diagnostic model was further constructed upon the hub genes and its predictive efficacy was validated. Following this, clinical samples and animal tissues were used to confirm the expression levels of the hub genes. Finally, this study explored the influence of these necroptosis genes on CHD and their effect on immune cell infiltration levels, which may offer a novel diagnostic approach for CHD.

Materials and Methods

Data and Pre-Processing

Our study proceeded according to the process shown in Figure 1. GSE20680¹⁸ and GSE20681,¹⁸ which constitute the gene expression data (Table S1) for the CHD and control groups, respectively, were obtained from the Gene Expression Omnibus database. Both datasets had human origins, and both used the GPL4133 microarray platform. GSE20680 consists of 139 samples, which includes 52 control samples and 87 samples from patients with CHD. GSE20681 consists of 198 samples, including 99 control samples and 99 samples from patients with CHD. All samples were obtained from whole blood cells isolated from patients undergoing cardiac catheterization. The two datasets were merged for further research, and the R package sva¹⁹ was employed to address the issue of batch impact across various datasets and to conduct log₂ normalization.

Figure 1 Workflow.

Constructing a Landscape of Necroptosis Genes

First, 17 necroptosis genes were acquired from the literature²⁰ and their expression was examined in all samples. The specific gene names are shown in Table S2. Analysis of the gaps between the two groups was performed using the limma package.²¹ The ggplot2 and complex heat map packages were used to plot volcano maps and expression heat maps of the above genes for all samples. For the CHD and control groups, box plots and violin combination plots were produced using the R package ggpubr.²² The locations of the 17 genes in the chromosomes were mapped using the R package Rcircos,²³ which also provided the chromosome data. Downloads from the ENSEMBL²⁴ database were used to determine the loci of the genes on the chromosomes.

Correlation Analysis Between Necroptosis Genes

To explore the associations between necroptosis gene expression in each patient, Spearman correlation coefficient was computed between the genes. Correlation scatter plots between pairs of genes that satisfied the specified criteria were plotted, and correlation curves were fitted using the R package ggplot2. Graph-edge histograms and density curves were obtained using the R package ggExtra.²⁵

Constructing a Diagnostic Model Based on Necroptosis Genes

As necroptosis genes play vital biological roles in the progression of CHD, constructing a diagnostic model based on them is possible. In this study, hub genes were selected due to their significant changes in differential expression analyses across different datasets, confirmed by Venn diagram intersection analysis. First, we screened 17 necroptosis genes using least absolute shrinkage and selection operator (LASSO), implemented the method using the R package glmnet,²⁶ and selected the best lambda value. Genes with regression coefficients besides zero were retained. Subsequently, the top 10 important genes were screened using a random forest method.²⁷ Finally, intersection analysis of the genes selected for LASSO regression, the top 10 important genes screened by the random forest method and differentially expressed necroptosis genes between patients with CHD and controls. Four hub genes were identified from the overlapping regions of the three gene sets in the Venn diagram to construct disease diagnostic models, which were visualized using a nomogram.²⁸ The anticipated precision of the diagnostic model was validated by plotting the calibration and decision curve analysis (DCA), where calibration curves were plotted by the R package rms²⁹ and DCA by the ggDCA³⁰ R package.

Constructing Protein Interaction Networks

As genes are closely interlinked, particularly those responsible for common biological activities, a protein-protein interaction (PPI) network was constructed to determine the relationships between the necroptosis-related genes.

The PPI network was created by inputting the mentioned genes with the confidence level set to a default value of 0.4 and utilizing the STRING³¹ database. Subsequently, the PPI was analyzed in depth using Cytoscape³² software to count each node’s network attributes and identify hub nodes according to the degrees of the nodes (hub nodes were defined as those with 10 nodes having a top-10 rating) using the plugin Cytohubba.³³ These nodes, which exhibit a strikingly high degree of connection, merit additional study because they could be crucial in the control of the entire biological process; therefore, we conducted further prediction studies on 10 hub nodes based on the miRNet2.0 database (https://www.mirnet.ca/miRNet/home.xhtml)³⁴ for miRNAs and transcription factors of hub nodes, respectively. Cytoscape was used to export, analyze, and visualize the prediction findings.

Single-Sample Gene Set Enrichment Analysis (ssGSEA) Scores of the Samples

A modified variant of GSEA, ssGSEA, determines an enrichment score that measures the level of enrichment of a gene set in a single sample within a given data set³⁵ Since the samples in the CHD group had different amounts of necroptosis gene expression and heterogeneity was prevalent among the patients, the samples were scored with ssGSEA based on 17 necrotic genes to explain the heterogeneity between individual samples in the CHD group using the ssGSEA score.

Seventeen necroptosis genes were used as background genes to assess each sample in the CHD group using the GSVA package. The ssGSEA median score was employed as the cutoff value to separate the CHD group into subgroups with high and low scores. The ggpubr package was used to visualize variations in the expression of each gene between the two groups.

Analysis of Expression Differences and Functional Enrichment

To further reveal the heterogeneity within the CHD groups, both functional enrichment and differential expression analyses were performed for the high- and low-score groups, and heat and volcano maps were used to display the findings of the differential expression analysis. Genes were used for GSEA when they had adjusted p < 0.05, and absolute log₂fold change (FC) values >0.5, which indicated that they were significantly differentially expressed. Gene Ontology (GO) enrichment analysis is used extensively in functional enrichment studies of genes in numerous dimensions and at various levels. The biological process (BP), molecular function (MF), and cellular component (CC) levels are typically used.³⁶ Functional annotations of all genes that exhibited substantial differential expression was performed using the R package Clusterprofiler^37,38 to identify statistically significant biological processes. Enrichment results were selected to visualize the top five in each group via bubble and circle plots, and an adjusted p < 0.05, which was the cutoff point for determining the significance of the enrichment analysis.

GSEA is an automated approach to assess the statistical disparities between two biological states in a given set of predetermined genes. It is frequently used to evaluate alterations in pathway activity and biological processes within expression datasets obtained from samples.³⁹ The gene expression profiling dataset was used to investigate how the two groups differed in terms of biological processes. The reference gene set “c2.cp. v7.2. symbols.gmt” was retrieved from the MSigDB database⁴⁰ and then enhanced and displayed using the GSEA function within the R package ClusterProfiler. Statistical significance was set at p < 0.05. Additionally, differences in gene expression patterns between the two groups were visualized.

Immune Infiltration Analysis

The immunological microenvironment is made up of an interwoven network of interstitial tissues, immune cells, inflammatory cells, and various cytokines and chemokines. Analysis of immune cell infiltration is critical in fields such as disease research and treatment. The amount of immune cell infiltration between the two groups was explored using ssGSEA with the R-package GSVA.³⁵ The signature genes for 28 immune cells were sourced from the literature.⁴¹ This background set of genes was analyzed individually for ssGSEA. The results of immune cell infiltration are shown using heat maps and box plots. The associations between immune cells in various illness states were plotted using the R package corrplot.⁴²

To guarantee the validity of the results, the degree of immune cell infiltration was also measured using the R package CIBERSORT.⁴³ The number of immune cells per sample was determined using the LM22 background gene set provided by CIBERSORT, thereby indicating the extent of infiltration, and the findings were displayed using box plots. Furthermore, lollipop plots were generated with the R program ggplot2 to show the associations between hub genes and immune cells.

In vivo Hub Gene Validation

Patient Selection and Clinical Procedures

We enrolled 102 patients who underwent coronary angiography from January 2022 to December 2022 at the Yan’an Hospital in Kunming. Inclusion criteria for the CHD group were as follows: maximum stenosis ≥ 70% in the primary epicardial coronary artery or ≥ 50% in two arteries by invasive angiography. Inclusion criteria for the control group were as follows: patients with 0% stenosis in the coronary arteries by invasive angiography. Exclusion criteria were as follows: valvular disease, nonischemic cardiomyopathy, severe hepatorenal insufficiency, neoplastic disease, systemic infections, or inflammatory conditions. The inclusion and exclusion criteria for this study were based on the corresponding criteria for clinical studies from the GSE20680 and GSE20681 dataset sources.^44,45 This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Blood Samples

Peripheral blood (5 mL) was used to extract peripheral blood mononuclear cells (PBMCs) using Lymphatic Separation Solution (STEMCELL Technologies, Vancouver, BC, Canada) in accordance with the manufacturer’s instructions. Subsequently, 1 mL of TRIzol reagent (Invitrogen, Waltham, MA, USA) was added to the cell precipitate, which was fully dissolved and stored at −80 °C. All participants completed signed informed consent forms after the study received approval from Medical Ethics Committee of Kunming Yan’an Hospital (No. 2022–096-01).

RNA Extraction and Quantitative PCR (qPCR)

The TRIzol method was used for total RNA extraction. A NanoDrop2000 instrument (Thermo Fisher Scientific, Waltham, MA, USA) was used to determine the quantity and quality of the RNA samples. The 5× All-In-One MasterMix (AccuRT Genomic DNA Removal Kit; Applied Biological Materials, Vancouver, BC, Canada) was used for reverse transcription. mRNA expression was detected using the BlasTaq 2× qPCR MasterMix (AppliedBiological Materials). The relative mRNA expression normalized to glyceraldehyde-3-phosphate dehydrogenase was measured using the 2^−ΔΔCt method. The primers used are listed in Table S3.

AS Mouse Model Construction

Twelve eight-week-old mice (Viton Lever, Beijing, China), comprising six APOE^−/− mice and six C57BL/6J wild-type mice were divided into AS model and control groups, respectively. For 12 weeks in a row, the AS model group was fed a high-fat diet, whereas the control group received a regular diet. At the end of the 12^th week, the mice were executed via cervical dislocation, and aortic samples were collected. Laboratory animal studies were conducted under a license granted (2022084) by the Animal Ethics and Welfare Committee (AFWC) of Yan’an Hospital Affiliated to Kunming Medical University, according to its standards for the care and use of animals.

Immunohistochemistry

Fresh mouse aortic root specimens were paraffin-embedded and sectioned (4 μm) after overnight fixation in 4% neutral formaldehyde. After dewaxing, hydration, and thermal antigen recovery, 3% hydrogen peroxide and 3% bovine serum albumin were used to suppress endogenous peroxidase. Subsequently, primary and secondary antibodies were incubated in sequence. Finally, 3.3′-diaminobenzidine (P0202, Beyotime) was used for color visualization and hematoxylin re-staining. The expression levels of hub genes were detected using ImageJ 1.48v and visualized using GraphPad Prism Version 9.0.0. The following primary antibodies were used: anti- TLR3 (Ab62566, Abcam, Cambridge, UK), anti-MLKL (phospho S345) (Ab196436, Abcam), anti-HMGB1 (GB11103, Servicebio, Wuhan, China), and anti-NDRG2 (GB111833, Servicebio).

Statistical Analysis

Statistical analyses were performed using SPSS version 23.0 (IBM). Quantitative variables with normal distribution were expressed as mean ± SD, non-normally distributed as median (25, 75th percentiles. Significant differences were verified using an independent-sample t-test or Wilcoxon rank-sum test. Qualitative variables were expressed as counts and rates, and significant differences were verified using the chi-square test or Fisher’s exact test. The cutoff for statistical significance was set at p < 0.05.

Results

Landscape of Necroptosis Genes

In total, 1218 genes were identified, of which 860 were upregulated and 358 were downregulated (Figure 2A). We generated a heat map showing the expression of 17 necroptosis genes in the normal and CHD groups (Figure 2B). Box and scatter plots were used to visualize the similarities and differences in the transcription levels of all necroptosis genes in the CHD and control groups (Figure 2C). The findings revealed that MLKL (p < 0.01), TLR3 (p < 0.05), and HMGB1 (p < 0.05) were upregulated and NDRG2 (p <0.05) was downregulated in the CHD group compared to the control group. Furthermore, we analyzed the chromosomal localization of necroptosis genes to generate a panorama of genes on chromosomes (Figure 2D). These results indicated that certain genes had extremely close chromosomal positions, implying genomic linkage and potential similarity in their expression patterns at the transcriptomic level.

Figure 2 Landscape of necroptosis genes. (A) Volcano map of necroptosis genes. Red dots in the graph represent upregulated genes and blue dots represent downregulated genes. The gray dots represent genes without statistical significance. (B) Heat map of necroptosis genes. Rows stand for genes, columns stand for samples, and samples as well as different types of genes are marked using different-colored blocks. Red in the graph represents high expression values while blue represents low expression values. (C) Combined box and scatter plots of necroptosis genes. Genes are plotted on the X-axis, while gene expression values are plotted on the Y-axis, sample groupings are distinguished by different colors and faceted according to gene type. The middle line of the box plot is the median value, the upper box line is the upper quartile, and the lower box line is the lower quartile. *p < 0.05; **p < 0.01. (D) Chromosome positioning loops. The outer circles are chromosomes and the inner circles are labeled with the positions of these genes on the chromosomes. Abbreviations: ns, not statistically significant.

Correlation Between Necroptosis Genes

We assessed the association between these genes in both groups to further analyze the relationship between necroptosis gene expression in the CHD and control groups. The correlation heat map is shown in Figure 3A. The results displayed that the necroptosis genes had varying correlations in the CHD and control groups. For instance, the negative correlation between NDRG2 and TLR2 was stronger in the CHD group than in the control group (−0.29945 vs −0.213958), suggesting that the expression of necroptosis genes undergoes changes during CHD. Further, we selected gene pairs with high correlation coefficients (p < 0.001) in the CHD group and plotted the correlation scatter plots and fitted correlation curves (Figure 3B–D). The findings showed positive and statistically significant correlations between MLKL and TLR4 (r=0.527, p < 0.001) and between MLKL and RIPK3 (r=0.492, p < 0.001). In contrast, NDRG2 and TLR2 levels displayed a statistically significant negative correlation (r=−0.299, p < 0.001).

Figure 3 Correlation between necroptosis genes. (A) Correlation analysis of necroptosis genes between the CHD and control groups. Red represents positive correlation; blue represents negative correlation. * p < 0.05; ** p < 0.01. (B–D) Gene pairs with high correlation coefficients. Each point in the graph represents a patient sample, the straight line is the correlation fit curve, the shaded part is the confidence interval, and the outer part of the graph is the density curve. (B) MLKL-TLR4 correlation. (C) RIPK3-MLKL correlation. (D) TLR2-NDRG2 correlation. Abbreviations: CHD, coronary heart disease.

Necroptosis Gene Diagnostic Model

As necroptosis genes play a crucial biological role in the progression of CHD, we built a diagnostic model of CHD utilizing all necroptosis genes. First, 17 necroptosis genes were screened using the LASSO regression method. Subsequently, the optimal lambda values were set resulting in the identification of 12 genes (RIPK1, RIPK3, MLKL, TLR3, TNFRSF1A, NR2C2, HMGB1, CXCL1, USP22, ALDH2, EZH2, and NDRG2) that were preserved (Figure 4A and B). Furthermore, to narrow down the candidate genes, a random forest method was used to identify essential trait variables, and the role of the genes in disease were ranked. The results revealed that the top 10 genes, ranked by importance to the disease, were NDRG2, TLR3, TNFRSF1, ALDH2, RIPK1, MLKL, HMGB1, TLR2, NR2C2, and USP22 (Figure 4C and D). Finally, necroptosis gene expression profiles and the two approaches mentioned above were combined to screen for genes, and the intersection of the three genes was considered a candidate gene for the construction of the diagnostic model. The results showed that TLR3, MLKL, HMGB1, and NDRG2 are important necroptosis genes involved in CHD progression and could be appropriate candidates for the diagnosis model (Figure 4E).

Figure 4 Construction of the diagnostic model. (A) LASSO regression curve. This figure shows the convergence screening process of LASSO regression for 17 genetic features, the log lambda value is plotted on the X-axis, while the regression coefficient is on the Y-axis, and the different color lines represent different features. (B) Lambda value selection curve. This graph is used to select the best lambda value for the regression model, and usually the lowest point, that is, the dashed line in the graph, is selected as the best lambda value. (C) Random forest model construction, the X-axis shows the number of selected forests, and the error rate is plotted along the Y-axis. (D) Random forest screening of significant genes. The Y-axis shows the gene names, and the gene importance values are plotted along the X-axis. (E) The Venn diagram integrates three gene sets for Lasso, RF Top, and DE-NAGs, with four genes overlapped. The four genes were included in the construction of the diagnostic model. Abbreviations: LASSO, least absolute shrinkage and selection operator; RF Top 10, top 10 genes ranked by random forest screening for gene importance; DE-NAGs, differentially expressed necrotic apoptosis genes.

Subsequently, to visualize the contribution of the four hub genes to the disease, we illustrated the results using a nomogram (Figure 5A). The results revealed that each of the four genes had a significant effect on the prediction model with the following logistic regression equation:

Figure 5 Validation of a diagnostic model for necroptosis genes. (A) Nomogram plot with the predictors on the left and the scale on the right. * p < 0.05. (B) Calibration curve, the results show a good fit between the predicted and actual values of the model. Nonadherence represents proportion of people treated as prescribed or recommended. (C) DCA. The risk threshold is plotted on the X-axis, and the net benefit rate is plotted on the Y-axis. None (black line) indicates 0 net benefit rate, all (gray line) indicates that all samples received the intervention, and nomogram (blue line) indicates the model curve. Abbreviations: DCA, decision curve analysis.

The efficiency of the model predictions was further assessed using calibration and DCA. The calibration curve findings demonstrate that the prediction model provided a close fit to the ideal line (Figure 5B). The DCA revealed that the prediction model constructed with the hub gene for necroptosis (blue curve) was associated with a superior level of benefit to patients compared to the model constructed with all genes (gray curve). This suggests that necroptosis hub gene-based predictive model has diagnostic efficacy with respect to CHD and may benefit patients within a certain range.

Interaction Networks of Necroptosis Genes

As genes are commonly associated with each other, and this connection is strengthened for genes regulating the same biological process, a PPI network was built using the STRING database to further analyze the interactions between necroptosis genes. Additionally, 17 necroptosis genes were grouped into two clusters using STRING’s built-in cluster analysis, and Cytoscape was used for visualization (Figure 6A). Genes with high node degrees are pivotal to networks. Therefore, we visualized and highlighted the top five genes in the network based on the node degree. These were RIK1, RIK3, MKLK, TLR3, and TLR4 (Figure 6B). To fully address the genetic origins and regulatory networks of these hub genes, the miRNet database was used to project their miRNAs and transcription factors (Figure 6C). These findings indicated that these genes contained both unique and shared miRNAs or transcription factors. Consequently, they may be exposed to the same regulatory mechanisms, ultimately exhibiting comparable biological roles.

Figure 6 PPI network. (A) PPI network of 17 necroptosis genes. The colors represent the expression of genes, the darker the color the higher the expression. Red and blue indicated high- and low-expressed genes, respectively. The different shapes of the genes represent different clusters. Four of the genes were not linked to other genes, so they are not shown in the network. (B) Five hub genes. The sub-networks of genes in the PPI network were extracted with the cytoHubba tool in Cytoscape, and the redder the color of the node, the greater the degree of the node in the original network. The top five genes were RIK1, RIK3, MLKL, TLR3, and TLR4. (C) The TFs-miRNA prediction network of hub genes. Green inverted triangles are transcription factors, cyan circles are miRNAs, and other colors are necroptosis genes. Abbreviations: PPI, protein-protein interaction network; TFs, transcription factors.

ssGSEA Scores Based on Necroptosis Genes

Necroptosis genes are crucial to CHD progression; therefore, we used the expression of necroptosis genes to score ssGSEA in samples from patients in the CHD group. The median ssGSEA score was also set as the cutoff value to divide the samples into high- and low-score clusters. Uniform manifold approximation and projection (UMAP) analysis revealed that the samples in the CHD group exhibited distinct subgroups characterized by high and low scores (Figure 7A). Differential gene expression analysis was performed between these two groups, totaling 2080 genes after screening for statistical significance thresholds, with 587 upregulated and 1493 downregulated genes. Volcano plots visualizing the differentially expressed genes are shown in Figure 7B. The top 40 genes in terms of differentially changed ploidy are shown by their expression patterns on a heat map (Figure 7C). The risk factor plot depicted the expression of 17 necroptosis genes after dividing the disease sample into two clusters using the median ssGSEA score as the threshold value (Figure 7D). The combined box and violin plot graphs visualized the parallels and variations between the two types of samples with respect to the expression levels of all necroptosis genes (Figure 7E). As shown in Figure 7, 14 of the 17 necroptosis genes, which did not include PGAM5, HMGB1, EZ2H2, and NDRG2, demonstrated significant differences (p < 0.05) between the two groups, indicating that necroptosis genes may be the crucial characteristics to identify CHD and warrant further investigation.

Figure 7 Necroptosis scores of the samples. (A) CHD group samples could be distinguished into high- and low-expression groups of necroptosis scores. (B) Volcano plot of differentially expressed genes. The log2 (Fold Change) is plotted on the X-axis and -log10 (adj value) is plotted on the Y-axis. Each point represents a gene: blue is for downregulated genes, red shows upregulated genes, and gray represents genes with no significant change in expression. (C) Heat map. The genes with the top 40-fold change in behavioral differences are listed as samples. The groupings and gene types to which the samples belong are marked using different color blocks. Blue in the graph represents low expression values and red represents high expression values. (D) Heat map of necroptosis scores. The CHD group was divided into high- and low-score groups using the median necroptosis score as the cutoff value. The heat map below represents the expression of necroptosis genes between the two groups. (E) Combined box and violin plot. Genes are plotted on the X-axis, while gene expression values are plotted on the Y-axis. The sample groupings are distinguished by different colors, the middle line of the box plot is the median value, the upper box line is the upper quartile, and the lower box line is the lower quartile. The statistical test used was the Wilcoxon rank sum test. Abbreviations: CHD, coronary heart disease; ns, not statistically significant. Notes: **p < 0.01; ***p < 0.001.

Biological Variations Between Groups

To delve deeper into potential molecular mechanisms or signaling pathways associated with differentially expressed genes, we conducted GO and Kyoto Encyclopedia of Genes and Genomes enrichment analyses specifically for genes with absolute logFC values > 0.5 (Figure 8A and B, Table S4). The GO enrichment analysis results revealed 185 enhanced GO terms in total, most of which were related to myocardial contraction and blood circulation, namely the muscle system process, regulation of ion transmembrane transport, and muscle contraction. KEGG enrichment analysis revealed that the DEGs were mostly linked to calcium signaling and neuroactive ligand receptor interactions.

Figure 8 Biological differences between groups. (A) GO enrichment results. The gene ratio, that is, the total number of genes/differentially expressed genes enriched to the term, is plotted on the X-axis, whereas the names of the GO terms are plotted on the Y-axis. The size of the dot represents the number of genes enriched to the term, and the color indicates the corrected p-value; the smaller the p, the closer to red it is. (B) Circle plot of GO and KEGG enrichment results. The red dots represent upregulated genes, while blue represents downregulated genes. The bar graph of the inner ring represents the z-score of the enriched entries. (C) Results of GSEA analysis. Each color represents one pathway. Only the top two most significantly upregulated pathways and the top three most significantly downregulated pathways are shown here. (D) This figure shows the difference in expression of molecules related to major histocompatibility complex signaling pathways between high and low groups of necroptosis score in a heat map. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; GSEA, Gene Set Enrichment Analysis.

Subsequently, we performed GSEA on each gene to determine their biological roles. These findings (Table S5) revealed that interleukin (IL)-6 and major histocompatibility complex (MHC)-related immune and inflammatory signaling pathways were activated in the CHD group. In contrast, voltage-gated potassium and calcium channels were inhibited in this group (Figure 8C). Furthermore, we visualized the expression patterns of MHC signaling pathway molecules in the two groups. These findings demonstrate that MHC-related signaling molecule levels were markedly increased in the high-score group (Figure 8D). This result suggests that necroptosis genes might influence the progression of CHD through their involvement in immune-inflammatory processes.

Immune Infiltration Analysis

Enrichment analysis revealed the differences in immune-inflammatory responses between the two groups. To further examine the level of immune infiltration in the two groups, we calculated the scores of 28 immune cell types for all samples using the ssGSEA method, visualized them using violin plots, and tested them for statistical significance (Figure 9A). The findings showed that the scores of three of the 28 immune cell types, namely activated B cells, activated dendritic cells, and neutrophils, were substantially different in the CHD group (Figure 9A). Additionally, the scores of the 28 immune cells in the CHD group were obtained using ssGSEA, and the disparities between the necroptosis gene scores of the 28 immune cells in the high- and low-score groups were measured using violin plots. The findings suggested that, with the exception of activated B, activated CD8 T, CD56 bright natural killer, CD56 dim natural killer, effector memory CD4 T, immature B, type 17 T helper, type 1 T helper, and T follicular helper cells, the remaining 19 immune cells exhibited significant differences between the two groups. Specifically, the high score group had significantly higher levels of neutrophils (p < 0.001), effector memory CD8 T cell (p < 0.001), and monocyte (p < 0.001) infiltration (Figure 9B). Furthermore, to assess the association between immune cells, the correlation coefficient was computed independently for both the CHD and control groups (Figure 9C and D). These findings suggest a positive correlation between immature dendritic cells, macrophages, mast cells, myeloid-derived suppressor cells, memory B cells, monocytes, and neutrophils. In particular, monocytes and neutrophils showed an interesting association in the CHD group (cor = 0.4047, p = 1.005E-08), similar to that of monocytes and myeloid-derived suppressor cells (cor = 0.537, p = 2.665E-15).

Figure 9 ssGSEA immune infiltration assessment. (A) Violin plot of the distribution of immune scores in normal and disease groups. The X-axis shows 28 immune cell types, while the level of immune infiltration is plotted on the Y-axis. Each color represents a sample subgroup, and the statistical test used was the Wilcoxon rank sum test. (B) Violin plot of immune score between groups with high and low necroptosis scores. The X-axis shows 28 immune cell types, whereas the level of immune infiltration is plotted on the Y-axis. Each color represents one sample grouping. The statistical test used was the Wilcoxon rank sum test. (C) Heat map of the correlation between 28 immune cells in control samples. (D) Heat map of the correlation of 28 immune cells in the CHD group. Abbreviations: ssGSEA, single-sample gene set expression analysis; CHD, coronary heart disease.

Similarly, for all samples, 22 immune cell scores were established using the CIBERSORT method, and the immune score results for all samples are shown in Figure 10A. To investigate the association between individual immune cells within the CHD group, we conducted a correlation study on immune cells. The findings revealed that neutrophils were correlated with T cells, whereas resting natural killer (NK) cells and neutrophils had a strong negative correlation (Figure 10B). In addition, a combination of box and scatter plots was used to visualize and statistically evaluate the variations in immune cells between the CHD and control group. The findings displayed that naïve B cells (p < 0.05) and monocytes (p < 0.05) were substantially more prevalent in the CHD group. The distribution of neutrophils was higher in the CHD group than in the control group; however, the difference was not statistically significant (p = 0.13) (Figure 10C). The scores of 22 immune cells were estimated using the CIBERSORT method, and the differences between the two groups were visualized using a mix of box and scatter plots. The findings displayed that eight of the twenty-two immune cells were statistically different between the high- and low-score groups; these included neutrophils (p < 0.001), monocytes (p < 0.05), and resting mast cells (p < 0.01), which were considerably increased in the high-score groups; whereas plasma cells (p < 0.01), CD8 T cells (p < 0.01), M0 macrophages (p < 0.05), M2 macrophages (p < 0.05), and resting mast cells (p < 0.01) were significantly downregulated in the high-score group (Figure 10D). This suggests that necroptosis genes may be linked to the activation of monocytes and neutrophils, depletion of CD8+ T lymphocytes, and suppression of plasma cells and M2 macrophages.

Figure 10 CIBERSORT immune infiltration assessment. (A) Plot of the percentage of immune cells in each sample. The vertical coordinate is the percentage of immune cells in each sample, and different colors represent different samples. (B) Correlation analysis between immune cells. (C) Combined box plot and scatter plot graphical analysis of immune scores in CHD and control groups. The X-axis is for the 22 immune cells, while the level of immune infiltration is plotted on the Y-axis. Each color represents a sample subgroup, and the Wilcoxon rank sum test was used. (D) Combined box and scatter plot graph analysis of immunity scores between high- and low-score groups for necroptosis scores. The X-axis is for the 22 immune cells, whereas the Y-axis shows the level of immune infiltration. Each color represents a subgroup, and the Wilcoxon rank sum test was used. p < 0.05 was set as the cutoff for significant differences in the level of immune infiltration between the two groups.

To analyze the association between hub genes and immune cells, we measured the links between hub genes and the degrees of all forms of immune cell infiltration and visualized the results of hub gene-immune cell pair correlations with a lollipop plot (Figure 11A–D). We found a positive correlation between MLKL-neutrophils, MLKL-monocytes, MLKL-resting mast cells, NDRG2-CD8 T cells, NDRG2-naïve CD4 T cells, TLR3-CD8 T cells, TLR3-resting NK cells resting, and TLR3-activated CD4 memory T cells. Further, there was a negative correlation between MLKL-CD8 T cells, MLKL-plasma cells, TLR3-plasma cells, TLR3-macrophages M0, H-plasma cells, HMGB1-macrophages M2, NDRG2-monocytes, NDRG2-neutrophils, and NDRG2-plasma cells. The exact magnitudes of the correlations are shown in Figure 11. This suggests that these crucial necroptosis genes may directly influence the infiltration of specific immune cells, which warrants further research.

Figure 11 Correlation analysis between hub genes and immune cells. (A) MLKL-22 immune cell relationship pairs; (B) TLR3-22 immune cell relationship pairs; (C) HMGB1-22 immune cell relationship pairs; (D) NDRG2-22 immune cell relationship pairs.

In vivo Hub Gene Validation

A total of 102 patients were enrolled in this study (56 and 46 patients in the CHD and control groups, respectively). The common clinical features of the patients are presented in Table S6.

Compared with the control group, PCR revealed that the relative mRNA expression of TLR3 (p < 0.05), MLKL (p < 0.05), and HMGB1 (p < 0.05) was elevated in the CHD group, whereas the relative expression of NDRG2 (p < 0.01) was decreased (Figure 12).

Figure 12 The qPCR results of hub genes between CHD and controls. Abbreviations: qPCR, quantitative polymerase chain reaction. Notes: *p < 0.05; **p < 0.01.

Hematoxylin and eosin staining showed that APOE^−/− mice had a larger area of aortic root plaque formation than C57BL/6 mice (Figure 13A). The expression of TLR3, p-MLKL, HMGB1, and NDRG2 in the aortic root sections of mice were detected using immunohistochemistry (Figure 13B). Compared to those in the control group, the positive regions of TLR3, p-MLKL, and HMGB1 were higher in the AS group; however, the positive area of NDRG2 was lower (Figure 13C).

Figure 13 The immunohistochemical results of hub genes between CHD and controls.(A) HE stains showed the control group in C57BL/6J mice and atherosclerosis group in APOE−/− mice. (B) Immunohistochemical staining of p-MLKL, TLR3, HMGB1, and NDRG2 in the aortic root. (C) Analysis of p-MLKL-positive area, TLR3-positive area, HMGB1-positive area, and NDRG2-positive area in immunohistochemistry stains, n=4. Abbreviations: CHD, coronary heart disease; HE, Hematoxylin and eosin. Notes: *p < 0.05; **p < 0.01.

Discussion

In clinical practice, invasive angiography and coronary computed tomography are currently the main methods for diagnosing CHD. Due to the high cost, radiation risk, and especially the invasive nature of the former, the application of these methods is constrained in the screening and diagnosis of CHD, especially in cases of subclinical AS. Therefore, there is an urgent need to identify effective peripheral blood biomarkers. AS is a lipid-driven, maladaptive vascular wall disease with persistent inflammation.⁴⁶ Necroptosis, a unique proinflammatory type of cell death,⁴⁷ has been shown in recent studies to have a crucial role in the evolution of atherogenesis.¹⁵ UMAP analysis revealed that the CHD group could be well categorized into high- and low-score groups in our study (Figure 7A), suggesting that necroptosis is an important differentiator of CHD. Hence, this study aimed to identify whether necroptosis genes could serve as biological markers of CHD.

In this study, four diagnostic candidate genes, TLR3, MLKL, HMGB1, and NDRG2 were identified. TLR3 is a key pattern-recognition receptor in the immune system that detects ligands for the sustained release of endogenous RNA from inflammatory tissue and necrotic cells in atherosclerotic lesions, thereby promoting and exacerbating the inflammatory responses that ultimately lead to the further exacerbation of AS^48,49 Cen et al highlighted the atherogenic and pro-inflammatory effects of TLR3, as evidenced by the significant upregulation of TLR3 expression in PBMCs from patients with AS.⁵⁰ Additionally, the activation of TLR3 in mice has been demonstrated to promote foam cell formation.⁵⁰ This finding is consistent with our study, which shows significant elevation of TLR3 expression in the peripheral blood of patients in the CHD cohort, as well as in AS plaques in mice. The TLR3-specific inhibitor SMU-CX24 was tested in their experiments for its anti-inflammatory and anti-atherosclerotic actions in animal models.⁵⁰ This provides a promising strategy for the development of therapies for AS.

MLKL is a major executor of necroptosis; its oligomerization has recently been shown to determine the kinetics and threshold of necrotic cell death.⁸ Therefore, the importance of MLKL in necroptosis cannot be overstated. Studies have shown that MLKL actively contributes to the development of the inflammatory process and the establishment of the necrotic core of atherosclerotic plaques.^51,52 Consistent with our findings, Karunakaran et al showed elevated MLKL expression in patients with AS, and Rasheed et al showed that the inhibition of MLKL resulted in a decrease in necroptosis and necrotic cores in plaques. However, there was no overall reduction in plaque area. Further investigation revealed a novel role of MLKL in protecting splenic endothelial cells from the excessive recruitment of circulating leukocytes, thereby exerting an anti-atherogenic effect.⁵³ However, the dual role of MLKL in AS requires further investigation.

At 50 years since the discovery of HMGB1, a chromatin structural protein, various novel functions continue to be uncovered. In AS, the metabolic and biochemical environment generated by lipid entry into the vessel wall drives the production of receptors for advanced glycosylation end products, which contain HMGB1.^54,55 Furthermore, HMGB1 may be released as a result of various modes of cell death, including apoptosis, pyroptosis, necrosis, and ferroptosis.⁵⁶ HMGB1 contributes to the development of inflammation via a damage-associated molecular pattern (DAMP).⁵⁶ Recent studies have shown that sympathetic activation and HMGB1, an alarmin released from the heart or brain by necrosis after myocardial infarction or stroke, might play a key role.⁵⁷ In the present study, HMGB1 was highly expressed in the peripheral blood of patients with CHD and in AS plaques of APOE^−/− mice. Therefore, HMGB1 is actively involved in AS and deserves further investigation.

NDRG2, the downstream regulatory gene of N-myc, has received much attention in oncology as a tumor suppressor or stress response gene.^58–60 However, little has been reported in the cardiovascular field. In our study, the expression of NDRG2 was decreased in the peripheral blood of patients with CHD, as well as in the atherosclerotic plaques found in mice. The results of the study by Liu et al⁶¹ were inconsistent with our conclusions, possibly for the following reasons: first, the samples examined in the two studies were different; Liu et al examined RAW264.7 and human aortic plaque tissue, whereas our study examined peripheral blood from patients with CHD and aortic plaque tissue from mice. Additionally, the sample sizes of both studies were small and their representativeness may be limited. Further in-depth research is required to investigate the expression of NDRG2 in CHD and its effects.

We conducted GO/KEGG enrichment analyses and GSEA to gain insight into the biological mechanisms and signaling pathways associated with the necroptosis genes. Our findings indicate that they potentially contribute to myocardial contraction and circulation, neuroendocrine regulation, ion channels, and signal transduction. In particular, IL-6- and MHC-related immune and inflammatory signaling pathways were activated in the CHD group. Further analysis of the MHC signaling pathway revealed a significant upregulation of MHC-related signaling molecules in the high-score group, indicating that necroptosis genes might impact the progression of CHD by immune inflammation.

Necroptosis is regarded as a regulated mode of inflammatory cell death. TLR3, MLKL, and HMGB1 play various roles in necroptosis and play proinflammatory roles in CHD. TLR3, a member of the Toll-like receptor family, recognizes pathogenesis-associated molecular patterns and DAMPs within vascular plaques,⁵⁰ triggering necroptotic cell death within the plaque. MLKL is an executioner of necroptosis and aids in the formation of the necrotic core of the plaque^51,52 by forming pores in the plasma membrane or inducing the opening of ion channels in the plasma membrane,^11,56 which impairs the integrity of cells, triggering cell death. HMGB1, a key immune mediator released in various models of cell death and injury in atherosclerotic plaques, plays a key role in mediating cell death-induced inflammation as a DAMP. However, the role of NDRG2 in CHD remains to be elucidated, and human genetic evidence suggests that IL-6 is not only directly involved in the atherosclerotic process but also has a possible causal relationship with AS.^62,63 A recent study by Ridker et al provides data from a Phase 2 clinical trial (rescue trial) in patients at high risk for AS, in which the administration of zitivizumab, a fully human monoclonal antibody directed against IL-6, decreased a variety of systemic inflammation biomarkers related to atherosclerotic thrombosis.⁶⁴ Therefore, IL-6 receptor blockage appears to provide a novel strategy for the management of CHD. To the best of our knowledge, there is no research on the function of any hub genes in CHD via the IL-6 signaling pathway, which deserves further in-depth attention.

Immune cells are activated as a result of continuous inflammatory stimuli or failure of inflammatory abatement, resulting in chronic inflammation, which is a defining feature of AS. In this study, we determined the landscape of immune cell infiltration in the peripheral blood of patients with CHD using ssGSEA and CIBERSORT. The results were as follows: B cells, neutrophils, monocytes, and activated dendritic cells were more abundantly found in the peripheral blood of patients with CHD. Neutrophils and monocytes exhibited greater infiltration in the high-score group, whereas CD8+ T cells, M0/M2 macrophages, plasma, and resting mast cells infiltrated less.

In this study, the levels of neutrophils and monocytes were increased in both the CHD and high-score groups, and these two cells were highly correlated in the CHD group. Neutrophils interact with monocytes; for example, they attract monocytes by secreting chemotactic molecules and activate macrophages by extruding nuclear material to create neutrophil extracellular traps (NETs).⁶⁵ NETs induce plaque erosion and platelet accumulation, ultimately resulting in thrombosis.⁶⁶ In general, neutrophils exert pro-atherogenic actions;⁶⁷ in both mice⁶⁸ and humans,⁶⁹ AS severity is correlated with an increase in the blood monocyte pool. The inhibition of monocyte recruitment and targeting of monocyte proliferation are considered promising therapeutic targets^70,71 Recent studies have shown that various risk factors for AS, such as hypercholesterolemia, aging, stress, and inflammation, can promote hematopoietic stem cell proliferation, leading to systemic mononucleosis,⁷² which promotes the progression of AS.⁷³

Immune cells are typically adaptable and exhibit different functional subtypes and phenotypes depending on the stage of AS development, tissue location, and microenvironment.⁷⁴ Accordingly, they can exhibit either pro- or anti-atherogenic effects. First, one example of plasticity is the polarization of macrophages into the respective pro- and anti-inflammatory M1 and M2 types. Studies have demonstrated that M2 macrophages contribute to plaque stability by secreting collagen and improving the elimination of apoptotic cells.⁷⁴ The number of M2 macrophages within the plaque generally declines as the lesion progresses, whereas the number of M1 macrophages gradually rises. The relative proportions of macrophage subpopulations, rather than absolute numbers, serve as a better indicator of AS lesion progression.⁷⁵ Recent single-cell studies on human⁷⁶ and mouse⁷⁷ AS plaques have identified additional macrophage populations with different inflammatory properties.⁷⁸ The description of macrophage polarization as M1-M2 has become simplistic and ill-defined, necessitating further refinement of our knowledge through research.⁷⁹ Furthermore, distinguishing between pro-AS and anti-AS B cell subsets is crucial. However, one major limitation is that most subsets have similar surface markers, making it challenging to isolate the unique impacts of different B cell subtypes.⁸⁰ In addition, results on CD8+ T cells, as assessed by ssGSEA and CIBERSORT between the high- and low-score groups, seem to give contradictory results; however, this is not the case, as various subtypes of CD8+ T cells exist, exerting synergistic or completely opposite effects. Effector CD8+ T cells release cytotoxic particles, especially perforin and granzymes, which induce lysis of the target cells,⁸¹ thereby promoting macrophage death and necrotic core formation. Therefore, it is reasonable to speculate that there was higher infiltration of effective CD8+ T cells in the high-score group as determined by ssGSEA. The number of CD8+ T cells was dramatically reduced in the high-score group, as determined by CIBERSORT. This could be due to T-cell depletion in advanced AS lesions. Progressive adaptation and possible loss of T-cell function occur in chronic, sustained inflammatory responses.⁸² A previous study demonstrated a decrease in CD8+ T-cell levels as atherosclerotic lesions increased in MHCI-deficient mice.⁸³

In this study, we found the potential diagnostic value of necroptosis genes in CHD based on the expression analysis in public databases and our proprietary data, but there are certain limitations. First, the clinical sample of our study was limited; to evaluate the diagnostic efficacy of candidate genes, further expansion of the sample or large-scale community risk screening for CHD is necessary. Second, we built a diagnostic model and externally validated it using public database datasets, but we were unable to validate its diagnostic performance in our private patient cohort by plotting calibration curves and calculating the C-statistic. We intend to test and validate our diagnostic model in a broader patient cohort and explore its clinical diagnostic applications. For now, we advise caution regarding our model’s present evaluation and outcomes as a potential foundation for future study.

Conclusion

In our study, four hub genes with a common biological function, necroptosis, were used to construct a diagnostic model, which might contribute to improved clinical management with favorable diagnostic efficacy for CHD. Therefore, TLR3, MLKL, HMGB1, and NDRG2 may be potent peripheral blood biomarkers of CHD. Necroptosis-related genes may play a role in CHD via immune inflammatory responses and immune infiltration.