In this study, we utilize the single-cell database GSE159677 (26) [comprising carotid atherosclerotic plaque (AS) and normal carotid tissue plaque area (PA)] for analysis of copy number variation and annotation of cell types. To accurately annotate cell types, we visualized the results of dimensionality reduction and clustering, as demonstrated in Figures 1A,B. There were 8 distinct cell types identified, including endothelial cells characterized by the expression of VWF (27) and PECAM1 (28), VSMCs distinguished by expression of CALD1 (29)and TAGLN (26, 30), NKT cells labeled with NKG7 (31) and CTSW (17) markers, macrophages defined by C1QA (21) and C1QB (32), CD4⁺ T cells identified by IL7R (33) and LTB (34) expression, plasma cells expressing CD27 (35) and SDC1 (36), mast cells labeled with TPSAB1 (37) and TPSB2 (38) markers, and B cells labeled by MS4A1 (39) and CD79A (40). Using Seurat and ggplot2 to draw the cell population proportional stack map of PA and AS cells tSNE distribution and cell population proportion map, we found that the macrophage cluster was obviously enriched in AS (Figures 1C–E). Gene Ontology (GO) functional analysis based on the differentially expressed genes between AS and PA revealed that this macrophage population was associated with macrophage activation, antigen presentation (MHC II), myeloid cell migration, and other pathways (Figure 1F). Cell-cell interactions in AS and normal tissue samples showed that macrophages have a strong interaction with a variety of cell types (Figure 1G).

To further classify the aforementioned macrophage populations and identify the macrophage subset with the most significant group differences, we generated multiple volcano plots to visualize the clusters of macrophage subsets (M18, M7, M6, M4, M12) (Figure 2A). Using t-SNE distribution and analysis of cell population proportions in PA and AS cells, we identified a significant enrichment of heterogeneous SPP1⁺ FABP5⁺ macrophages within the AS population, which has the largest proportion (Figures 2B–D). To gain deeper insights, we utilized the microarray dataset GSE100927 (41), which includes AS and PA samples, to extract the top 25 highly variable characteristic genes of SPP1⁺ FABP5⁺ macrophages as the enrichment background gene set. Single-sample gene set enrichment analysis (ssGSEA) showed the enrichment score of SPP1⁺ FABP5⁺ macrophages in the AS group (Figure 2E). To clarify the correlation between this macrophage subpopulation and the progression of AS disease, we conducted a functional analysis using GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. The results showed that SPP1⁺ FABP5⁺ macrophages are characterized by lipid transport and storage, cell migration, and chemotaxis, all of which are associated with the pathogenesis of AS (Figures 2F,G). Furthermore, to identify cell populations significantly associated with SPP1⁺ FABP5⁺ macrophages, ssGSEA scores for each cell population were calculated from the microarray data, followed by correlation analysis using the corrplot and ggplot2 packages (Figures 2H–N). Interestingly, scatter plot analysis revealed the highest negative correlation between SPP1⁺ FABP5⁺ macrophages and VSMCs (Figure 2N). This finding contradicts our previous understanding that VSMC-macrophage interactions contribute to the progression of AS. We speculated that this negative correlation might be attributed to the inclusion of the PA cell population in our analysis, where VSMCs exhibit functions that are contrary to those of the SPP1⁺ FABP5⁺ macrophage population. Alternatively, VSMCs show varying functions at different stages of AS. The proliferation of VSMCs can stabilize AS during the early stage of AS (42, 43), which could also explain the negative correlation between VSMC and macrophages that we analyzed. Therefore, we then focus on VSMCs and identify pathogenic subpopulations that contribute to the progression of AS.

As shown in Figure 3A, differential gene expression analysis was performed in the volcano map to identify the top 5 up-regulated markers in each VSMC subset (V17, V16, V11, V10, V5, V3). Notably, we identified a highly specific subset of ITLN1⁺ VSMCs (V3) that expressed traditional VSMC markers while also exhibiting macrophage-associated markers, such as C1QA (21) and C1QB (44). This suggests that these ITLN1⁺ VSMCs represent a distinct foam cell population derived from VSMCs, strongly associated with AS progression. To further investigate this hypothesis, we analyzed the expression of foam cell markers, including CD36, TREM2, APOC1, CTSB, FABP5, and PLIN2 (45–49), across distinct VSMC subpopulations (Figure 3B). Remarkably, ITLN1⁺ VSMCs exhibited robust expression of foam cell markers. Furthermore, the up-regulation of ITLN1 was found to enhance the expression of foam cell-related and lipid metabolism-related genes, such as CD36, OLR1, SRA1, ABCG1, and PPARG (26, 50–53) (Figures 3C,D). Additional analysis revealed differential gene expression between ITLN1_high and ITLN1_low VSMCs. Using the Enhanced Volcano package, we found that PLIN2, a key marker of foam cells, was significantly enriched in ITLN1_high VSMCs (Figure 3E). The pie chart in Figure 3F further underscores the strong correlation between ITLN1⁺ VSMCs and foam cell-related genes. In conclusion, our findings suggest that the ITLN1⁺ VSMC subset represents a foam cell population derived from VSMCs, highlighting its critical role in AS progression.

The tSNE distribution and cell population proportion analyses are performed to distinguish the differences between ITLN1⁺ foam cells and other VSMC subpopulations. The results revealed distinct cell subsets, including CFD⁺ fibromyocytes, C11orf96⁺ VSMCs, and ITLN1⁺ foam cells, with clear differentiation in the AS group compared to the PA group, as shown in Figures 4A–C. For each cell subset, 25 highly variable characteristic genes were selected as background gene sets for enrichment analysis. Notably, in Figures 4D–I, the gene signature of ITLN1⁺ foam cell displayed significant enrichment in AS compared to PA, highlighting a pronounced difference. Further functional analysis of ITLN1⁺ foam cells, as shown in Figures 4J–K, revealed their close association with macrophage activation, polarization, migration, and chemotaxis through GO and KEGG pathway enrichment analyses, implicating their pivotal role in the pathogenesis of AS. Additionally, as shown in Figure 4L, correlation analysis of cell populations was conducted using the same approach as described above. The resulting correlation pie chart revealed a strong relationship (r = 0.90) between ITLN1⁺ foam cells and SPP1⁺ FABP5⁺ macrophages, suggesting significant cross-talk interactions between these two cell clusters in AS.

To elucidate the interactions between these two cell clusters, further analysis was conducted. Based on the preprocessed t-SNE object, cell-cell interactions in AS and normal tissue samples were analyzed using the CellChat package (http://www.cellchat.org/), with the CellChatDB human. Intercellular ligand-receptor pairs were systematically explored in samples from both PA and AS tissues. A relatively dense network of ligand-receptor pairs was observed between ITLN1⁺ foam cells and SPP1⁺ FABP5⁺ macrophages, providing compelling evidence for their close intercellular communication, as shown in Figures 5A,B. Furthermore, Figures 5C,D visualizes the cell-cell communication network between SPP1⁺ FABP5⁺ macrophages and various cell populations derived from VSMCs, as mentioned above. SPP1⁺ FABP5₊ macrophages exhibited a broader spectrum of ligandreceptor interactions with ITLN1⁺ foam cells, among which the MIF-(CD74 + CD44) and SPP1-CD44 ligand-receptor pairs were identified as key mediators facilitating the crosstalk between these distinct cell subsets. Additionally, as shown in Figures 5E–G, the CellChat analysis revealed the presence of ligand receptors, including VEGFB-VEGFR1, NAMPT-INSR, and VEGFA-VEGFR1 (54–56) which play a crucial role in AS, all participate in communication between these two types of cells, underscoring the importance of their interaction in AS.

To emphasize the role of ITLN1⁺ foam cells as a critical endpoint in the differentiation of VSMCs and their pivotal contribution to AS progression, we utilized scRNA-seq data to conduct pseudo-time trajectory analysis, which allowed us to examine the evolution of VSMC cells and the dynamic changes in pathogenic genes expression during differentiation, As the pseudo-time trajectory extended further backward, it signifies a more heterogeneous cellular state. By employing the Slingshot R package, we traced the evolutionary dynamics of single-cell transcriptomes. As shown in Figure 6A, the extension of the evolutionary curve reflects the progressive differentiation of VSMCs. The analysis considered six VSMC subpopulations as potential differentiation trajectories, and we selected lipid metabolism and AS-related genes, including APOC1, APOEA, FABP5, PLOD2, SERPINE1, SPP1, PLIN2, and FABP4 (57–63), to analyze the dynamic changes in their expression levels along the evolutionary trajectory. Dimensionality reduction was performed using the reducedDims function and differentially expressed genes (DEGs) along pseudotime were identified with the tradeSeq package. Marker gene expression trends were visualized using pseudotemporal heatmaps generated with the plotGeneCount function. As shown in Figure 6B, ITLN1⁺ foam cells were found to represent a highly evolved pathogenic VSMC subpopulation, characterized by elevated expression of genes linked to lipid synthesis and AS progression. Among the genes we selected, the expression levels of SPP1 and FABP5 are significantly higher in the evolutionary node of ITLN1⁺ foam cells compared to other pseudo-temporal evolutionary nodes, providing support for our hypothesis regarding the interaction between SPP1⁺ FABP5⁺ macrophages and ITLN1⁺ foam cells. These findings underscore the highly differentiated state of ITLN1⁺ foam cells and their significant role in lipid accumulation and plaque formation, highlighting their critical involvement in AS pathology.

The spatial co-localization of ITLN1⁺ foam cells and SPP1⁺ FABP5⁺ macrophages was analyzed using spatial transcriptomics (ST) datasets GSE243179 (64) (Figure 7A) and GSE159677 (26) (Figure 7B), derived from tissue sections of AS patients. These datasets provided high-resolution localization of specific gene expressions. We analyzed the gene-spot matrix and normalized the data, followed by visualizing the cell clusters on the slices. The analysis revealed a high degree of spatial similarity among ITLN1⁺ foam cells, SPP1⁺ FABP5⁺ macrophages, and SPP1⁺ APOC1⁺ macrophages. This finding was unexpected, as the SPP1⁺ APOC1⁺ macrophage subpopulation was previously identified as significantly smaller than SPP1⁺ FABP5⁺ macrophages. Moreover, the correlation heatmap (Figure 4L) demonstrated that SPP1⁺ APOC1⁺ macrophages had a lower correlation with ITLN1⁺ foam cells (r = 0.83) compared to SPP1⁺ FABP5⁺ macrophages (r = 0.90). As a result, SPP1⁺ APOC1⁺ macrophages were not emphasized in the analysis. Nevertheless, the significant correlation suggests that SPP1 may play a critical role as a macrophage-related gene. Collectively, the high expression of SPP1 likely enables macrophage subpopulations to interact with ITLN1⁺ foam cells, potentially contributing to the progression of AS (Figure 8).

The information regarding the dataset is illustrated in tabular form, as depicted in Table 1.

In order to identify different cell types, the Seurat package was used to remove low-quality cells, including those with more than 8,000 detected genes or fewer than 200, as well as those with mitochondrial gene expression exceeding 20% of total expression. After filtration, data normalization was performed using the NormalizeData function within the Seurat package, with a scaling factor of 10,000. Next, the FindVariableFeatures function was used to identify highly variable genes, followed by RunPCA for dimensionality reduction on the normalized dataset. Clustering analysis was then performed using the FindNeighbors and FindClusters functions. Dimensionality reduction and clustering results were visualized using UMAP and t-SNE techniques. Cell types were annotated using the classical marker genes of known cell types in combination with the SingleR packages. Marker genes for each cluster were defined based on their expression levels, where these genes were expressed in more than 30% of cells and had an average expression value exceeding 0.5.

To perform enrichment and analysis of distinct cell populations, the Seurat package was used to preprocess the scRNA-seq data and conduct cell clustering. The metadata of cells was extracted from the processed Seurat object, including the sample ID, group, cell types, and cell clusters. Using the dplyr package, cells were grouped, and the number of different cell populations in each sample was calculated. Subsequently, the proportions of each cell population were determined. The results of the cell composition proportions for each group were visualized using the ggplot2 package.

The dataset GSE100927 containing AS and PA data was used to identify the top 25 most variable feature genes from each cell subpopulation, thereby creating an enriched background gene set. The GSVA package was employed to perform single-sample gene set enrichment analysis (ssGSEA) on this dataset and calculate the ssGSEA score for the dataset. After calculating the score, the enrichment scores of these clusters in the AS and PA groups were compared to measure the difference in the proportion of specific cell populations in these two groups.

To clarify the biological roles of cell subpopulations, we conducted analyses using gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG), establishing a significance threshold of P < 0.05 to identify relevant pathways.

To assess the activity of various ligand-receptor pairs across distinct cell populations and the enrichment of signaling pathways, the CellChat analysis was conducted. Based on the preprocessed tSNE object, cell-cell interactions in AS and normal tissue samples were analyzed by running the CellChat package (http://www.cellchat.org/), with the interactive ligand-receptor database set as CellChatDB human.

We performed a pseudo-temporal trajectory analysis on scRNA-seq datasets using Slingshot R package to investigate the dynamic changes in gene expression in different cell states. Slingshot enables the inference of developmental trajectories by identifying continuous, ordered cell states, and ordering cells along a pseudotime axis that represents their progression through a biological process. Using the reducedDims function in Slingshot, we apply methods such as PCA, t-SNE to project the cells into a lower-dimensional space, facilitating the identification of distinct cell groups. These projections were then processed by Slingsho to infer the pseudotime trajectory. To identify the differentially expressed genes (DEGs) along the pseudotime continuum, we use the tradeSeq package, which is designed for modeling gene expression along continuous trajectories. To better understand how gene expression changes over pseudotime, we use the plotGeneCount function to create heatmaps visualizing the expression of marker genes or DEGs across the pseudotime continuum.

In order to demonstrate the spatial interactions between cell subpopulations, the spatial transcriptomics (ST) datasets GSE243179 and GSE159677 from two AS patients were incorporated. The Seurat package was used to process both ST and Visium sample data, generating gene-spot matrices. Gene expression data were normalized using the SCTransform function. Dimensionality reduction and clustering were performed using principal component analysis (PCA), followed by mapping cell clusters onto tissue images with the SpatialDimPlot function. Signature scoring was conducted using the AddMetaData function in the Seurat package, based on scRNA-seq or ST signatures, with default parameters. Finally, spatial feature expression plots were generated using the SpatialFeaturePlot function in the Seurat package.