Introduction

Atherosclerosis (AS) is a chronic inflammatory disease characterized by the deposition of fats within the walls of arteries. During the progression of the disease, plaques form within the arteries, restricting blood flow and leading to serious complications such as heart attacks and strokes.¹ AS is a life-threatening manifestation of cardiovascular disease, primarily affecting large and medium-sized arteries, including the coronary arteries, carotid arteries, aorta, and cerebral arteries. Epidemiological studies indicate that in 2019, over 17 million people worldwide died from cardiovascular diseases, with 85% of these deaths attributed to stroke and other conditions caused by AS.² As of 2020, nearly 2 billion people globally were estimated to have carotid atherosclerosis.³ In summary, AS-related cardiovascular diseases impose a significant global health and economic burden.⁴ It is widely recognized that the development of AS is associated with dyslipidemia and hypercholesterolemia, conditions that can be triggered by factors such as age, genetics, diet, lifestyle, and underlying health conditions like hypertension, diabetes, and obesity.⁵ Currently, the primary pharmacological treatments for AS are statins or PCSK9 inhibitors, aimed at controlling low-density lipoprotein (LDL) cholesterol levels in the blood.⁶ Although these therapies are effective in reducing LDL levels, the incidence of major adverse cardiovascular events (MACEs) remains unmitigated. Therefore, gaining a deeper understanding of the pathophysiological mechanisms of AS and developing more effective diagnostic and therapeutic approaches are critical areas of ongoing research.

Cellular senescence is a biological process associated with the deterioration of cellular structure and function, and it is a major cause of age-related diseases.⁷ Currently, cellular senescence is widely recognized as an important target for the prevention of cardiovascular diseases. AS is an age-related disease, and its progression may involve various senescent cardiovascular cells, including cardiomyocytes, endothelial cells, vascular smooth muscle cells, fibroblasts/myofibroblasts, and T cells.⁸ Additionally, multiple forms of cell death are involved in the progression of AS, including autophagic cell death, apoptosis, pyroptosis, necroptosis, and ferroptosis.⁹ Ferroptosis is an iron-dependent form of programmed cell death. The main mechanism of ferroptosis is the catalysis of highly expressed polyunsaturated fatty acids on the cell membrane by ferrous iron or lipoxygenase, leading to lipid peroxidation and cell death.¹⁰ Numerous studies have found that ferroptosis can participate in the onset and progression of AS through different pathways.¹¹ In advanced AS, cell ferroptosis induced by functional defects can accelerate endothelial dysfunction and the formation of vulnerable plaques. Additionally, dysregulated intracellular iron ions can damage macrophages, vascular smooth muscle cells, and endothelial cells, affecting many AS risk factors or pathological processes, such as oxidative stress, inflammation, and lipid metabolism disorders.¹² Recent studies have shown that ferroptosis and cellular senescence can act synergistically to regulate the development of atherosclerosis (AS). For example, activation of the ferroptosis signaling pathway in vascular smooth muscle cells (VSMCs) induces the loss of nicotinamide adenine dinucleotide (NAD⁺) and promotes cellular senescence. It also enhances the release of senescence-associated secretory phenotypes (SASPs).¹³ Conversely, senescent cells may disrupt iron homeostasis, leading to iron accumulation. This can increase cellular resistance to ferroptosis and alter ferroptotic sensitivity.¹⁴ These findings suggest a potential positive feedback loop between ferroptosis and senescence during AS progression. Understanding this interplay is important for developing new ferroptosis-based therapies for aging-related cardiovascular diseases. However, whether these two processes intersect functionally in AS and which factors regulate this crosstalk remain unclear. To date, no systematic studies have addressed this question. Meanwhile, peripheral biomarkers are becoming increasingly important as non-invasive tools for monitoring coronary artery disease (CAD). They show great promise in early detection, prognosis evaluation, and subtype classification.¹⁵ Some studies also suggest that senescent and dying cells can recruit myeloid cells, such as immunosuppressive macrophage subsets, to mediate immune suppression directly or indirectly.¹⁶ Given the close link between AS and immune therapy response, optimizing immunotherapy strategies for AS has been a research focus.¹⁷

In recent years, bioinformatics has provided new approaches for the treatment of atherosclerosis. Machine learning has been widely applied across various research fields as a powerful tool for identifying important diagnostic markers.¹⁸ Single-cell technology is an ideal tool for the comprehensive and unbiased analysis of cellular heterogeneity. High-throughput, high-resolution single-cell sequencing technology can be used to analyze the complex cellular and molecular composition of atherosclerotic plaques. It can also be used to examine the state, origin, dynamic transformation trajectories, intercellular communication, and molecular regulatory mechanisms of cells within plaques.¹⁹ The combination of these technologies facilitates a comprehensive and in-depth understanding of the mechanisms of AS onset and progression. This approach provides new directions for the diagnosis and treatment of AS and the study of its immunological mechanisms.

Methods

Data Acquisition and Processing

The genome-wide expression data used in this study were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Specifically, the GSE10097 and GSE132651 datasets, both based on the GPL96 platform, were downloaded using the GEOquery package in R. The original expression matrices and sample annotation files were retrieved with the getGEO function. Quantile normalization was applied to standardize the data. The combined datasets included 82 atherosclerosis cases and 41 controls (In this study, AS refers exclusively to atherosclerosis and is not related to the term “asymptomatic.”). The GSE10097 dataset was used as the training set. It was generated from transcriptome sequencing and contains 104 arterial tissue samples, including 69 Atherosclerosis samples and 35 non-Atherosclerosis controls. The GSE132651 dataset served as the validation set. It includes 19 subjects younger than 50 years old (mean age 37 years, range 22–49), who were divided into normal endothelial function (n = 6) and impaired endothelial function (n = 13) groups based on coronary endothelial reactivity testing. Single-cell RNA sequencing data were obtained from the GSE159677 dataset on the GPL18573 platform, which includes carotid artery tissue samples from six Atherosclerosis patients. This study was conducted in accordance with national ethical guidelines. According to Article 32, Items 1 and 2 of the “Administrative Measures for Ethical Review of Life Science and Medical Research Involving Humans” (issued February 18, 2023, China), ethical approval is exempted because this study used only publicly available datasets and did not involve human intervention or the collection of personally identifiable information.

Ferroptosis-related genes were obtained from the FerrDb V2 database (http://www.zhounan.org/ferrdb/current/), including a total of 905 genes. Senescence-related genes were obtained from the GenAge database (https://genomics.senescence.info/genes/), with a total of 866 genes.

Screening of Differentially Expressed Genes

Differentially expressed genes (DEGs) in the GSE10097 dataset were identified using the “Limma” package in R software. Background correction and data normalization were performed using the NormExp and Normalize Between Arrays functions of the “Limma” package. The screening parameters for DEGs were set with an absolute value of log FC ≥ 0.5 and a P-value < 0.05. After identifying the DEGs, heatmaps and volcano plots of the DEGs were created using the “pheatmap” and “ggplot2” packages in R software.

Weighted Gene Co-Expression Network Analysis

Weighted gene co-expression network analysis (WGCNA) is used to identify gene modules with high correlation. In this study, WGCNA was performed using the “WGCNA” package in R software.²⁰ Hierarchical clustering was performed on the study samples to detect and remove outliers. A soft-thresholding power of β = 6 was selected to construct a scale-free network and calculate the topological overlap matrix (TOM). The correlation between each module and the disease group was then assessed. Modules with a high correlation coefficient were considered key modules closely related to the pathogenesis of atherosclerosis (AS). The threshold for module–phenotype correlation was set as a correlation coefficient > 0.5, with a significance level of P < 0.05. Genes from the key modules were intersected with senescence-related genes, ferroptosis-related genes, and differentially expressed genes (DEGs). As a result, a set of senescence- and ferroptosis-related DEGs (CF-DEGs) was obtained.

Consensus Clustering Analysis of Genes

Unsupervised clustering analysis of gene expression data from AS patients was performed based on the identified CF-DEGs using the “ConsensusClusterPlus” package in R. The optimal k value for the subtypes was determined using consensus matrices, consensus scores, and cumulative distribution function (CDF) analysis. To identify differentially expressed genes between different clusters, the “Limma” package in R was used for differential expression analysis, with parameters set to an absolute value of log FC ≥ 0.5 and P-value < 0.05. Differentially expressed genes between subtypes were obtained, define them as subtypes of DEGs related to ferroptosis and senescence, and volcano plots and heatmaps were generated using the “ggplot2” package. A protein-protein interaction (PPI) network of the differentially expressed genes between subtypes was then constructed based on the STRING database (https://www.stringdb.org), and Hub genes were calculated using the Cytohubba plugin in Cytoscape (3.8.2).

Immune Infiltration Analysis and Enrichment Analysis

The “CIBERSORT” package in R was used to perform immune infiltration analysis based on differentially expressed genes between subtypes. This approach quantified the relative abundances of 22 immune cell types across individuals with distinct immune profiles. Differences in immune cell proportions between subtypes were evaluated using the Wilcoxon rank-sum test, with P < 0.05 considered statistically significant. The results were visualized using the “ggplot2” package. Subsequently, functional enrichment analysis was conducted using the “clusterProfiler” package to investigate potential biological mechanisms. Significantly enriched pathways were identified based on an adjusted P-value < 0.05.

Screening and Validation of Diagnostic Biomarkers

Diagnostic genes for atherosclerosis were screened from ferroptosis- and aging-related genes using eight machine learning methods, including Decision Tree (DT), Extreme Gradient Boosting (XGBoost), C5.0, Neural Network (NNET), K-Nearest Neighbors (KNN), Lasso Regression (LASSO), Support Vector Machine (SVM), and Gradient Boosting Machine (GBM). The machine learning algorithms covered linear models (LASSO), nonlinear models (SVM, NNET), tree-based models (DT, C5.0, XGBoost, GBM), and instance-based algorithms (KNN) to enable a comprehensive comparison of model performance from multiple perspectives. Multiple evaluation metrics were used to assess the model performance. These included the root mean square error (RMSE) of residuals, residual reverse cumulative distribution, and the area under the receiver operating characteristic curve (AUC). A diagnostic model was constructed based on the selected genes using the “rms” package. The predictive ability of the nomogram model was evaluated by calibration curves. Decision curve analysis (DCA) was performed with the “ggDCA” package to assess the clinical utility of the model. The diagnostic value was further evaluated by the receiver operating characteristic (ROC) curve. To further validate the generalization ability of the model, an independent GEO dataset, GSE132651, was used as an external validation set for performance evaluation of the optimal model.

Correlation Analysis Between Hub Genes and Immune Cells

The intersection of Hub genes obtained from PPI analysis and diagnostic genes identified through machine learning was used to obtain key diagnostic genes for AS. The expression levels of these genes in disease and control groups were visualized. Finally, Pearson correlation analysis was performed to reveal the correlation between key diagnostic genes, immune cells, and CF-DEGs. The P-value filter for correlation tests was set to P < 0.01. The strong correlations between Hub genes and immune cells were visualized.

Single-Cell Analysis

The downloaded scRNA-seq dataset included calcified atherosclerotic core plaques and adjacent proximal carotid tissues from three patients who underwent carotid endarterectomy, totaling six tissue samples. The scRNA-seq data were converted into Seurat objects using the “Seurat” package in R software. Quality control was performed by removing cells with a mitochondrial gene proportion greater than 20% and those expressing fewer than 200 or more than 3000 genes. The “Find Variable Features” function was used to identify the top 2000 highly variable genes for principal component analysis (PCA), resulting in the identification of 15 principal components. The “RunUMAP” function was then used for dimensionality reduction and clustering analysis of the Seurat objects. The clustering results were annotated based on literature searchings.²¹ Finally, pseudotime analysis and cell communication analysis were conducted on the annotated cells using the “monocle3” and “CellChat” packages.

Experimental

Animals

The experimental protocol for animal studies was reviewed and approved by the Experimental Animal Ethics Committee of Henan University of Chinese Medicine (Approval Number: IACUC-202408005). All animal studies were conducted in accordance with the UK Animals (Scientific Procedures) Act 1986 and the Guidelines for the Care and Use of Animals by the Chinese Physiological Society. C57BL/6 and ApoE^−/− mice were purchased from Si bei fu (Beijing) Biotechnology Co., Ltd. The mice were randomly divided into two groups: the control group consisted of C56BL/6 mice fed a normal diet for 24 weeks (n = 3); the model group consisted of ApoE^−/− mice fed a high-fat diet for 24 weeks (n = 3). After the experiment, the mice were euthanized by injection of pentobarbital sodium, and their aortas were collected.

Atherosclerotic Lesion Analysis

To analyze the atherosclerotic lesions, the hearts were embedded in OCT to prepare cross-sections of the aortic root. The sections were stained with oil red for morphological analysis. The oil red “O” staining solution was prepared according to the instructions of the Oil Red O Staining Kit (Beyotime, Lot Number: Z907240905), and the staining operation was carried out according to the procedure. All images were obtained using a microscope. The lesion areas in the aortic root sections were quantified using Image-Pro Plus 6.0 for statistical analysis.

qPCR Analysis

The total RNA of the cells was extracted using the RNA extraction solution (Servicebio, Lot Number: CR2406057), and the concentration and purity of the RNA were detected. The cDNA was synthesized using the reverse transcription kit (Vazyme, Lot Number: 7E1981F4), and the real-time quantitative PCR experiment was carried out using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Lot Number: Q711-02). The experiment was repeated three times. GAPDH was used as the internal reference gene. The primer sequences of each gene are shown in Table 1. The relative expression levels of each gene were calculated using the 2-ΔΔCt method.

Table 1 Gene Primer Sequence

Statistical Analysis

Statistical analysis was carried out using GraphPad Prism 9.5 software. An unpaired t-test was used for the comparison between the two groups, and P < 0.05 was considered to indicate a statistically significant difference.

Results

Identification of Differentially Expressed Genes Related to Ferroptosis and Senescence

The study workflow is shown in Figure 1, and the rationale for each method used is explained in the Supplementary Material. First, the “limma” package was used to perform background correction and data normalization on the gene expression data from the GSE10097 dataset (Figure 2a and b). Differential expression analysis was then conducted on the normalized data, resulting in the identification of 2288 DEGs. A volcano plot and heatmap were further generated to visualize these DEGs (Figure 2c and d).

Figure 1 Flow chart of this study.

Figure 2 Screening of Differentially Expressed Genes. (a) Boxplot of the dataset before normalization; (b) Boxplot of the dataset after normalization; (c) Heatmap of Atherosclerosis-related differentially expressed genes. Red indicates the control group; blue indicates the disease group; (d) Volcano plot of Atherosclerosis-related differentially expressed genes.

Construction of Weighted Gene Co-Expression Network

Clustering analysis was performed on the GSE10097 dataset from GEO, and outlier samples were identified and excluded (Figure 3a). The soft threshold β was set to 6 when R² > 0.9 and the mean connectivity was relatively high (Figure 3b). Using the dynamic tree cut method, 11 target modules were identified (Figure 3c and d). Correlation analysis between the target modules and the AS group revealed that the “turquoise” and “yellow” modules were highly correlated with the disease group. A total of 3695 key genes were obtained from these modules (Figure 3e). The intersection of DEGs, WGCNA, and genes related to cellular senescence and ferroptosis yielded 23 CF-DEGs (Figure 3f), The genetic information is provided in Supplementary Table 1.

Figure 3 Weighted Gene Co-expression Network Analysis (WGCNA) (a) Dendrogram of 104 samples after removing 2 outliers; (b) Process of soft-threshold selection. The X-axis shows soft-threshold values, and the Y-axis shows the scale-free topology fit index. The final soft-threshold β was set to 6 for further analysis; (c) The upper part shows the hierarchical clustering tree of genes. The lower part shows gene modules, each represented by a different color; (d) Topological overlap matrix (TOM) of all genes. Light colors indicate low overlap; dark colors indicate high overlap. Darker colors reflect stronger interactions between genes; (e) Identification of Atherosclerosis-related gene modules. Each square shows the correlation (r-value and p-value) between Atherosclerosis and gene modules. Red indicates positive correlation; blue indicates negative correlation; (f) Venn diagram used to identify DEGs related to cellular senescence and ferroptosis.

Identification of Gene Subclusters Based on Ferroptosis- and Senescence-Related Genes

Using the expression profiles of 23 CF-DEGs, unsupervised clustering analysis was performed with a consensus clustering algorithm to identify subtypes associated with AS. The value of k was set from 1 to 9, and the results indicated that k=2 was the optimal parameter, dividing the 69 Atherosclerosis samples into Subtype A and Subtype B (Figure 4a–d). Differential analysis between the A and B subtypes was conducted, resulting in the identification of 421 ferroptosis- and senescence-related subtype DEGs. Volcano plots and heatmaps were generated using the “ggplot2” package to visualize these findings (Figure 4e and f).

Figure 4 Identification of Subclusters Based on Key Genes. (a) Results of consensus clustering analysis. The final number of clusters was set to k=2; (b) Cumulative distribution function (CDF) plot. The Y-axis represents the consensus index, and the X-axis represents the cumulative distribution; (c) and (d) Tracking plots showing classification results under different k values. Each colored block represents a different cluster; (e) Heatmap of DEGs between the two subtypes; (f) Volcano plot of DEGs between the two subtypes.

Immune Infiltration and Enrichment Analysis of DEGs

Immune infiltration analysis of the 421 key genes was performed using the “CIBERSORT” package. The results showed that dendritic cells had the highest proportion in all samples. Eight immune cell types showed significant differences between the two subtypes (P<0.05): activated dendritic cells (p=0.008), macrophages M0 (p=0.0007), resting NK cells (p=0.0088), plasma cells (p=0.04), naive CD4 T cells (p=0.0004), follicular helper T cells (p=0.036), gamma delta T cells (p=0.0028), and regulatory T cells (Tregs) (p=0.0028). These findings suggest that these eight immune cell types play a role in the pathogenesis of AS (Figure 5). GO enrichment analysis of the 421 genes indicated that AS lesions involve processes such as myeloid leukocyte activation, phagocytosis, leukocyte-mediated immunity, leukocyte migration, positive regulation of cytokine production, positive regulation of tumor necrosis factor superfamily cytokine production, secretory granule membrane, endocytic vesicle, specific granule, lysosomal lumen, secretory granule lumen, low-density lipoprotein particle, pattern recognition receptor activity, integrin binding, immune receptor activity, proteoglycan binding, glycosaminoglycan binding, and structural constituent of muscle (Figure 6). KEGG enrichment analysis identified pathways such as lysosome, Toll-like receptor signaling pathway, chemokine signaling pathway, leukocyte transendothelial migration, lipid and atherosclerosis, neutrophil extracellular trap formation, PPAR signaling pathway, Type I diabetes mellitus, efferocytosis, apoptosis, and NF-kappa B signaling pathway (Figure 6). GSEA analysis revealed that these gene sets are associated with the upregulation of the cellular senescence pathway and the NOD-like receptor signaling pathway, and the downregulation of cholesterol metabolism and ferroptosis pathways.

Figure 5 Immune Infiltration Analysis. (a) Immune cell infiltration percentages in each sample were estimated using the CIBERSORT algorithm. The X-axis shows different immune cell types, and the Y-axis represents the infiltration proportion; (b) Box plots comparing immune cell expression differences between Cluster A and Cluster B.

Figure 6 Gene Enrichment Analysis. (a) GO enrichment analysis of subtype DEGs, including Biological Process (BP), Cellular Component (CC), and Molecular Function (MF); (b) KEGG enrichment analysis of subtype DEGs; (c) GO enrichment circular plot. Blue dots represent downregulated genes; red dots represent upregulated genes; (d) GSEA plot of subtype DEGs. The top line graph shows the enrichment score of gene sets, where the peak indicates the location of core genes. The middle panel displays the rank positions of genes in the functional annotation gene sets, marked by vertical lines. The bottom area plot shows the distribution and signal-to-noise ratio of all genes.

Construction of PPI Network and Development and Validation of Machine Learning Models

To identify potential interactions among the 421 ferroptosis- and senescence-related subtype DEGs, the genes were imported into STRING for PPI analysis, and the data were further analyzed using Cytoscape software. Using the Cytohubba plugin, 20 key genes were identified: IL1B, CCL3, TLR1, ITGAX, TLR7, CXCR4, IRF8, FCGR3A, MMP9, CCR5, CCL4, ITGB2, CCR1, CD68, ITGAM, CD86, TYROBP, IL18, SPI1, and TLR2 (Figure 7a and b).Next, eight machine learning algorithms were applied to the 421 ferroptosis- and senescence-related subtype DEGs to identify diagnostic genes. LASSO regression analysis was found to be the most accurate in predicting the onset of AS (Figure 7c–e). A VENN diagram was created by intersecting LASSO analysis and PPI genes, resulting in the identification of two diagnostic biomarkers: IL1B and CCL4 (Figure 7f and g). An AS diagnostic nomogram was then constructed using the “rms” package (Figure 8a). The predictive ability of the model was assessed using a calibration curve, which showed minimal difference between the predicted and actual risk of disease, indicating that the nomogram model accurately predicts AS onset (Figure 8b). DCA analysis showed that the IL1B and CCL4 curves provided significant clinical benefit for patients when the risk threshold was between 0 and 0.75 (Figure 8c), and the validation set analysis supported these findings (Figure 8d). ROC analysis revealed that both IL1B (AUC=0.931), CCL4 (AUC=0.825) had AUC values above 0.6 (Figure 8e and f). The validation set results were consistent with those of the training set (Figure 8g–h), indicating a high predictive accuracy of these genes for AS onset. This further demonstrates that the nomogram model has good diagnostic capability for early-stage AS patients. The information for each machine learning model is provided in Supplementary Table 2.

Figure 7 Screening Process of Hub Genes. (a and b) Construction of the PPI network and identification of core genes; (c) Distribution of prediction residuals for each model. Red dots represent the Root Mean Square Error (RMSE), reflecting the overall model fitting error; (d) Cumulative distribution of prediction errors. A steeper and more left-shifted curve indicates lower error and more stable performance; (e) ROC curves showing the classification performance of each model, with AUC values labeled; (f) LASSO analysis identified 20 key genes; (g) Venn diagram used to screen and identify 2 hub genes.

Figure 8 Model Construction and Validation; (a) Nomogram for diagnosing atherosclerosis based on IL1B and CCL4, *** indicates that the variable is significantly associated with disease prediction in the model (p < 0.001); (b) Calibration curve. The X-axis shows predicted probability, and the Y-axis shows actual probability. The dashed line represents the ideal model; the closer the solid line is to the dashed line, the more reliable the model; (c and d) Decision Curve Analysis (DCA) of hub genes in the training and validation sets. The X-axis shows the threshold probability, and the Y-axis shows the net benefit; (e) ROC curve analysis of IL1B in the training set. The X-axis represents specificity, and the Y-axis represents sensitivity; (f) ROC curve analysis of CCL4 in the training set; (g) ROC curve analysis of IL1B in the validation set; (h) ROC curve analysis of CCL4 in the validation set.

Correlation Analysis Between Core Genes and Immune Cells

Further analysis of the Hub genes revealed that both IL1B and CCL4 were highly expressed in the disease group (Figure 9a and b). Correlation analysis between significantly different immune infiltrating cells and Hub genes showed that IL1B was significantly positively correlated with follicular helper T cells and significantly negatively correlated with gamma delta T cells. CCL4 was significantly positively correlated with follicular helper T cells and significantly negatively correlated with regulatory T cells (Tregs) (Figure 9c and d). Correlation analysis between Hub genes and CF-DEGs showed that IL1B was highly positively correlated with PTPN6 (r=0.72) and highly negatively correlated with ZEB1 (r=−0.64). Similarly, CCL4 was highly positively correlated with PTPN6 (r=0.89) and highly negatively correlated with ZEB1 (r=−0.74) (Figure 9e–i).

Figure 9 Correlation Analysis Between Hub Genes and Immune Cells. (a) Box plot showing the expression level of CCL4 in the disease group and control group; (b) Box plot showing the expression level of IL1B in the disease group and control group; (c) Correlation analysis between CCL4 and immune cells; (d) Correlation analysis between IL1B and immune cells; (e–i) Correlation analysis between hub genes and CF-DEGs. Each square displays the correlation coefficient (r-value). Blue indicates a positive correlation, and red indicates a negative correlation. Darker colors represent stronger correlations, Cross mark represents p > 0.05, indicating no statistical significance.

Single-Cell Analysis Results

Since the previous analyses were based on transcriptomics and did not describe the cellular diversity of the AS vascular microenvironment at the single-cell level, we analyzed the importance of Hub genes using published single-cell transcriptome data (GSE159677). After constructing the Seurat object from the dataset, the normalized expression profiles were displayed (Figure 10a), showing a positive correlation between nCountRNA and nFeatureRNA (Figure 10b). Principal component analysis and clustering analysis of the dataset identified 15 clusters (Figure 10c and d). After annotation using the “SingleR” package, six different cell types were identified: Endothelial Cells (EC), Fibroblasts, Macrophages, Monocytes, T cells, and Vascular Smooth Muscle Cells (VSMC) (Figure 10e). The study also found that IL1B and CCL4 were highly expressed in the disease group, consistent with the transcriptome sequencing data analysis (Figure 10f). Further analysis of IL1B and CCL4 expression across the 12 cell types showed that IL1B was highly expressed in Macrophages and Monocytes, while CCL4 was significantly expressed in Macrophages, Monocytes, and T cells (Figure 10g–h). Cell communication analysis revealed important interactions between immune cells and stromal cells during AS development, suggesting mutual regulation between these cells and the AS vascular microenvironment (Figure 11a–c). Next, pseudotime analysis was conducted using the “monocle3” package to map the developmental trajectories of single cells (Figure 11d and e). Arrows indicate the starting points of the timelines for each cell group. We observed changes in IL1B and CCL4 expression during the dynamic transition of macrophages and monocytes, suggesting that these genes may be involved in the fate determination of these cell groups (Figure 11f).

Figure 10 Analysis of Single-Cell RNA Sequencing Data. (a)Quality control metrics of single-cell RNA-seq data, including the number of detected genes per cell (nFeature_RNA), the number of unique molecular identifiers (UMIs) per cell (nCount_RNA), and the percentage of mitochondrial gene expression per cell (percent.mt); (b) Correlation analysis between nCount_RNA and percent.mt, as well as between nCount_RNA and nFeature_RNA; (c) Elbow plot of principal component analysis (PCA) indicating the optimal number of principal components for downstream analysis; (d) Cells were clustered into 14 distinct groups; (e) Cells were visualized and grouped into six major cell types using Uniform Manifold Approximation and Projection (UMAP); each color represents a different cell type cluster; (f) Expression distributions of IL1B and CCL4 genes in disease and control datasets; (g) Spatial distribution of IL1B and CCL4 gene expression across the identified clusters; (h) Expression patterns of IL1B and CCL4 across different cell types.

Figure 11 Cell chat and pseudo temporal analysis. (a) Cell-cell communication network diagram; thicker lines indicate stronger interaction weights between cell types; (b) Detailed intercellular communication network among six cell subpopulations; the thickness of the lines represents the strength of interaction between cell types; (c) Hierarchical communication network of the CCL signaling pathway; thicker lines denote stronger communication intensity via the CCL signaling pathway between cell types; (d) Pseudotime trajectory plot generated using the Monocle3 algorithm. Light gray circles represent distinct terminal states along the developmental trajectory, while black circles indicate branch points, suggesting potential divergent cell fate decisions, The numbers inside the circles denote pseudotime values, where smaller numbers correspond to earlier differentiation stages and larger numbers correspond to later stages; (e) Trajectory of macrophage and monocyte differentiation into VSMC cells predicted by the Monocle3 algorithm, Numbers indicate pseudotime, with smaller values representing earlier differentiation stages and larger values indicating later stages; (f) Dynamic expression changes of IL1B and CCL4 genes across different cell types.

Experimental Verification of the Expression of Hub Genes in the Animal Model

In this study, we evaluated the atherosclerotic lesions in APOE^−/− mice by Oil Red O staining, and quantified the ratio of the plaque area to the total cross-sectional area of the aortic root. The results showed that, compared with the control group, the plaque area in the aortic root of mice in the model group increased significantly, as shown in Figure 12a and b. The results of the qPCR experiment demonstrated that, compared with the control group, the expression levels of IL1B and CCL4 genes in the model group were significantly upregulated (P < 0.05), as shown in Figure 12c and d). These results are consistent with the findings of bioinformatics analysis.

Figure 12 Results of experimental verification. (a and b) Lesions in aortic root cross sections were determined by Oil Red O staining; (c) The relative mRNA expressions of IL1B; (d) The relative mRNA expressions of CCL4. Data are presented as mean ± SEM, *p < 0.05, **p < 0.01.

Discussion

Currently, the number of AS patients worldwide is increasing annually, with a trend toward younger age groups, and the complications it causes place a significant burden on patients.²² The pathogenesis of atherosclerosis is complex and involves multiple cell types, including endothelial cells (ECs), vascular smooth muscle cells (SMCs), adventitial fibroblasts, macrophages, and other immune cells. Key factors in the development of atherosclerosis include endothelial dysfunction, leukocyte adhesion, foam cell formation, and SMC phenotypic transformation.^23,24 Notably, the pathological progression of AS is closely linked to cellular senescence and ferroptosis. Cellular senescence is a cell state triggered by endogenous or exogenous stimuli, characterized by stable cell cycle arrest and a complex senescence-associated secretory phenotype (SASP). Once senescent cells accumulate in tissues, they can accelerate the progression of age-related diseases such as atherosclerosis, osteoarthritis, and cancer.²⁵ Senescent endothelial cells reduce NO production and impair the endothelial barrier function due to disruptions in intercellular adhesion and tight junctions in AS.²⁶ The extent of endothelial coverage in AS lesions is a key protective factor for plaque stability, and the erosion of endothelial cells can lead to atherosclerotic thrombosis on the lesion surface. Additionally, increased endothelial inflammation is considered a harmful process in the chronic inflammatory environment of arteries. Senescent endothelial cells can exacerbate the inflammatory response by increasing the expression of adhesion molecules such as VCAM-1 and ICAM-1, which enhances the recruitment of leukocytes.^27,28 Ferroptosis is a distinct form of iron-dependent, lipid peroxidation-driven programmed cell death. In recent years, it has increasingly been recognized as an important process mediating the onset and progression of various cardiovascular diseases, including atherosclerosis, myocardial ischemia-reperfusion injury, and septic cardiomyopathy.²⁹ Furthermore, AS patients often exhibit changes in peripheral blood immune components during disease onset, reflecting the disease’s severity. Therefore, exploring immune-related biomarkers and potential therapeutic targets is a research focus in this field.¹⁵

This study identified 23 ferroptosis- and senescence-related DEGs. Based on these genes, unsupervised clustering analysis was performed using a consensus clustering algorithm to identify AS-related subtypes, resulting in 421 subtype DEGs. Immune infiltration analysis of these 421 DEGs revealed differences in eight immune cell types between the two subtypes: activated dendritic cells, Macrophages M0, resting NK cells, plasma cells, naive CD4 T cells, follicular helper T cells, gamma delta T cells, and regulatory T cells (Tregs). Previous studies have shown that both the innate and adaptive immune systems play key roles in driving AS-related chronic inflammation, with significant heterogeneity among leukocyte subtypes in the arterial wall, which regulate inflammation in atherosclerosis formation.³⁰ This study also identified two diagnostic biomarkers through machine learning methods based on the 421 ferroptosis- and senescence-related subtype DEGs, providing insights into AS-related immune genes. Currently, many immune infiltration algorithms have revealed differences in immune cell infiltration between disease groups and controls in AS onset. Jing Wang³¹ using bioinformatics methods, studied the role of immune cells in the formation of unstable plaques and explored diagnostic biomarkers. Their findings indicated that M1 macrophages are an important factor in unstable plaque formation, and CD68, PAM, and IGFBP6 can serve as diagnostic markers for identifying unstable plaques. In future research, predictive biomarkers should not only involve the composition of immune infiltration and the characteristics of the inflammatory response but also reveal the heterogeneity of immune composition, spatial distribution, and function.

GO analysis of the differentially expressed genes indicated that AS lesions involve biological processes such as myeloid leukocyte activation, phagocytosis, leukocyte-mediated immunity, leukocyte migration, positive regulation of cytokine production, positive regulation of tumor necrosis factor superfamily cytokine production, pattern recognition receptor activity, and immune receptor activity. KEGG enrichment analysis showed that pathways such as the Toll-like receptor signaling pathway, chemokine signaling pathway, leukocyte transendothelial migration, lipid and atherosclerosis, PPAR signaling pathway, apoptosis, and NF-κB signaling pathway are associated with AS onset. This suggests that the mechanisms of AS are closely related to biological processes such as fatty acid metabolism, inflammatory response, and immune regulation. For example, PPARs are a class of ligand-activated transcription factors that regulate genes related to lipid metabolism. PPARα, one of its subtypes, is highly expressed in the heart and vascular walls and regulates target genes closely related to lipid metabolism and inflammation. Activation of the PPARα pathway can reduce autophagy-dependent ferroptosis caused by mitochondrial DNA damage and alleviate hyperlipidemia-induced vascular calcification.³² Studies have also shown that senescent immune cells in AS undergo cell cycle arrest, morphological changes, and phenotypic alterations in their abundance and secretory profiles, including cytokines, chemokines, matrix metalloproteinases, and Toll-like receptor expression. Currently, the clearance of senescent cells is considered a key target for preventing or treating AS.³³ GSEA analysis showed that these gene sets are associated with the upregulation of the cellular senescence pathway and the NOD-like receptor signaling pathway, and the downregulation of cholesterol metabolism and ferroptosis pathways. Since the production of free radicals, increased fatty acid supply, and enhanced lipid peroxidation are key inducers of ferroptosis, they may influence ferroptosis by regulating lipid metabolism and fatty acid metabolism, thereby mediating the onset of AS.^10,34 Senescence is associated with elevated levels of inflammatory cytokines in the aorta, and impaired vascular mitochondrial function accelerates the formation of AS. Of course, this requires further study.³⁵ Research has found that macrophages, vascular endothelial cells, and vascular smooth muscle cells are in a senescent state at atherosclerotic lesion sites. Additionally, bacterial infections and the accumulation of lipopolysaccharides have been observed in AS patients, along with a SASP. Therefore, some studies suggest that LPS derived from Gram-negative bacteria may exacerbate AS lesions by inducing and enhancing the senescence of specific vascular cells and SASP-associated inflammatory characteristics in AS lesions.³⁶

In recent years, machine learning has been widely applied in the clinical field and is considered an important tool in healthcare.³⁷ This study identified IL1B and CCL4 as ferroptosis- and senescence-related biomarkers for AS, and validation showed that these biomarkers have good diagnostic value for AS. Previous studies have confirmed that IL-1β is associated with acute and chronic inflammation and plays an immunoregulatory role in the development of AS. Anti-inflammatory treatment with canakinumab, targeting the IL-1β innate immune pathway, has been shown to significantly reduce the recurrence rate of cardiovascular events compared to placebo.³⁸ In APOE^−/− mice, upregulation of IL-1β expression increased the area of atherosclerotic lesions in the ascending aorta and aortic root.³⁹ However, it is important to note that most previous studies on IL-1β have focused on its regulation of immune or inflammatory cells. Recent research has found that IL-1β can regulate iron-sulfur cluster homeostasis by inducing acetylation of the mitochondrial inner membrane protein NNT in tumor cells, thereby inhibiting ferroptosis and mediating resistance to immunotherapy.⁴⁰ This finding provides insights into potential AS treatments. Additionally, studies have shown that inhibiting the secretion of IL-1β can improve doxorubicin-induced senescence in cardiac fibroblasts and improve poor prognosis in heart disease.⁴¹ IL-1β, as a key biomarker of vascular calcification, may also be involved in the senescence of VSMCs. Linzi Han and colleagues found that IL-1β induces VSMC senescence and promotes VSMC calcification by activating the NF-κB/p53/p21 signaling pathway, which may be associated with AS plaque formation.⁴² CCL4, also known as macrophage inflammatory protein MIP-1β, belongs to the CC chemokine family and can be secreted by monocytes, B cells, T cells, NK cells, dendritic cells, neutrophils, fibroblasts, endothelial cells, and epithelial cells.⁴³ Elevated levels of CCL4 have been reported in the peripheral blood of AS patients. Inhibiting CCL4 levels can reduce the activity of metalloproteinases-2 and −9 in macrophages, suppress the production of TNF-α and IL-6, and decrease the activation of endothelial cells and macrophages, suggesting that CCL4 has potential as a new therapeutic target for AS.⁴⁴ Some studies hypothesize that CCL4 may also play a crucial role in aging-related vascular dysfunction. CCL4 promotes inflammation and cellular senescence by stimulating ROS production, leading to impaired cell function, and may be a potential therapeutic target for vascular protection during aging.⁴⁵ Studies have shown that iron levels significantly affect the expression of MIP-1β (CCL4) in murine macrophages. Specifically, ferrous sulfate suppresses CCL4 expression, while the iron chelator desferrioxamine promotes its upregulation. This suggests that CCL4 expression may be regulated by iron homeostasis and potentially linked to the ferroptosis pathway. Therefore, it is plausible that CCL4, as an immune-inflammatory chemokine, may contribute to ferroptosis-related inflammatory responses in atherosclerosis by modulating macrophage activity.⁴⁶ Through the qPCR experiment, we found that the expression levels of IL1B mRNA and CCL4 mRNA were significantly upregulated in the model group. Therefore, it is speculated that IL1B and CCL4 have the potential to be used as biomarkers for predicting AS. At present, there are few studies related to CCL4 and ferroptosis. In the future, the therapeutic strategies for AS can be explored by studying the association between ferroptosis and immune cells during the onset of AS.

Correlation analysis between Hub genes and immune cells revealed that the IL1B gene is significantly associated with T cells.follicular.helper and T cells.gamma.delta, while the CCL4 gene is significantly associated with T cells.follicular.helper and T cells.regulatory (Tregs). T cells.follicular.helper (Tfh) have recently been defined as a new subset of CD4+ T cells, primarily expressing surface molecules such as CXC chemokine receptor 5 and inducible costimulatory molecule, and secreting the cytokine interleukin-21. The key transcription factor for Tfh cells is B-cell lymphoma 6. Tfh cells are found in the germinal centers of lymphoid tissues and peripheral blood, where they promote B-cell maturation and differentiation, germinal center formation, and antibody production, playing an important role in the pathogenesis of various autoimmune diseases.⁴⁷ T cells.gamma.delta (γδT cells) are lymphocytes with multiple functions in innate and adaptive immune responses, pathogen defense, antigen presentation, and inflammation regulation, participating in the early development of atherosclerosis.⁴⁸ Duc M Vu revealed the pathogenic role of γδT cells in early AS formation through animal experiments, with a mechanism that may involve IL-17 production and neutrophil induction. Therefore, γδT cells are promising targets for AS intervention, potentially reducing plaque accumulation and inhibiting inflammatory responses.⁴⁹ Studies have shown that T cells.regulatory (Tregs) account for 5–10% of the total CD4+ T lymphocytes in peripheral blood and have been proven to have a protective effect against atherosclerosis, making them a new target for cardiovascular disease and atherosclerosis.⁵⁰ Monika Sharma⁵¹ also found that in several mouse models of AS regression, an increase in Tregs was a common feature of plaque regression. In summary, the Hub genes are closely associated with these immune cells, leading us to speculate that CCL4 and IL1B may influence the onset and progression of AS by regulating these cells. Correlation analysis between CCL4 and IL1B with CF-DEGs revealed that these genes are highly correlated with PTPN6 and ZEB1.

The emergence of single-cell biology has opened a new chapter in understanding the biological processes of disease development and in diagnostics, monitoring, and treatment. This technology can identify novel cell populations that play critical roles in the progression of AS and assist in exploring new therapeutic strategies.⁵² The infiltration of immune cells is closely related to the progression of AS and its response to immunotherapy. Studies using single-cell sequencing technology have explored the immune heterogeneity of atherosclerosis and identified interferon-induced CD8+ T cells and macrophage subpopulations as cells that influence AS progression and poor prognosis, which is of significant value for developing new AS immunotherapies.⁵³ In this study, 6 different cell types were identified through single-cell analysis, including ECs, Fibroblasts, Macrophages, Monocytes, T cells, and VSMC, providing a deeper understanding of the cellular diversity within atherosclerotic plaques and surrounding tissues. It was also found that the Hub gene IL1B is highly expressed in Macrophages and Monocytes, while CCL4 is significantly expressed in Macrophages, Monocytes, and T cells. Additionally, this study visualized the communication between different types of cells, confirming a close relationship between immune cells and stromal cells in the onset and progression of AS. Using the “Monocle3” package, the pseudotime algorithm was applied to single-cell RNA-Seq data to order and arrange individual cells based on their states during differentiation and other biological processes, revealing the connection between gene expression changes and cell fate determination.⁵⁴ Through this method, it was observed that IL1B and CCL4 expression levels change during the dynamic transition of macrophages and monocytes. These results suggest that Hub genes may influence the onset and progression of AS by regulating immune cell function.

This study has several limitations. First, the available AS-related RNA-seq and scRNA-seq datasets are limited, with small sample sizes. In the future, larger AS datasets will need to be analyzed to validate the findings derived from GSE10097, GSE132651, and GSE159677. Second, it should be noted that the key effectors involved in ferroptosis and cellular senescence pathways are mostly regulated at the post-translational level. Therefore, the post-translational levels of these molecules are crucial for maintaining cellular function and should be considered. This study is limited by the lack of corresponding proteomics data, as we only screened ferroptosis- and cellular senescence-related genes based on RNA-seq datasets. Additionally, at the beginning of the study, important diagnostic genes were screened solely from candidate genes related to ferroptosis and senescence. Consequently, only genes with significant diagnostic value were retained to establish the predictive model, while genes closely related to ferroptosis and senescence but with lower diagnostic value were excluded.

Conclusions

In conclusion, this study has provided new diagnostic targets for predicting the onset of AS. The diagnostic values of IL1B and CCL4 for AS have also been confirmed through the combination of experimental verification. The study also predicts characteristic genes related to cellular senescence and ferroptosis in the development of AS, which holds significance for understanding the mechanisms of AS and exploring therapeutic strategies. Future research should validate the clinical applicability of these diagnostic biomarkers and investigate the roles of IL1B and CCL4 in the development of AS, assessing their potential as biomarkers and therapeutic targets for AS.