Introduction

Atherosclerosis (AS) is a chronic inflammatory disease marked by disruptions in lipid metabolism and immune response, resulting from lipid-laden macrophages accumulating in arterial walls.1 This disease is characterized by the formation of plaques within arterial walls, which can lead to life-threatening events such as heart attacks, and strokes.2 Stroke, a consequence of AS with high mortality and disability rates worldwide, affects 30 million people and poses a significant threat to human health. Approximately 795,000 people in the United States experience a new or recurrent stroke annually, 87% of which are ischemic strokes.3 Ischemic strokes, are caused by advanced atherosclerotic plaque (AAS), surface erosion, and rupture, often without warning,4 lead to significant adverse ischemic events (IE). Hence, early detection and prediction of AS plaque instability are crucial. With the deepening of AS-related research, we noticed that the biological properties of macrophages are key to determining lesion size, composition, and instability in AS plaque.5 There are significant differences in immune cell infiltration between early atherosclerotic plaque (EAS) and advanced atherosclerotic plaque. Macrophages in AS plaques demonstrate significant variability and adaptability, modulating their microenvironment through mechanisms such as elevated lipid concentrations, deregulated cytokine functions, hypoxia, apoptosis, and necroptosis.6 Furthermore, oxidative stress and inflammation are interconnected processes that establish a strong feed-forward cycle, accelerating the progression of atherosclerotic plaque.7 A previous study has indicated that excessive mitochondrial reactive oxygen species (mtROS), a byproduct of oxidative stress, inflict damage on mitochondrial DNA (mtDNA), proteins, and lipids. These damages have been observed to escalate during AS progression, influencing whether the system activates or remains unresponsive and ultimately affecting the stability of AS plaques.8 Additionally, the risk of plaque rupture is significantly influenced by inflammatory responses, which are driven by interactions among lipoproteins, monocytes, and various macrophage subtypes. It’s crucial to acknowledge a significant dilemma in current AS research, namely, the evident heterogeneity within the atherosclerotic plaque microenvironment. This heterogeneity primarily stems from the diverse biological functions of macrophage subtypes.9 Simply classifying them as M1/M2 types cannot adequately elucidate the mechanism of action of macrophage subtypes in AAS and is likely to produce deviation.10 Therefore, it is urgent to conduct in-depth research on more effective classification characteristics of AS plaque subtypes and macrophage subtypes, in order to accurately identify patients who will suffer from adverse ischemic events.

Conventional imaging techniques, focusing on plaque appearance and size, are insufficient for predicting rupture risk and acute thrombotic events.11 Macrophage heterogeneity within AS plaques has become a focal point of recent research.12 High-throughput sequencing technologies such as transcriptome analysis (bulk-seq) and single-cell RNA sequencing (scRNA-seq) have provided deep insights into the complexity and dynamic nature of atherosclerotic lesions. In contrast to bulk RNA-seq, which focuses primarily on identifying differentially expressed genes (DEGs) in different groups, scRNA-seq offers diverse perspectives on heterogeneity, cellular interactions, and transcriptome changes, crucial for understanding disease progression. It also enables the identification of distinct cell populations, including various immune cell subtypes within the disease, and offers an in-depth understanding of their roles in disease progression.13–15 However, scRNA-seq’s high cost, small sample size, and limited clinical follow-up information may bias results. Bulk-seq complements single-cell analysis, providing a comprehensive view of gene expression patterns across large cell populations within atherosclerotic tissue. Integrating data from both single-cell and bulk-seq allows for a more holistic understanding of the pathological processes in AS. It also aids in identifying key biomarkers and therapeutic targets for predictive diagnostics, targeted prevention, and personalized treatment strategies, in line with the Predictive, Preventive, and Personalized Medicine (PPPM) model. Partial scRNA-seq analysis has elucidated the transcriptomic profiles of immune cell subpopulations in AS.16,17 However, the prognostic significance of macrophage subtypes and their cellular interactions during AS progression, especially in unstable plaque subtypes, remain poorly understood. Therefore, examining the metabolic features and phenotypic changes in dysregulated macrophage subpopulations linked to AS plaque subtypes, along with identifying more robust diagnostic biomarkers, holds significant clinical importance. Such research could provide a potential treatment for AS by targeting the plaque macrophage microenvironment strategy.

Applying machine learning algorithms, such as weighted correlation network analysis (WGCNA), least absolute shrinkage and selection operator (LASSO), support vector machine (SVM-RFE), and random forest algorithm (RF), to screen for signature genes can further elucidate the characteristic distribution links between macrophage-related biomarkers and inflammatory immune cells, particularly in macrophage subtypes.18 Additionally, a diagnostic model based on these biomarkers has been developed to detect unstable plaques and mitigate AS. The efficacy of this model was assessed using ten machine learning algorithms. This research aims to identify a new characteristic of AS plaque subtype and enhance our understanding of macrophage subtypes in AS, offering insights into their distinct roles and origins, finding reliable biomarkers for diagnosing AS plaque progression, and revealing new molecular mechanisms for AS detection and treatment.

Materials and Methods Data Collection

The gene expression profiling data and corresponding clinical information for AS were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. Seven GSE datasets (GSE159677, GSE155512, GSE21545, GSE41571, GSE43292, GSE28829, and GSE163154) were selected for further analysis. scRNA-seq data for AS were obtained from GSE159677 and GSE155512, including six carotid atherosclerotic samples and three proximally adjacent vascular samples. The “Seurat” package was utilized to integrate samples from these two datasets, both derived from Illumina sequencing platforms. GSE21545 comprised 223 AS samples with corresponding prognosis information, while GSE41571 included 11 AS samples, with prognosis information available for six samples. GSE28829 included 13 EAS and 16 AAS plaque samples, while GSE43292 comprised 32 EAS and 32 AAS plaque samples. Due to methodological and design similarities between GSE21545 and GSE41571, as well as between GSE28829 and GSE43292, we merged these two pairs of datasets separately for expression analysis. Batch effects were directly adjusted using the combat function of the “sva” package. GSE163154 contained 16 EAS and 27 AAS plaque samples. Probe expression matrices were converted to gene expression matrices using platform annotation files for GSE21545, GSE41571, GSE43292, GSE28829, and GSE163154 datasets. Principal component analysis (PCA) was conducted to assess the effectiveness of the batch effect correction using the combat function. Furthermore, 585 immune-related marker genes (IRGs) were obtained from a previously published article,19Supplementary Table 1 provides detailed characteristics of the seven datasets.

The Landscape of Immune Cell Infiltrations and Analysis in Bulk Dataset

Single-sample gene set enrichment analysis (ssGSEA) was conducted using the “GSVA” package to quantify the relative infiltration scores of 28 types of immune cells in the GSE21545−41571 dataset.20,21 Kaplan-Meier survival curve were employed to analyze the differences in ischemic event-free survival among AS groups with high and low macrophage infiltration scores. Additionally, the hallmark gene set (h.all.v7.5.1.symbols.gmt) was obtained from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb/index.jsp). Spearman correlation was utilized to assess the relationship between macrophage infiltration scores and various hallmark gene sets. The “limma” package was used to identify DEGs between high and low macrophage infiltration score groups in AS, with a filter criterion of |logFC| > 0.4 and adjusted p-value < 0.05. Visualization was performed using the “ggstatsplot”, “IOBR”, and “ggplot2” packages.

Identification of AS Plaque Subtypes and It’s Biofunction

Based on the expression profiles of 130 DEIRGs in the GSE21545-41571 dataset, we used unsupervised k-means clustering analysis in the “ConsensusClusterPlus” package to divide AS plaques into three clusters. The number of repetitions was set to 1000 to ensure the stability of clustering. Survival analyses and PCA were performed to evaluate prognostic values and the distribution within three clusters. Meanwhile, we performed an ssGSEA algorithm to quantify the infiltration scores of the various immune cells among three clusters, and the hallmark gene set was used to assess the biological function differences among three clusters. Furthermore, Gene Set Variation Analysis (GSVA) was applied using the “GSVA” and “limma” packages to evaluate the pathway activities in these clusters (adjusted p-value < 0.05 and |t-value| > 1 were selected). The gene expression of immune checkpoints was quantified between these clusters. Immune checkpoints were collected from the literature.22

Weighted Gene Co-Expression Network Analysis

WGCNA was used to construct highly correlated co-expressed gene modules associated with cluster 3 plaque subtype. Briefly, the adjacency matrix was established based on the soft power value ß and the similarity matrix was calculated by Pearson correlation analysis between each gene. Then, with a threshold soft power of 8 in the GSE21545-41571 dataset, the adjacency matrix was transformed into a topological overlap matrix (TOM) and the corresponding dissimilarity (1-TOM). Next, genes with similar expression profiles were classified into modules to further build a hierarchical clustering dendrogram (minModuleSize = 80, mergeCutHeight = 0.25). Finally, we calculated the gene significance (GS) and module membership (MM) through WGCNA. This study focused on gene modules significantly associated with the cluster 3 plaque subtype.

Define DEMGs and Enrichment Analysis

Differentially Expressed Macrophage-Related Genes (DEMGs) were identified through the intersection of WGCNA gene modules correlated with the cluster 3 plaque subtype and the previously identified 130 DEIRGs. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were conducted on the DEMGs. Enrichment analyses were performed based on the “clusterProfiler” package with a p-value < 0.05 as the cutoff criterion. Visualization was achieved using the “circlize”, “ComplexHeatmap”, and “ggplot2” packages.

ScRNA-Seq Data Processing

The “Seurat” package was used to analyze the scRNA-seq data. First, quality control was conducted by filtering out cells satisfying the following criteria: cells with > 200 and < 4000, and a percentage of mitochondrial and ribosomal genes < 10%. Next, a global-scaling normalization method (“LogNormalize”) was applied to equalize the total gene expression in each cell, with a set scale factor of 10,000. The top 2000 variably expressed genes were returned for downstream analysis using the “vst” method. The ScaleData and RunPCA functions to obtain the number of principal components (PCs), and the top 30 PCs were further used for downstream integration. Afterwards, the data were processed by homogenization using the “harmony” package.23 Clustering analysis was conducted utilizing edge weights between individual cells, generating a shared nearest-neighbor graph via the Louvain algorithm, seamlessly integrated into the FindNeighbors and FindClusters functions. The cell cluster/sub-cluster analysis with the resolution set to 0.4 and 1. The identified clusters were visualized using the Uniform Manifold Approximation and Projection (UMAP) method. The cell clusters were labeled as different major cell types based on DEGs and classical marker genes. For sub-clustering analysis, we performed the second-round PCA reduction, harmony homogenization and UMAP projection separately on cells within myeloid cell. The second-round dimension reduction and clustering revealed 10 distinct myeloid subtypes. Cell types were visualized using UMAP plots, displayed by the “scRNAtoolVis” package. We obtained DEGs of different macrophages subtypes using the FindAllMarkers function, and the Wilcoxon rank sum test algorithm was used to calculate the marker genes (logFC > 0.25; p-value < 0.05; min.pct > 0.1). Furthermore, we use the CellMarker2 database (http://bio-bigdata.hrbmu.edu.cn/CellMarker/index.html) and the “SingleR” package to assist with annotation cell types. The “ggplot2”, “ggalluvial”, and “pheatmap” package were used to display the cellular composition and correlation between normal and AS groups.

Evolutionary Trajectory Analysis

To investigate dynamic biological processes, such as interconversion and evolutionary trajectories of macrophages. We first extracted macrophage subtype objects for trajectory analysis using the “Monocle2” package. 2500 significantly pseudotime-dependent genes were identified by the differentialGeneTest function and used for cell ordering. Then we used the “DDRTree” method to reduce cell dimensions and calculated the cell differentiation state types using the reduceDimension function. Finally, we used the plotCellTrajectory function to visualize the differentiation trajectory of DEMGs score/macrophage subtypes/ in the three-state differentiation process (Pre-branch, Cell fate 1, Cell fate 2) for the macrophage subtypes, and the plot_pseudotime_heatmap function was used to visualize the expression patterns of 43 DEMGs in different macrophage subtypes. Branched Expression Analysis Modeling (BEAM) analysis was used to identify DEGs associated with differentiated macrophage states (q-value < 1e-200 were selected), which were then visualized using the plot genes branched heatmap function. Based on the enrichment scores of Cell fate 2’s DEGs calculated using the “GSVA” package, we conducted survival analysis. Furthermore, we conducted enrichment analysis to clarify the biological functions of DEGs associated with differentiated macrophage states (p-value < 0.05). CytoTRACE scores were calculated using the “CytoTRACE” package (version 0.3.5) to verify the linear transitions inferred by Monocle. The CytoTRACE algorithm can infer developmental initiation stages of macrophage subtypes from scRNA-Seq data based on stemness scores.24 The “ggridges” package was used to analyze the distribution of macrophage subtypes across different groups/states along the pseudotime axis.

Identification of DEMGs Score and Pathways Activities in Different Cell Types/States

For scRNA-seq data, we calculated DEMGs enrichment scores at a single-cell resolution using the AddModuleScore function in the “Seurat” package. Additionally, the “ggpubr” package was used for visualizing DEMGs scoring across different states generated with Monocle. We also assessed whether fifty hallmark gene sets were enriched in myeloid subtypes using the RRA method of the “irGSEA” package. Signature genes of hypoxia, angiogenesis, cholesterol homeostasis, and fatty acid metabolism were obtained from fifty hallmark gene sets in MSigDB. Inflammatory genes included TNF, CCL2, CCL3, CCL4, CXCL10, S100A8, S100A9, CXCL1, and MPO. And then, the AUCell algorithm with default settings was used to infer these three AS-related pathways (hypoxia, inflammatory, and fatty acid metabolism) score for myeloid subtypes embedding the UMAP plot. Furthermore, the Ucell algorithm in the “Ucell” package was used to evaluate the enrichment scores of the AS-related pathways in different groups/states, as pseudo-trajectory time progressed.

SCENIC Analysis

The Single-Cell Regulatory Network Inference and Clustering (SCENIC) method was employed to define regulons, which are complexes consisting of a transcription factor (TF), its potential targets, and their activity levels. Macrophage subtypes in different states confirmed using Monocle were further included in SCENIC, and they were then sorted based on states. This process begins with a count matrix representing gene abundances across all macrophage subtypes and unfolds in three steps. Initially, co-expression modules are deduced using “GRNBoost2”, a regression approach targeting individual genes. Further, “cisTarget” is utilized to remove indirect targets from these modules by identifying cis-regulatory motifs. The final step involves quantifying regulon activity by calculating an enrichment score for the target genes of the regulon (AUCell). We employed the “pySCENIC” package, a Python-based tool of the SCENIC pipeline, to analyze TF activity. This study utilized its command-line version, databases for cis-target (hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather), and TF motifs (motifs-v10nr_clust-nr.hgnc-m0.001-o0.0), along with specific command-line parameters. All 1892 TFs with available motifs in the database were included as input. Lastly, the “CalcRSS” function was employed to pinpoint regulons with high Regulon Specificity Scores (RSS), facilitating the identification of macrophage subtypes/states-specific regulons based on RSS ranking, and regulon modules are identified according to connection specificity index (CSI).25 Further details of the SCENIC analysis are described in a previous study.26

Metabolism Pathways Analysis at the Single-Cell Resolution

To assess metabolic activity at the single-cell level, “scMetabolism” package, a novel computational framework for single-cell metabolic activity quantification, was used.27 This method calculates a metabolic activity score for each cell across various metabolic pathways using a single-cell matrix and the Vision algorithm. This package is preloaded with KEGG pathways and Reactome entries. The Vision algorithm also calculates the metabolic score. We compared metabolic activities across different groups/states in various pathways to identify significant differences. For this analysis, we employed the “VISION” method using KEGG metabolic gene sets. Visualization of the results was achieved using the DotPlot.metabolism and BoxPlot.metabolism functions.

Intercellular Communication Analysis

The “CellChat” package was used to infer and visualize the intercellular communication networks between macrophage subtypes and other immune cells at the single-cell levels. This method can predict major signaling inputs and outputs for macrophage subtypes and how those ligand–receptor interactions and specific signaling pathways for functions using network analysis and pattern recognition approaches. The “CellChat” was conducted as described previously,28 and the “circlize” package was used for visualization.

Deconvolution of Cell Type and Gene Expression in AS

The BayesPrism model was employed to construct a matrix from the scRNA-seq of AS as a reference to deconvolute the immune cell abundances of macrophage subtypes in the bulk-seq datasets of AS plaque.29 The scRNA-seq reference were obtained from GSE159677-155512. Default parameters were used for the deconvolution analysis. The BayesPrism model showed the best overall performance and robustness compared to eight other deconvolution methods when deconvolving granular immune lineages, as demonstrated in previously research.30 Spearman correlation of macrophage subtypes abundances was evaluated using the “corrplot” package. Kaplan-Meier survival curves were employed to analyze the differences in ischemic event-free survival between groups with high and low macrophage subtypes abundances in AS plaques. Visualization was performed using the “ggstatsplot” and “ggplot2” packages.

To Explore the Diagnostic Value of Biomarkers in AS by Multiple Machine Learning Algorithms

We used Venn diagram to intersect marker genes identified through scRNA-seq datasets and DEMGs in the bulk-seq. Thus, we obtained biomarkers highly associated with the macrophage subtype in the cluster 3 plaque subtype. The diagnostic value of biomarkers in AS was determined via integrated analysis using three machine learning algorithms: LASSO, RF, and SVM-RFE. Moreover, ROC curve analysis using the “pROC” package to assessed the diagnostic value of biomarkers in multiple GEO datasets. Subsequently, the diagnosis biomarkers were used to construct a diagnostic model in AS, and the model’s performance was evaluated by 10 machine learning algorithms using the “mlr3verse” package. Additionally, the reliability of diagnostic model was validated using the GSE21545-41571 and GSE163154 datasets. The expression of biomarkers in different groups and their correlation with 28 types of immune cells were calculated and visualized using the “ggpubr” and “circlize” packages. Kaplan–Meier survival analysis curves for high versus low biomarker gene expression groups in AS plaques were generated using the “survival” package.

Mendelian Randomization Analysis in AS

Additionally, for biomarkers exhibiting an area under curve (AUC) < 0.6, we executed a two-sample Mendelian Randomization (MR) analysis using the “TwoSampleMR” and “coloc” packages. This analysis was aimed to investigate whether the genetic determinants (eQTLs) associated with these low-AUC biomarkers play a causal role in the advancement of AS. eQTLs linked to AS were identified using the publicly accessible genome-wide association study (GWAS) datasets from the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/). We extracted eQTLs related to the biomarker as the exposure factor from public databases.31 Significant SNPs were selected as instrumental variables by setting a p-value < 5×10−8 and performing population stratification. F-statistics were calculated for the selected eQTLs to assess the strength of instrumental variables. Only eQTLs with F-statistics greater than 10 were retained for further analysis.32 In case of multiple instruments, the inverse variance weighting (IVW) method was used, which uses information on all instruments. Additionally, assessments for reverse causality and colocalization analysis were conducted to determine the pleiotropy between FLT3LG’s eQTL and AS shared the same genetic variation.

Mice Atherosclerotic Experiments

Twelve male ApoE−/− mice, aged 20 weeks, were obtained from Wuhan Servicebio Technology Co., Ltd., for animal experiments. After acclimating for four weeks with a standard chow diet (SCD) (SCD obtained from Jiangsu Xietong Pharmaceutical Bio-engineering Co., Ltd., Nanjing, China) and water, six ApoE−/− mice were switched from SCD to an atherogenic high-cholesterol diet (HCD) (20% kcal protein, 40% kcal fat, and 40% kcal carbohydrate; Beijing Hfk Bioscience Co. Ltd., Beijing, China) for a 16-week long feeding period. Six additional mice, serving as a negative control group, were fed SCD for 16 weeks. In the final stage of the study, all mice were sacrificed by intraperitoneal injection of sodium pentobarbital (60 mg/kg) and euthanized to obtain aortic samples.

Quantitative Real-Time Polymerase Chain Reaction and Western Blot

Biomarker expression was detected by qPCR to validate our bioinformatics results. Total RNA extraction and qPCR were conducted as previously described.33 Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control gene. The relative expression level of the target gene was calculated using the 2−ΔΔCt method. The qPCR forward and reverse primer sequences are as follows: (1) GAPDH, CCTCGTCCCGTAGACAAAATG, and TGAGGTCAATGAAGGGGTCGT; (2) CD36, GGAACTGTGGGCTCATTGCT, and CAACTTCCCTTTTGATTGTCTTCTC; (3) FCER1G, CCGTGATCTTGTTCTTGCTCCT, and GGCTCGGAGAGAATTAGAAGTGG; (4) VAMP8, CTGCCTTGGGTGGAAACAGA, and TTGTTTCGGAGGTGGTCCAG. Western blotting (WB) was performed as previously described17 with primary antibodies against HSP70 (GB12241, Servicebio, China), CD36 (GB112562, Servicebio, China), and FCER1G (BS-13167R, Bioss, China). Western blots were quantified by densitometric analysis with ImageJ software V1.53a.

Multiplex Immunofluorescence Staining

Tissue slices of mouse arterial tissues were obtained from samples set in paraffin. A multiplex immunofluorescence staining assay was conducted as previously described.17 Primary antibodies included anti-CD68 (1:3000, GB113109, Servicebio, China), anti-FCER1G (1:200, BS-13167R, Bioss, China). Slides were counterstained with DAPI (1:5, G1012, Servicebio, China) at room temperature for 10 min. Multispectral images were acquired and performed using a fluorescence microscope (Nikon Eclipse C1, Nikon, Japan) and ImageJ software V1.53a.

Statistical Analysis

The bioinformatics analysis was conducted with R software (version 4.2.1). Survival differences between groups were assessed using Kaplan-Meier curves and Log rank tests. p-value < 0.05 was considered statistically significant. Correlation analysis was performed by Pearson and Spearman correlation. P-value < 0.05 was deemed statistically significant in all cases. Statistical analyses and graphs were performed using GraphPad Prism (version 8 GraphPad Software, La Jolla, CA, USA).

Results Identification and Analysis of High Infiltration Macrophage-Score DEGs in AS

We combined the gene expression matrices GSE21545 and GSE41571 after removing batch effects and standardization, resulting in a total of 229 AS samples with prognostic information related to ischemic events. Supplementary Figure 1 demonstrates the absence of apparent batch effects between the two pairs of datasets after batch removal. Applying a deconvolution method based on ssGSEA, we evaluated the infiltration scores of 28 immune cell types in 229 AS samples (Figure 1A). AS samples were then categorized into high and low macrophage score groups based on the median value of macrophage infiltration scores. As anticipated, survival analysis showed significantly lower ischemic event-free survival in the high macrophage score group compared to the low group (Figure 1B, p-value < 0.05). Further, we identified 1438 up-DEGs in the high macrophage score group (Figure 1C). Utilizing the hallmark gene set from MsigDB, we conducted Spearman correlation analysis to assess the relationship between high macrophage infiltration scores and various hallmark gene sets. A significant positive correlation was found between high macrophage infiltration scores and multiple immune signaling pathways, including inflammation response, pyroptosis, and autophagy (Figure 1D and E). Conversely, a negative correlation was noted with epithelial-mesenchymal transition (EMT) and cholesterol homeostasis. Previous studies have shown that as inflammatory macrophages are continuously activated, the microenvironment within atherosclerotic plaques experiences increasingly intense and complex inflammatory responses.34–36 Combine our results above suggest a complex interplay of immune-inflammatory reactions associated with macrophage infiltration. To clarify the immune-related functions of up-DEGs, we intersected them with 585 IRGs identified in a prior study.19 This process yielded 130 distinct differentially expressed immune-related genes (DEIRGs) (Figure 1F).

Figure 1 Identification of 130 DEIRGs in high macrophage infiltration AS samples. (A) Infiltration scores of 28 immune cell types in 229 AS samples. (B) Kaplan–Meier analysis showed the differences between AS samples with high and low macrophage infiltration scores. (C) Volcano plot of mRNA expression between high and low macrophage infiltration score groups in AS. (D) Scatter plot of correlation between macrophage infiltration score and immune-related pathways with the linear model. (E) Heatmap of the correlations between macrophage infiltration score and five pathways. (F) Venn plot of 130 DEIRGs based on up-DEGs from high macrophage infiltration score groups and IGRs.

Distinguishing Characteristics and Functional Attributes Indicative of a Poorer Prognosis in AS Plaques Subtypes

To further unravel potential mechanisms and subtype characteristics of AS plaques, we applied unsupervised clustering to 130 previously identified DEIRGs. This analysis revealed distinct clustering of 131 carotid atherosclerotic plaques into three clusters: cluster 1, cluster 2, and cluster 3. Intriguingly, DEIRGs demonstrated significantly increased expression in the cluster 3 plaque subtype. Notably, most plaques linked to ischemic events were categorized within this cluster. Furthermore, DEIRGs showed reduced expression in cluster 2, whereas cluster 1 exhibited expression levels intermediate to both clusters (Figure 2A). This finding has not been observed in prior similar studies.16,17 Survival analysis demonstrated significant prognostic differences among the three subtypes, with cluster 3 plaque subtype exhibiting markedly lower ischemic event-free survival (p-value < 0.01) (Figure 2B). PCA results further indicated distinct separation among the three clusters, suggesting diverse features and functions among three plaque subtypes (Figure 2C). Through ssGSEA of the immune landscape, we observed increased infiltration of 28 immune cell types in the cluster 3 plaque subtype. In contrast, the lowest infiltration was noted in cluster 2, with cluster 1 exhibiting intermediate levels (Figure 2D). This variation in immune cell infiltration corresponded with the DEIRGs expression patterns across three AS plaque subtypes, suggesting a potential link between the cluster 3 subtype and more prevalent immune inflammatory responses. Moreover, by employing hallmark gene sets, we evaluated the biological functions of the different plaque subtypes. Consistently, the cluster 3 subtype showed increased expression in several immune-inflammatory signatures, such as apoptosis, autophagy, ferroptosis, pyroptosis, complement activation, fatty acid metabolism, hypoxia, and inflammation (p-value < 0.05). Conversely, no significant differences or lower expression were observed in signatures related to cholesterol homeostasis and EMT (Figure 2E). Given the evident biological functional differences between cluster 3 and the other two clusters, we designated cluster 3 as the Differentially Expressed Macrophage-Related Genes (DEMGs) cluster, while clusters1 and 2 were defined as non-DEMGs clusters. GSVA analysis revealed that the DEMGs cluster was implicated in various immune cell-related signaling pathways, while non-DEMGs clusters were associated with EMT-related signaling pathways, such as vascular smooth muscle contraction and cell adhesion aggregation (Figure 2F and Supplementary Figure 2). Additionally, we assessed the expression of immune checkpoint-related genes in the DEMGs and Non-DEMGs clusters, and as anticipated, multiple immune checkpoint genes were upregulated in the DEMGs cluster, suggesting a more intense immune cell response and inflammatory gene expression in this cluster (Figure 2G).

Figure 2 The consensus clustering of 130 DEIGRs indicated prognostic differences among AS plaque subtypes. (A) The mRNA expression pattern of 130 DEIGRs were observed across the three AS plaque subtypes. (B) Kaplan‐Meier analysis revealed poorer prognosis in the cluster 3 plaque subtype compared to clusters 1 and 2 in AS. (C) PCA analyses showed that the three clusters were fine isolated. (D) Three AS plaque subtypes have different Immune cell infiltration pattern. (E) Three AS plaque subtypes had different biofunction through using the ssGSEA algorithm of well-known signatures. (F) Top 10 KEGG pathways ranked by GSVA score differences between non-DEMGs cluster vs DEMGs cluster. (G) The differences of expression patterns of immune checkpoint-related genes between non-DEMGs cluster vs DEMGs cluster. *p < 0.05; **p < 0.01; ***p < 0.001; and ****p < 0.0001. “ns” indicates not significant (p > 0.05).

Identification of DEMGs in AS Plaque Subtypes Through WGCNA

The WGCNA was employed to determine the association of the 1438 up-DEGs with the cluster 3 plaque subtype. Illustrated in Figure 3A, we set the optimal soft threshold for WGCNA at 8. The black and brown gene modules showed a positive correlation with the cluster 3 plaque subtype, with correlation coefficients of 0.28 and 0.35, respectively (Figure 3B and C, both p-values < 0.05). Moreover, the gene significance of the black and brown modules within cluster 3 exhibited correlation coefficients of 0.6 and 0.37, with p-values < 0.01 (Figure 3D). These results suggest a strong positive correlation between the genes in the black and brown modules and the cluster 3 plaque subtype. Gene modules within the black and brown categories of the WGCNA, which intersect with the previously identified 130 DEIRGs, were defined as DEMGs (Figure 3E). This intersection yielded 43 DEMGs that exhibited a high correlation with the cluster 3 plaque subtype. Consequently, these genes were designated as pivotal in the differential macrophage expression associated with the cluster 3 plaque subtype. GO enrichment analysis of these 43 DEMGs indicated their involvement in positive immune regulation functions, such as the positive regulation of inflammatory effector process, positive regulation of macrophage cytokine production, regulation of platelet activation, and regulation of inflammatory responses (Figure 3F). KEGG analysis revealed that these DEMGs are implicated in immune-inflammatory and lipid metabolism-related signaling pathways, including platelet activation, digestion and absorption of fats, cholesterol metabolism, lipid and atherosclerosis, PI3K-Akt signaling pathway, NF-κB signaling pathway, inflammatory factor interactions, and cell autophagy (Figure 3G). Both GO and KEGG enrichment analyses collectively suggest the participation of these DEMGs in the formation of cluster 3 subtype plaques and their involvement in intricate immune-inflammatory activities within the plaque.

Figure 3 Identification of DEGs correlated with the cluster plaque subtype based on WGCNA. (A) Selection of the soft threshold powers when the scale-free fit index reaches 0.8 and the mean connectivity for various soft threshold powers. (B) Establishment of co-expressed gene modules based on a hierarchical clustering algorithm. (C) Correlations between the gene modules and three AS plaque subtypes. The number at the top of each square indicates the correlation coefficients and the p-values shown in the brackets below (red indicates a positive correlation, blue indicates a negative correlation, and the depth of the color indicates the degree of correlation). (D) Scatter plot of correlation of MM with GS in black and brown modules, respectively. (E) Venn plot of 43 DEMGs based on WGCNA genes and 130 DEIRGs. (F) Circle plots of GO annotation (BP, CC, and MF) of 43 DEMGs in AS. (G) KEGG pathways enriched by 43 DEMGs in AS.

A Single-Cell Transcriptomic Atlas of Paired Human Carotid Artery and Carotid Atherosclerotic Plaque

To provide a more comprehensive understanding of the immune cell composition and potential intercellular connectivity patterns within the cluster 3 plaque subtype, we integrated data from a total of nine single-cell samples sourced from two datasets, GSE159677 and GSE155512. These samples included six carotid atherosclerotic plaques and three adjacent normal vascular tissues. Following data processing and quality control procedures, a total of 54,701 single cells were subjected to subsequent analysis. The major cell-type annotation was performed based on the “SingleR” package and classical gene markers to identify cell types. Ultimately, we identified eight major cell types, including B cells, myeloid cells, endothelial cells, fibroblasts, smooth muscle cells, T cells, NK cells, and mast cells (Figure 4A). The major cell-type atlas showed a significant increase in myeloid cells within the AS group (Figure 4B), and the quality control results for the nine single-cell samples indicated no apparent heterogeneity between samples. The mRNA quantity and mitochondrial ratio in each sample were controlled within a reasonable range (Figure 4C). Additionally, we observed significant changes in cell distribution in the AS compared to the normal group: myeloid cells and T cells increased markedly, whereas endothelial cells and fibroblasts showed a notable decrease in distribution proportions (Figure 4D and E). To assess the impact of gene changes in major cell types between the AS and normal groups, hierarchical clustering was performed based on the expression changes of major cell type genes relative to the normal group. The results clearly indicated an aggregation of these eight major cell types in the AS group (Figure 4F). Moreover, there were similar transcriptional features between the AS and normal groups, further confirming the accuracy of our cell-type annotation. Furthermore, using the 43 DEMGs, we calculated enrichment scores for eight major cell types employing the “AddModuleScore” function. Notably, the DEMGs score exhibited a significant elevation in myeloid cells (Figure 4G). This finding aligning with previous studies,36,37 suggests that the polarization of subtype cells in myeloid cells may be a crucial factor contributing to plaque instability.

Figure 4 Single-cell profile between AS and normal groups. (A) Heatmap plot visualizes the classical markers for each cell type. (B) UMAP plot visualizing eight major cell types identified from scRNA-seq in AS vs normal groups. (C) Violin plot of quality for nine samples in scRNA-seq data. (D) Bar graph showing the percentages of eight major cell types in AS vs normal groups. (E) The percentage of eight major cell types according to their origins from AS tissue. (F) The hierarchical clustering heatmap shows the correlation between eight major cell types in AS vs normal groups. (G) Violin plot shows the DEMGs score for eight major cell types.

Identification of Macrophages Subtypes Related to DEMGs Score in AS at a Single Cell Resolution

We further isolated 11,584 myeloid cells for secondary clustering and annotation. Through literature searches and the identification of classical gene markers, we successfully annotated a total of ten myeloid subtypes. These include classical monocytes, DC1 cells, DC2 cells, DC3 cells, SPP1+/VCAN+ macrophages (Mac-C1-SPP1/VCAN), SPP1+ macrophages (Mac-C2-SPP1), IL1B+ macrophages (Mac-C3-IL1B), TREM2+ macrophages (Mac-C4-TREM2), FLT3LG+ macrophages (Mac-C5-FLT3LG), and monocyte-like macrophages (Monolike Mac). As illustrated in Figure 5A, we employed multiple classical gene markers to annotate these ten myeloid subtypes. Notably, the classical gene markers for SPP1+/VCAN+ macrophages and SPP1+ macrophages exhibited significantly higher expression in the AS group, suggesting these two macrophages may be AS-induced. Previous literature has reported the unique functions of SPP1+ macrophages, including slowed cholesterol clearance, enhanced vascularization, and glycosylation effects.16,38,39 Interestingly, the gene expression profile of Mac-C1-SPP1/VCAN not only include the classical gene markers of SPP1+ macrophage but also includes inflammatory, angiogenesis and lipid generation-related genes, such as EREG, VCAN, and AQP9. These results implied that Mac-C1-SPP1/VCAN may represent a novel macrophage subtype with more intense inflammatory and angiogenesis reactions based on SPP1+ macrophage subtype. Additionally, the single-cell atlas indicates that Mac-C1-SPP1/VCAN is closer in distance to Mac-C3-IL1B and Mac-C5-FLT3LG, suggesting potential functional similarities and closer cell connections among these three subtypes (Figure 5B). Concurrently, as depicted in Supplementary Figure 3, the distribution of Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG gradually increases in AS. However, the distribution patterns of these subtypes in nine single-cell samples suggest a decrease in these three macrophage subtypes in AS. We attribute this result to these three macrophage subtypes being relatively rare cell types in AS4, AS5, and AS6 samples, with an uneven proportion of each AS sample compared to normal groups. The higher number of AS samples may have statistically diluted the proportion of these rare macrophage subtypes (Figure 5C). Based on the previously identified 43 DEMGs, we conducted a subsequent enrichment scoring for these ten myeloid subtypes (Figure 5D). Interestingly, we observed a significant upregulation of DEMGs scores in Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG, suggesting that these three macrophage subtypes may contribute to the increased instability of plaques. Single-cell pathway enrichment analysis was conducted on fifty hallmark gene sets, focusing on the ten identified subtypes of myeloid cells. The results indicate that Mac-C1-SPP1/VCAN is significantly upregulated in angiogenesis, while oxidative phosphorylation and oxygen responses are significantly downregulated. Mac-C3-IL1B is significantly upregulated in angiogenesis, NF-κB signaling pathway, inflammatory response, and IL6-JAK-STAT3 signaling pathway. Mac-C5-FLT3LG is significantly upregulated in NF-κB signaling pathway and apoptosis (Figure 5E). These findings suggest that these three macrophage subtypes may primarily exert biological functions such as angiogenesis, regulation of immune inflammatory responses, and apoptosis. We also conducted enrichment scoring for three AS-related pathways embedded in the UMAP plot: fatty acid metabolism, inflammation, and hypoxia. Interestingly, the high-enrichment regions of these three pathways correspond one-to-one with the regions of the three macrophage subtypes mentioned above (Figure 5F).

Figure 5 Characterization of macrophage subtypes in AS. (A) Bubble plot shows the average expression of the classical markers for 10 myeloid subtypes in AS vs normal groups. (B) UMAP plot visualizing 10 myeloid subtypes identified from scRNA-seq. (C) The percentages of 10 myeloid subtypes across nine samples. (D) Violin plot shows the DEMGs score for 10 myeloid subtypes. (E) Enrichment scores of fifty hallmark gene sets for 10 myeloid subtypes in AS vs normal groups. (F) UMAP plot showing the enrichment regions of three AS-related pathways using the AUCell algorithm.

Evolutionary Trajectory Analysis of Macrophages Subtypes in AS

Due to the proximity observed on the UMAP plot among three subtypes of macrophages, namely Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG. We conducted further pseudo-trajectory analysis to investigate the dynamic immune states and cellular differentiation trajectories of these macrophage subtypes. Simultaneously, we utilized two algorithms, “CytoTRACE” and “Monocle2”, to mutually validate and infer the developmental initiation stages of the macrophages. Impressively, two distinct activation paths emerged, suggesting a three-state differentiation process for the macrophage subtypes, namely state1, state2, and state3 (Figure 6A). Both algorithms indicated that state2 represents the initiation stage of development, characterized by generally high stemness scores and the earliest development onset. State2 then differentiates into state1 and state3 via node1. We defined these three states inferred from the trajectory distribution as Pre-branch, Cell fate 1, and Cell fate 2. Notably, significant differences in the distribution of macrophage subtypes were observed under different activation paths. Specifically, Pre-branch was mainly dominated by Monolike Mac and Mac-C4-TREM2, Cell fate 1 was mainly dominated by Mac-C1-SPP1/VCAN and Mac-C2-SPP1, and Cell fate 2 was mainly dominated by Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG (Figure 6A). Ridgeline plot indicates that among the macrophage subtypes, those in the Pre-branches show the earliest activation in development time. In contrast, the subtype under Cell Fate 1 exhibits the latest activation. The subtype (Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG) associated with Cell Fate 2, on the other hand, displays activation at a mid-to-late stage in the development timeline (Figure 6B). Consistent with expectations, pseudo-trajectory analysis of the AS-related pathways revealed that as pseudo-trajectory time progressed, AS group showed increasing expression scores in four pathways: angiogenesis, fatty acid metabolism, cholesterol homeostasis, and hypoxia. In contrast, the normal group exhibited the opposite trend. Interestingly, Cell fate 2 also demonstrated the highest expression scores in these four pathways, showing a similar expression trend to the AS group (Figure 6C). This finding suggests a key role of these three macrophage subtypes within Cell fate 2 in AS formation. Furthermore, we assessed the enrich trends of DEMGs scores across the three stages. Surprisingly, the DEMGs scores gradually increased along the activation trajectory of Cell fate groups and were significantly higher in Cell fate 2 than in Cell fate 1/Pre-branch (p < 0.0001) (Figure 6D). Subsequently, we investigated the transcriptional changes associated with the differentiation state of DEMGs and observed that, except for DPP4, PTGS1, LAPTM5, and TMC6, the gene expression trends of other DEMGs gradually upregulated with pseudo-trajectory time change, mirroring the expression trend of Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG macrophage subtypes (Figure 6E). Additionally, our SCENIC analysis of macrophage subtypes revealed that their regulons can be categorized into five distinct modules, as shown in Figure 6F. Each specific macrophage subtype exhibited unique TF regulons. Notably, TFs such as E2F3/6/7, ATF2/4, JUNB, IRF1/7, FOS, and ELF1known for their involvement in cell cycle regulation, cellular stress responses, IFN signaling pathways, and immune modulation, were found to be highly activated in Cell fate 2 (Figure 6G). This suggests a nuanced regulatory landscape in these macrophage subtypes, emphasizing their functional diversity in cellular processes. Collectively, tracing gene fluctuations along the bifurcated trajectory, Cell fate 2 appears to play a key role in pathways closely associated with AS, such as angiogenesis, lipid metabolism, and hypoxia. Additionally, DEMGs score exhibit the highest in Cell fate 2, suggesting selective high expression in Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG macrophage subtypes. These macrophage subtypes in Cell fate 2 may be key cell types leading to plaque rupture in AS, resulting in decreased ischemic event-free survival.

Figure 6 Activation trajectory of macrophage subtypes in AS. (A) The development trajectory of macrophage subtypes is colored by cell subtypes, cell states, pseudotime orders, and CytoTRACE scores. (B) Ridgeline plot shows the distribution of the pseudotime-ordered macrophage subtypes among three states. (C) Two-dimensional plots show the expression pattern for signature genes related to hypoxia, angiogenesis, cholesterol homeostasis, and fatty acid metabolism in different groups/states, along with the pseudotime. (D) Right: DEMGs score is enriched along the evolutionary trajectory of macrophage subtypes. Left: The box plot shows the comparison of DEMGs scores among three states. *p < 0.05; **p < 0.01; ***p < 0.001; and ****p < 0.0001. (E) Heatmap shows the expression pattern of DEMGs, gradually upregulated in Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG, along with the pseudotime. (F) Macrophage subtypes can be classified into five modules according to the cooperation of different transcription factors. (G) Heatmap shows the activity of TF regulons among different macrophage subtypes/states.

Identification of the Function and Prognostic Significance of Cell fate 2 Group in AS

Previously, we found significant association between Cell fate 2 and DEMGs scores. Thus, we utilized BEAM analysis to identify the DEGs associated with Cell fate 2. The GO Term results reveal that DEGs with Cell fate 2 are predominantly expressed in pathways related to immune-inflammatory responses, including immune cell receptor signaling pathways, positive regulation of the immune response, and the NF-κB signaling pathway (Figure 7A). Concurrently, Cell fate 1’s DEGs were found to be linked to cell recognition, cell signal transduction, and ATP synthesis. Further survival analysis demonstrated that high enrichment scores of Cell fate 2’s DEG could predict the occurrence of ischemic events in the AS patient cohort, suggesting that these macrophages within the Cell fate 2 can markedly reduce the prognosis of cardiovascular disease patients (Figure 7B). We further investigated the tissue distribution of these macrophages within the Cell fate 2 using a large pair of samples. Intriguingly, in two independent datasets of AAS, a higher enrichment score of Cell fate 2 was observed in AAS (Figure 7C). Additionally, through bulk-seq analysis, significant expression differences were observed among the 41 DEGs of Cell fate 2 across different AS plaque subtypes. Specifically, in the cluster 3 plaque subtype, Cell fate 2’s DEGs exhibited high expression, whereas they exhibited low expression in the cluster 2 subtype and intermediate expression in the cluster 1 subtype (Figure 7D). This finding suggests that these macrophages within Cell fate 2 are likely crucial components of the cluster 3 plaque subtype. Subsequently, we conducted a metabolic pathway analysis for these three states of macrophage subtypes. The results indicated significant upregulation of fatty acid synthesis, fatty acid elongation, and fatty acid degradation pathways in the Cell fate 2 (Figure 7E and F). Furthermore, GO enrichment analysis revealed that Mac-C1-SPP1/VCAN has relevant functions in both inflammation and lipid metabolism (Figure 7E and G). We implied that the inflammatory response of these macrophages within the Cell fate 2 group is one of the factors contributing to the overactivation of lipid metabolism. Existing research indicates a correlation between excessive fatty acid elongation and the development of AS, where fatty acids in the body can transform into pro-inflammatory metabolites, triggering inflammation and exacerbating vascular wall damage.6 In summary, these results support the likelihood that macrophage subtypes within Cell fate 2 accumulate in lesions and play a crucial role in the progression and prognosis of atherosclerosis. Additionally, these results suggest a close association between cluster 3 plaque subtype and the macrophage subtypes within Cell fate 2.

Figure 7 Identification of DEGs with Cell fate 2 and its prognostic significance. (A) Heatmap exhibits the expression dynamics of DEGs with differentiated macrophage states. (B) Kaplan‐Meier analysis shows that high enrichment scores of Cell fate 2’s DEGs lead a poorer prognosis in AS. (C) Box plot shows the difference between high and low-enrichment scores of Cell fate 2’s DEGs in the GSE28829-43292 dataset and GSE163154 dataset. ****p < 0.0001. (D) The mRNA expression patterns of Cell fate 2’s DEGs in three AS plaque subtypes. (E) Heatmap shows the differences in metabolic pathways among different macrophage subtypes/states. (F) Box plot shows the differences in three fatty-metabolic pathways among different macrophage states. (G) Top 10 GO-BP terms of Mac-C1-SPP1/VCAN in AS.

Transcriptome Characterization of Subtype Macrophages in the AS Plaque Subtype

To investigate macrophage subtype proportion changes in AS plaque subtypes, we employed the BayesPrism model trained using referenced literature on scRNA-seq of AS. Subsequently, we inferred proportional cell type fractions and gene expression levels in each cell type from the bulk-seq samples of AS plaque. Figure 8A depicts the proportions of various cell types (at the single-cell level) in 131 bulk-seq samples from carotid artery plaques. A heatmap of intercellular correlations reveals a significant positive correlation (p-value < 0.05) between Mac-C1-SPP1/VCAN and Mac-C3-IL1B, in line with expectations (Figure 8B). Notably, a positive correlation is observed between Mac-C1-SPP1/VCAN and NK cells and mast cells, suggesting potential immunoregulatory and cytotoxic functions, consistent with prior findings (Figure 7G). Interestingly, survival analysis based on the proportion of Mac-C1-SPP1/VCAN indicates a marked reduction in the ischemic event-free survival for high cell proportion (p-value < 0.0001) (Figure 8C). Furthermore, the clustering heatmap of AS plaque subtypes reveals that Mac-C1-SPP1/VCAN, Mac-C3-IL1B, and Mac-C5-FLT3LG exhibit high infiltration proportions in cluster 3 plaques, low infiltration proportions in cluster 2 plaques, and moderate infiltration proportions in cluster 1 plaques, consistent with our earlier findings (Figures 8D and 7D). Subgroup statistics based on AS plaque subtypes for these three macrophage subtypes indicate that, except for Mac-C5-FLT3LG cells showing no statistically significant difference, Mac-C1-SPP1/VCAN and Mac-C3-IL1B cell proportions are highest in cluster 3 plaques (p-value < 0.0001) (Figure 8E). In summary, we have identified the transcriptional features of these three macrophage subtypes in cluster 3 plaques and discovered that the high infiltration proportions of two macrophage subtypes (Mac-C1-SPP1/VCAN and Mac-C3-IL1B) are a crucial factor contributing to plaque instability. Finally, we selected 17 marker genes by intersecting DEMGs from bulk-seq analysis and DEGs from these three macrophage subtypes in scRNA-seq analysis for subsequent validation (Figure 8F).

Figure 8 The cell type percentage from bulk-RNA seq of AS plaque subtypes deconvoluted by BayesPrism. (A) Bar plots show the percentage of each cell type in the bulk-seq of AS plaque deconvolution by BayesPrism. (B) Correlation of the infiltration of pairwise 17 cell types in the bulk-seq of AS plaque deconvoluted by BayesPrism. *p < 0.05 was considered statistically significant. (C) The Kaplan–Meier analysis shows overall survival analyses for low and high percentages of Mac-C1-SPP1/VCAN. (D) Three AS plaque subtypes have different macrophage subtype infiltration patterns. (E) Violin plot shows the comparison of